AtomsBase.jl is a common interface for representing atomic structures in Julia. DFTK directly supports using such structures to run a calculation as is demonstrated here.
using DFTK
using AtomsBuilder
In this example we construct a bulk silicon system using the bulk function
from AtomsBuilder. This function
uses tabulated data to set up a reasonable starting geometry and lattice for bulk silicon.
system = bulk(:Si)
FlexibleSystem(Si₂, periodicity = TTT):
cell_vectors : [ 0 2.715 2.715;
2.715 0 2.715;
2.715 2.715 0]u"Å"
Atom(Si, [ 0, 0, 0]u"Å")
Atom(Si, [ 1.3575, 1.3575, 1.3575]u"Å")
By default the atoms of an AbstractSystem employ the bare Coulomb potential.
To employ pseudpotential models (which is almost always advisable for
plane-wave DFT) one employs the pseudopotential keyword argument in
model constructors such as model_DFT.
For example we can employ a PseudoFamily object
from the PseudoPotentialData
package. See its documentation for more information on the available
pseudopotential families and how to select them.
using PseudoPotentialData # defines PseudoFamily
pd_lda_family = PseudoFamily("dojo.nc.sr.lda.v0_4_1.standard.upf")
model = model_DFT(system; functionals=LDA(), temperature=1e-3,
pseudopotentials=pd_lda_family)
Model(lda_x+lda_c_pw, 3D):
lattice (in Bohr) : [0 , 5.13061 , 5.13061 ]
[5.13061 , 0 , 5.13061 ]
[5.13061 , 5.13061 , 0 ]
unit cell volume : 270.11 Bohr³
atoms : Si₂
atom potentials : ElementPsp(Si, "/home/runner/.julia/artifacts/326db5c901e2681584ec5c06fc17f6c96e516ff9/Si.upf")
ElementPsp(Si, "/home/runner/.julia/artifacts/326db5c901e2681584ec5c06fc17f6c96e516ff9/Si.upf")
num. electrons : 8
spin polarization : none
temperature : 0.001 Ha
smearing : DFTK.Smearing.FermiDirac()
terms : Kinetic()
AtomicLocal()
AtomicNonlocal()
Ewald(nothing)
PspCorrection()
Hartree()
Xc(lda_x, lda_c_pw)
Entropy()
Alternatively the pseudopotentials object also accepts a Dict{Symbol,String},
which provides for each element symbol the filename or identifier
of the pseudopotential to be employed, e.g.
path_to_pspfile = PseudoFamily("cp2k.nc.sr.lda.v0_1.semicore.gth")[:Si]
model = model_DFT(system; functionals=LDA(), temperature=1e-3,
pseudopotentials=Dict(:Si => path_to_pspfile))
Model(lda_x+lda_c_pw, 3D):
lattice (in Bohr) : [0 , 5.13061 , 5.13061 ]
[5.13061 , 0 , 5.13061 ]
[5.13061 , 5.13061 , 0 ]
unit cell volume : 270.11 Bohr³
atoms : Si₂
atom potentials : ElementPsp(Si, "/home/runner/.julia/artifacts/966fd9cdcd7dbaba6dc2bf43ee50dd81e63e8837/Si.gth")
ElementPsp(Si, "/home/runner/.julia/artifacts/966fd9cdcd7dbaba6dc2bf43ee50dd81e63e8837/Si.gth")
num. electrons : 8
spin polarization : none
temperature : 0.001 Ha
smearing : DFTK.Smearing.FermiDirac()
terms : Kinetic()
AtomicLocal()
AtomicNonlocal()
Ewald(nothing)
PspCorrection()
Hartree()
Xc(lda_x, lda_c_pw)
Entropy()
We can then discretise such a model and solve:
basis = PlaneWaveBasis(model; Ecut=15, kgrid=[4, 4, 4])
scfres = self_consistent_field(basis, tol=1e-8);
n Energy log10(ΔE) log10(Δρ) Diag Δtime --- --------------- --------- --------- ---- ------ 1 -7.921681503863 -0.69 5.5 839ms 2 -7.926133575373 -2.35 -1.22 1.0 770ms 3 -7.926833524972 -3.15 -2.37 2.0 152ms 4 -7.926861274900 -4.56 -2.99 2.9 199ms 5 -7.926861645407 -6.43 -3.38 2.2 162ms 6 -7.926861670500 -7.60 -3.84 1.4 141ms 7 -7.926861679389 -8.05 -4.17 1.8 146ms 8 -7.926861681752 -8.63 -4.90 1.6 148ms 9 -7.926861681834 -10.08 -5.00 2.2 160ms 10 -7.926861681869 -10.46 -5.69 1.1 137ms 11 -7.926861681872 -11.48 -6.12 2.2 162ms 12 -7.926861681873 -12.59 -7.10 1.8 152ms 13 -7.926861681873 -13.52 -7.41 2.9 202ms 14 -7.926861681873 -14.75 -7.75 1.1 141ms 15 -7.926861681873 + -Inf -8.16 1.6 196ms
If we did not want to use AtomsBuilder we could of course use any other package which yields an AbstractSystem object. This includes:
using AtomsIO
system = load_system("Si.extxyz");
Run the LDA calculation:
pseudopotentials = PseudoFamily("cp2k.nc.sr.lda.v0_1.semicore.gth")
model = model_DFT(system; pseudopotentials, functionals=LDA(), temperature=1e-3)
basis = PlaneWaveBasis(model; Ecut=15, kgrid=[4, 4, 4])
scfres = self_consistent_field(basis, tol=1e-8);
n Energy log10(ΔE) log10(Δρ) Diag Δtime --- --------------- --------- --------- ---- ------ 1 -7.921730371012 -0.69 5.8 187ms 2 -7.926131411206 -2.36 -1.22 1.0 133ms 3 -7.926833365786 -3.15 -2.37 2.1 154ms 4 -7.926861259195 -4.55 -3.00 3.1 241ms 5 -7.926861647203 -6.41 -3.39 2.1 159ms 6 -7.926861671123 -7.62 -3.84 1.4 139ms 7 -7.926861679640 -8.07 -4.18 1.8 172ms 8 -7.926861681772 -8.67 -4.90 1.4 140ms 9 -7.926861681856 -10.08 -5.18 1.8 150ms 10 -7.926861681871 -10.82 -5.84 1.8 148ms 11 -7.926861681873 -11.74 -6.41 2.0 155ms 12 -7.926861681873 -13.00 -7.52 1.8 155ms 13 -7.926861681873 -14.15 -8.00 3.5 197ms 14 -7.926861681873 + -14.75 -8.78 2.0 153ms
The same could be achieved using ExtXYZ
by system = Atoms(read_frame("Si.extxyz")),
since the ExtXYZ.Atoms object is directly AtomsBase-compatible.
using AtomsBase
using Unitful
using UnitfulAtomic
# Construct a system in the AtomsBase world
a = 10.26u"bohr" # Silicon lattice constant
lattice = a / 2 * [[0, 1, 1.], # Lattice as vector of vectors
[1, 0, 1.],
[1, 1, 0.]]
atoms = [:Si => ones(3)/8, :Si => -ones(3)/8]
system = periodic_system(atoms, lattice; fractional=true)
# Now run the LDA calculation:
pseudopotentials = PseudoFamily("cp2k.nc.sr.lda.v0_1.semicore.gth")
model = model_DFT(system; pseudopotentials, functionals=LDA(), temperature=1e-3)
basis = PlaneWaveBasis(model; Ecut=15, kgrid=[4, 4, 4])
scfres = self_consistent_field(basis, tol=1e-4);
n Energy log10(ΔE) log10(Δρ) Diag Δtime --- --------------- --------- --------- ---- ------ 1 -7.921728831917 -0.69 5.6 218ms 2 -7.926138025114 -2.36 -1.22 1.0 143ms 3 -7.926836807033 -3.16 -2.37 2.0 147ms 4 -7.926864499287 -4.56 -2.99 3.0 195ms 5 -7.926865060810 -6.25 -3.37 2.1 156ms 6 -7.926865082027 -7.67 -3.84 1.2 128ms 7 -7.926865089458 -8.13 -4.06 1.9 141ms
At any point we can also get back the DFTK model as an
AtomsBase-compatible AbstractSystem:
second_system = atomic_system(model)
FlexibleSystem(Si₂, periodicity = TTT):
cell_vectors : [ 0 5.13 5.13;
5.13 0 5.13;
5.13 5.13 0]u"a₀"
Atom(Si, [ 1.2825, 1.2825, 1.2825]u"a₀")
Atom(Si, [ -1.2825, -1.2825, -1.2825]u"a₀")
Similarly DFTK offers a method to the atomic_system and periodic_system functions
(from AtomsBase), which enable a seamless conversion of the usual data structures for
setting up DFTK calculations into an AbstractSystem:
lattice = 5.431u"Å" / 2 * [[0 1 1.];
[1 0 1.];
[1 1 0.]];
Si = ElementPsp(:Si, pseudopotentials)
atoms = [Si, Si]
positions = [ones(3)/8, -ones(3)/8]
third_system = atomic_system(lattice, atoms, positions)
FlexibleSystem(Si₂, periodicity = TTT):
cell_vectors : [ 0 5.13155 5.13155;
5.13155 0 5.13155;
5.13155 5.13155 0]u"a₀"
Atom(Si, [ 1.28289, 1.28289, 1.28289]u"a₀")
Atom(Si, [-1.28289, -1.28289, -1.28289]u"a₀")