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.921671305758 -0.69 5.8 2 -7.926164107424 -2.35 -1.22 1.0 250ms 3 -7.926837685413 -3.17 -2.37 1.6 325ms 4 -7.926861531114 -4.62 -3.04 2.9 298ms 5 -7.926861649114 -6.93 -3.43 1.6 250ms 6 -7.926861671688 -7.65 -3.84 1.8 288ms 7 -7.926861680394 -8.06 -4.32 1.2 240ms 8 -7.926861681757 -8.87 -4.87 2.0 258ms 9 -7.926861681862 -9.98 -5.28 2.0 265ms 10 -7.926861681871 -11.04 -5.79 1.8 270ms 11 -7.926861681873 -11.76 -6.43 1.8 258ms 12 -7.926861681873 -13.07 -7.40 2.0 268ms 13 -7.926861681873 -14.75 -8.18 3.1 311ms
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.921719603968 -0.69 5.6 2 -7.926168395527 -2.35 -1.22 1.0 251ms 3 -7.926836444110 -3.18 -2.37 1.9 343ms 4 -7.926861501842 -4.60 -3.01 2.5 290ms 5 -7.926861632765 -6.88 -3.33 1.8 256ms 6 -7.926861666827 -7.47 -3.73 1.8 253ms 7 -7.926861680698 -7.86 -4.40 1.2 278ms 8 -7.926861681828 -8.95 -5.03 2.1 267ms 9 -7.926861681851 -10.64 -5.13 1.9 262ms 10 -7.926861681871 -10.70 -5.87 1.1 237ms 11 -7.926861681873 -11.87 -6.87 2.0 291ms 12 -7.926861681873 -13.28 -7.34 3.0 305ms 13 -7.926861681873 + -14.75 -8.31 2.2 270ms
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.921699595531 -0.69 5.8 2 -7.926169244709 -2.35 -1.22 1.0 256ms 3 -7.926839506380 -3.17 -2.37 1.6 288ms 4 -7.926864878291 -4.60 -3.00 2.8 344ms 5 -7.926865032982 -6.81 -3.29 1.9 284ms 6 -7.926865076868 -7.36 -3.70 1.6 296ms 7 -7.926865091899 -7.82 -4.42 1.4 253ms
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₀")