AtomsBase integration
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 AtomsBuilderFeeding an AtomsBase AbstractSystem to DFTK
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₂
pseudopot. family : PseudoFamily("dojo.nc.sr.lda.v0_4_1.standard.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.921736482861 -0.69 5.6 253ms
2 -7.926133562027 -2.36 -1.22 1.0 218ms
3 -7.926834409913 -3.15 -2.37 2.0 171ms
4 -7.926861244988 -4.57 -3.03 3.0 217ms
5 -7.926861658425 -6.38 -3.50 1.8 173ms
6 -7.926861676882 -7.73 -4.10 1.9 166ms
7 -7.926861678287 -8.85 -4.09 2.0 175ms
8 -7.926861677756 + -9.28 -3.99 1.0 147ms
9 -7.926861680553 -8.55 -3.56 1.0 157ms
10 -7.926861681675 -8.95 -4.00 1.0 150ms
11 -7.926861681802 -9.90 -4.16 1.0 159ms
12 -7.926861681802 + -13.06 -4.17 1.0 156ms
13 -7.926861681822 -10.69 -4.26 1.0 149ms
14 -7.926861681864 -10.38 -4.88 1.0 155ms
15 -7.926861681864 -13.30 -4.84 1.0 150ms
16 -7.926861681866 -11.66 -5.12 1.0 158ms
17 -7.926861681859 + -11.11 -4.73 1.0 150ms
18 -7.926861681872 -10.88 -5.69 1.0 156ms
19 -7.926861681872 -12.26 -5.80 1.4 162ms
20 -7.926861681873 -12.73 -5.91 1.0 150ms
21 -7.926861681873 -13.23 -6.19 1.0 162ms
22 -7.926861681873 -14.35 -6.41 1.0 150ms
23 -7.926861681873 -14.75 -6.52 1.0 155ms
24 -7.926861681873 + -14.15 -6.50 1.0 149ms
25 -7.926861681873 -14.45 -6.65 1.0 157ms
26 -7.926861681873 -14.21 -7.02 1.0 151ms
27 -7.926861681873 -14.45 -7.36 1.0 159ms
28 -7.926861681873 + -14.75 -7.25 1.2 160ms
29 -7.926861681873 + -15.05 -8.31 1.0 150msIf we did not want to use AtomsBuilder we could of course use any other package which yields an AbstractSystem object. This includes:
Reading a system using AtomsIO
Read a file using AtomsIO, which directly yields an AbstractSystem.
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.921712523041 -0.69 5.5 209ms
2 -7.926133346057 -2.35 -1.22 1.0 175ms
3 -7.926833249452 -3.15 -2.37 2.0 168ms
4 -7.926861287193 -4.55 -2.99 2.9 212ms
5 -7.926861647510 -6.44 -3.39 2.0 178ms
6 -7.926861670692 -7.63 -3.85 1.6 160ms
7 -7.926861679625 -8.05 -4.20 1.4 160ms
8 -7.926861681785 -8.67 -5.12 1.8 161ms
9 -7.926861681857 -10.14 -5.16 2.6 198ms
10 -7.926861681868 -10.94 -5.53 1.0 156ms
11 -7.926861681873 -11.35 -6.43 1.0 150ms
12 -7.926861681873 -13.13 -6.99 2.6 189ms
13 -7.926861681873 -14.05 -7.31 2.0 177ms
14 -7.926861681873 + -15.05 -7.52 1.1 161ms
15 -7.926861681873 + -15.05 -8.77 1.0 159msThe same could be achieved using ExtXYZ by system = Atoms(read_frame("Si.extxyz")), since the ExtXYZ.Atoms object is directly AtomsBase-compatible.
Directly setting up a system in AtomsBase
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.921701522490 -0.69 5.5 305ms
2 -7.926133862367 -2.35 -1.22 1.0 155ms
3 -7.926836177660 -3.15 -2.37 2.0 209ms
4 -7.926864701205 -4.54 -2.98 3.0 242ms
5 -7.926865050391 -6.46 -3.34 1.9 195ms
6 -7.926865079422 -7.54 -3.76 1.4 161ms
7 -7.926865091139 -7.93 -4.26 1.2 163msObtaining an AbstractSystem from DFTK data
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₀")