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"Å")
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 --- --------------- --------- --------- ---- 1 -7.921682118187 -0.69 5.9 2 -7.926161830555 -2.35 -1.22 1.0 3 -7.926837794412 -3.17 -2.37 1.8 4 -7.926861520074 -4.62 -3.04 2.9 5 -7.926861644079 -6.91 -3.41 1.9 6 -7.926861669404 -7.60 -3.78 1.6 7 -7.926861680431 -7.96 -4.33 1.2 8 -7.926861681752 -8.88 -4.81 2.0 9 -7.926861681859 -9.97 -5.22 1.5 10 -7.926861681871 -10.91 -5.90 1.6 11 -7.926861681873 -11.86 -7.07 2.1 12 -7.926861681873 -13.47 -7.42 3.5 13 -7.926861681873 -15.05 -8.29 1.4
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 --- --------------- --------- --------- ---- 1 -7.921699031998 -0.69 5.6 2 -7.926167595907 -2.35 -1.22 1.0 3 -7.926837017823 -3.17 -2.37 1.9 4 -7.926861520391 -4.61 -3.03 2.6 5 -7.926861642562 -6.91 -3.39 1.9 6 -7.926861669757 -7.57 -3.79 1.9 7 -7.926861680378 -7.97 -4.33 1.2 8 -7.926861681766 -8.86 -4.86 2.0 9 -7.926861681862 -10.02 -5.29 1.8 10 -7.926861681871 -11.01 -5.86 2.0 11 -7.926861681873 -11.92 -6.99 1.8 12 -7.926861681873 -13.48 -7.67 3.2 13 -7.926861681873 -15.05 -8.20 2.9
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 --- --------------- --------- --------- ---- 1 -7.921684822302 -0.69 5.6 2 -7.926171465626 -2.35 -1.22 1.0 3 -7.926840765557 -3.17 -2.37 1.8 4 -7.926864901222 -4.62 -3.04 2.6 5 -7.926865055580 -6.81 -3.41 2.2 6 -7.926865081823 -7.58 -3.81 1.9 7 -7.926865091371 -8.02 -4.29 1.2
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₀")
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₀")