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 AtomsBuilder

Feeding 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₂, periodicity = TTT):
    cell_vectors      : [       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"Å")

By default the atoms of an AbstractSystem employ the bare Coulomb potential. To employ pseudpotential models (which is almost always advisable for plane-wave DFT) one employs the pseudopotential keyword argument in model constructors such as model_DFT. For example we can employ a PseudoFamily object from the PseudoPotentialData package. See its documentation for more information on the available pseudopotential families and how to select them.

using PseudoPotentialData  # defines PseudoFamily

pd_lda_family = PseudoFamily("dojo.nc.sr.lda.v0_4_1.standard.upf")
model = model_DFT(system; functionals=LDA(), temperature=1e-3,
                  pseudopotentials=pd_lda_family)
Model(lda_x+lda_c_pw, 3D):
    lattice (in Bohr)    : [0         , 5.13061   , 5.13061   ]
                           [5.13061   , 0         , 5.13061   ]
                           [5.13061   , 5.13061   , 0         ]
    unit cell volume     : 270.11 Bohr³

    atoms                : Si₂
    pseudopot. family    : PseudoFamily("dojo.nc.sr.lda.v0_4_1.standard.upf")

    num. electrons       : 8
    spin polarization    : none
    temperature          : 0.001 Ha
    smearing             : DFTK.Smearing.FermiDirac()

    terms                : Kinetic()
                           AtomicLocal()
                           AtomicNonlocal()
                           Ewald(nothing)
                           PspCorrection()
                           Hartree()
                           Xc(lda_x, lda_c_pw)
                           Entropy()

Alternatively the pseudopotentials object also accepts a Dict{Symbol,String}, which provides for each element symbol the filename or identifier of the pseudopotential to be employed, e.g.

path_to_pspfile = PseudoFamily("cp2k.nc.sr.lda.v0_1.semicore.gth")[:Si]
model = model_DFT(system; functionals=LDA(), temperature=1e-3,
                  pseudopotentials=Dict(:Si => path_to_pspfile))
Model(lda_x+lda_c_pw, 3D):
    lattice (in Bohr)    : [0         , 5.13061   , 5.13061   ]
                           [5.13061   , 0         , 5.13061   ]
                           [5.13061   , 5.13061   , 0         ]
    unit cell volume     : 270.11 Bohr³

    atoms                : Si₂
    atom potentials      : ElementPsp(:Si, "/home/runner/.julia/artifacts/966fd9cdcd7dbaba6dc2bf43ee50dd81e63e8837/Si.gth")
                           ElementPsp(:Si, "/home/runner/.julia/artifacts/966fd9cdcd7dbaba6dc2bf43ee50dd81e63e8837/Si.gth")

    num. electrons       : 8
    spin polarization    : none
    temperature          : 0.001 Ha
    smearing             : DFTK.Smearing.FermiDirac()

    terms                : Kinetic()
                           AtomicLocal()
                           AtomicNonlocal()
                           Ewald(nothing)
                           PspCorrection()
                           Hartree()
                           Xc(lda_x, lda_c_pw)
                           Entropy()

We can then discretise such a model and solve:

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.921697705881                   -0.69    5.6    253ms
  2   -7.926134119084       -2.35       -1.22    1.0    190ms
  3   -7.926833342142       -3.16       -2.37    2.0    201ms
  4   -7.926861279063       -4.55       -3.02    3.0    210ms
  5   -7.926861655723       -6.42       -3.46    2.1    188ms
  6   -7.926861674732       -7.72       -4.01    1.4    168ms
  7   -7.926861678825       -8.39       -4.12    1.6    166ms
  8   -7.926861681735       -8.54       -4.60    1.0    156ms
  9   -7.926861681829      -10.03       -4.93    1.2    154ms
 10   -7.926861681863      -10.47       -5.43    1.8    174ms
 11   -7.926861681872      -11.05       -5.77    1.4    159ms
 12   -7.926861681873      -12.49       -6.40    1.2    170ms
 13   -7.926861681873      -13.34       -6.29    2.2    183ms
 14   -7.926861681873      -13.88       -6.97    1.0    151ms
 15   -7.926861681873   +  -14.57       -7.69    1.4    164ms
 16   -7.926861681873   +  -14.75       -8.32    2.5    190ms

If we did not want to use AtomsBuilder we could of course use any other package which yields an AbstractSystem object. This includes:

Reading a system using AtomsIO

Read a file using AtomsIO, which directly yields an AbstractSystem.

using AtomsIO
system = load_system("Si.extxyz");

Run the LDA calculation:

pseudopotentials = PseudoFamily("cp2k.nc.sr.lda.v0_1.semicore.gth")
model  = model_DFT(system; pseudopotentials, 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.921709608728                   -0.69    5.6    247ms
  2   -7.926136056682       -2.35       -1.22    1.0    181ms
  3   -7.926834561108       -3.16       -2.37    2.0    173ms
  4   -7.926861292012       -4.57       -3.02    2.9    212ms
  5   -7.926861658115       -6.44       -3.47    2.1    176ms
  6   -7.926861675503       -7.76       -4.05    1.8    172ms
  7   -7.926861677875       -8.62       -4.06    2.0    178ms
  8   -7.926861677098   +   -9.11       -4.03    1.0    149ms
  9   -7.926861681663       -8.34       -4.32    1.0    157ms
 10   -7.926861681810       -9.83       -4.29    1.0    150ms
 11   -7.926861681832      -10.66       -4.34    1.0    157ms
 12   -7.926861681852      -10.70       -4.51    1.0    152ms
 13   -7.926861681857      -11.30       -4.60    1.0    157ms
 14   -7.926861681836   +  -10.68       -4.38    1.0    152ms
 15   -7.926861681842      -11.23       -4.39    1.0    159ms
 16   -7.926861681870      -10.55       -4.91    1.0    160ms
 17   -7.926861681870   +  -12.54       -4.89    1.0    153ms
 18   -7.926861681873      -11.57       -7.16    1.0    159ms
 19   -7.926861681873      -13.32       -7.42    3.8    221ms
 20   -7.926861681873      -15.05       -7.77    1.2    165ms
 21   -7.926861681873      -15.05       -8.09    1.1    160ms

The 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:
pseudopotentials = PseudoFamily("cp2k.nc.sr.lda.v0_1.semicore.gth")
model  = model_DFT(system; pseudopotentials, 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.921702335331                   -0.69    5.5    304ms
  2   -7.926137737936       -2.35       -1.22    1.0    156ms
  3   -7.926837816325       -3.15       -2.37    2.0    211ms
  4   -7.926864689107       -4.57       -3.02    3.1    247ms
  5   -7.926865066216       -6.42       -3.44    1.9    192ms
  6   -7.926865085915       -7.71       -3.97    1.6    172ms
  7   -7.926865089234       -8.48       -4.06    1.8    182ms

Obtaining 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₂, periodicity = TTT):
    cell_vectors      : [       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, pseudopotentials)
atoms     = [Si, Si]
positions = [ones(3)/8, -ones(3)/8]

third_system = atomic_system(lattice, atoms, positions)
FlexibleSystem(Si₂, periodicity = TTT):
    cell_vectors      : [       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₀")