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 Δtime --- --------------- --------- --------- ---- ------ 1 -7.921708560224 -0.69 6.1 2 -7.926165810615 -2.35 -1.22 1.0 343ms 3 -7.926836917725 -3.17 -2.37 1.8 427ms 4 -7.926861528980 -4.61 -3.01 2.8 386ms 5 -7.926861632232 -6.99 -3.33 2.0 325ms 6 -7.926861667071 -7.46 -3.73 1.4 314ms 7 -7.926861680728 -7.86 -4.43 1.4 339ms 8 -7.926861681831 -8.96 -5.06 2.1 397ms 9 -7.926861681859 -10.55 -5.22 2.1 338ms 10 -7.926861681872 -10.90 -5.97 1.2 305ms 11 -7.926861681873 -12.05 -6.94 2.4 369ms 12 -7.926861681873 -13.73 -7.52 2.8 417ms 13 -7.926861681873 -14.75 -8.31 2.5 351ms
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.921654140179 -0.69 5.8 2 -7.926165781197 -2.35 -1.22 1.0 305ms 3 -7.926835706901 -3.17 -2.37 1.5 326ms 4 -7.926861523784 -4.59 -3.02 2.8 427ms 5 -7.926861643868 -6.92 -3.39 1.8 310ms 6 -7.926861669831 -7.59 -3.80 2.0 339ms 7 -7.926861680205 -7.98 -4.30 1.5 311ms 8 -7.926861681766 -8.81 -4.87 2.2 392ms 9 -7.926861681851 -10.07 -5.15 1.8 329ms 10 -7.926861681870 -10.71 -5.76 1.6 319ms 11 -7.926861681873 -11.67 -6.66 2.0 338ms 12 -7.926861681873 -13.14 -7.15 2.6 385ms 13 -7.926861681873 -14.75 -7.94 1.8 326ms 14 -7.926861681873 + -14.75 -8.70 3.0 371ms
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.921686629323 -0.69 5.9 2 -7.926169278827 -2.35 -1.22 1.0 288ms 3 -7.926840648276 -3.17 -2.37 1.9 334ms 4 -7.926864907343 -4.62 -3.01 2.6 417ms 5 -7.926865045191 -6.86 -3.34 1.8 303ms 6 -7.926865077476 -7.49 -3.73 1.5 286ms 7 -7.926865091641 -7.85 -4.36 1.0 264ms
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₀")