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
In this example we construct a silicon system using the ase.build.bulk
routine
from the atomistic simulation environment
(ASE), which is exposed by ASEconvert
as an AtomsBase AbstractSystem
.
# Construct bulk system and convert to an AbstractSystem
using ASEconvert
system_ase = ase.build.bulk("Si")
system = pyconvert(AbstractSystem, system_ase)
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
To use an AbstractSystem in DFTK, we attach pseudopotentials, construct a DFT model, discretise and solve:
system = attach_psp(system; Si="hgh/lda/si-q4")
model = model_LDA(system; 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.921707077037 -0.69 5.6 2 -7.926168133060 -2.35 -1.22 1.0 252ms 3 -7.926837334351 -3.17 -2.37 1.6 235ms 4 -7.926861518265 -4.62 -3.03 3.0 262ms 5 -7.926861642237 -6.91 -3.38 1.9 248ms 6 -7.926861669582 -7.56 -3.78 1.9 219ms 7 -7.926861680503 -7.96 -4.34 1.2 203ms 8 -7.926861681759 -8.90 -4.85 2.1 228ms 9 -7.926861681859 -10.00 -5.23 1.9 246ms 10 -7.926861681871 -10.91 -5.86 1.8 219ms 11 -7.926861681873 -11.83 -6.98 2.0 234ms 12 -7.926861681873 -13.56 -7.34 3.5 284ms 13 -7.926861681873 -14.35 -8.26 1.2 207ms
If we did not want to use ASE we could of course use any other package which yields an AbstractSystem object. This includes:
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_LDA(system; 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.921703754529 -0.69 5.8 2 -7.926164546407 -2.35 -1.22 1.0 213ms 3 -7.926836478339 -3.17 -2.37 1.6 232ms 4 -7.926861454194 -4.60 -3.00 2.8 293ms 5 -7.926861630589 -6.75 -3.33 1.9 221ms 6 -7.926861665512 -7.46 -3.70 1.8 218ms 7 -7.926861680459 -7.83 -4.34 1.2 208ms 8 -7.926861681789 -8.88 -4.89 2.0 254ms 9 -7.926861681865 -10.12 -5.36 1.9 223ms 10 -7.926861681872 -11.16 -6.00 2.0 227ms 11 -7.926861681873 -12.12 -6.86 2.2 240ms 12 -7.926861681873 -13.59 -7.40 2.5 260ms 13 -7.926861681873 -14.75 -8.11 1.9 221ms
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:
system = attach_psp(system; Si="hgh/lda/si-q4")
model = model_LDA(system; 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.921713359642 -0.69 5.9 2 -7.926165461257 -2.35 -1.22 1.0 205ms 3 -7.926840189044 -3.17 -2.37 1.9 268ms 4 -7.926864919612 -4.61 -3.01 2.6 252ms 5 -7.926865039504 -6.92 -3.33 1.6 204ms 6 -7.926865077386 -7.42 -3.71 1.6 208ms 7 -7.926865091851 -7.84 -4.42 1.0 185ms
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