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₂, periodic = TTT):
bounding_box : [ 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"Å")
Si
Si
By default the atoms of an AbstractSystem employ the bare Coulomb potential. To make calculations feasible for plane-wave DFT we thus attach pseudopotential information, before passing the system to construct a DFT model, discretise and solve:
system = attach_psp(system; Si="hgh/lda/si-q4")
model = model_DFT(system; 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.921727585126 -0.69 6.0 252ms
2 -7.926164328327 -2.35 -1.22 1.0 151ms
3 -7.926839424585 -3.17 -2.37 1.9 171ms
4 -7.926861220504 -4.66 -3.05 2.2 187ms
5 -7.926861656653 -6.36 -3.44 1.8 160ms
6 -7.926861672948 -7.79 -3.88 1.8 185ms
7 -7.926861678462 -8.26 -4.06 1.5 157ms
8 -7.926861681780 -8.48 -4.83 1.0 148ms
9 -7.926861681852 -10.15 -5.13 2.6 184ms
10 -7.926861681871 -10.71 -6.08 1.4 172ms
11 -7.926861681873 -11.84 -6.27 2.6 201ms
12 -7.926861681873 + -13.88 -5.99 1.0 166ms
13 -7.926861681873 -13.46 -7.14 1.0 152ms
14 -7.926861681873 -15.05 -7.65 2.1 180ms
15 -7.926861681873 + -14.75 -8.67 2.1 176msIf we did not want to use ASE we could of course use any other package which yields an AbstractSystem object. This includes:
Reading a system using AtomsIO
using AtomsIO
# Read a file using [AtomsIO](https://github.com/mfherbst/AtomsIO.jl),
# which directly yields an AbstractSystem.
system = load_system("Si.extxyz")
# Now run the LDA calculation:
system = attach_psp(system; Si="hgh/lda/si-q4")
model = model_DFT(system; 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.921713900613 -0.69 6.0 257ms
2 -7.926162692150 -2.35 -1.22 1.0 160ms
3 -7.926839042500 -3.17 -2.37 2.0 184ms
4 -7.926861307622 -4.65 -3.04 2.1 193ms
5 -7.926861652433 -6.46 -3.40 2.0 186ms
6 -7.926861670101 -7.75 -3.80 1.4 158ms
7 -7.926861679151 -8.04 -4.12 1.1 162ms
8 -7.926861681796 -8.58 -5.07 1.6 160ms
9 -7.926861681845 -10.32 -5.04 2.8 195ms
10 -7.926861681872 -10.57 -6.16 1.0 148ms
11 -7.926861681873 -11.97 -6.44 2.9 199ms
12 -7.926861681873 -13.80 -6.74 1.2 169ms
13 -7.926861681873 -15.05 -7.48 1.6 170ms
14 -7.926861681873 + -14.45 -7.68 2.8 188ms
15 -7.926861681873 -14.57 -8.95 1.5 161msThe 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:
system = attach_psp(system; Si="hgh/lda/si-q4")
model = model_DFT(system; 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.921727092783 -0.69 5.9 329ms
2 -7.926167528302 -2.35 -1.22 1.0 170ms
3 -7.926842318741 -3.17 -2.37 1.9 219ms
4 -7.926864619908 -4.65 -3.03 2.2 238ms
5 -7.926865056398 -6.36 -3.38 1.8 195ms
6 -7.926865079084 -7.64 -3.75 1.5 212ms
7 -7.926865090198 -7.95 -4.12 1.5 170msObtaining 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₂, periodic = TTT):
bounding_box : [ 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₀")
Si
Si
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; psp=load_psp("hgh/lda/Si-q4"))
atoms = [Si, Si]
positions = [ones(3)/8, -ones(3)/8]
third_system = atomic_system(lattice, atoms, positions)FlexibleSystem(Si₂, periodic = TTT):
bounding_box : [ 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₀")
Si
Si