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

Feeding an AtomsBase AbstractSystem to 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.921728328688                   -0.69    5.9    224ms
  2   -7.926165571315       -2.35       -1.22    1.0    157ms
  3   -7.926838914180       -3.17       -2.37    1.9    203ms
  4   -7.926861197652       -4.65       -3.03    2.2    269ms
  5   -7.926861647829       -6.35       -3.39    2.1    174ms
  6   -7.926861668345       -7.69       -3.76    1.6    157ms
  7   -7.926861678631       -7.99       -4.07    1.2    150ms
  8   -7.926861681784       -8.50       -5.00    1.6    158ms
  9   -7.926861681858      -10.13       -5.17    3.0    215ms
 10   -7.926861681872      -10.86       -6.34    1.0    161ms
 11   -7.926861681873      -12.09       -6.49    3.0    205ms
 12   -7.926861681873      -14.21       -7.39    1.0    153ms
┌ Warning: Eigensolver not converged
  n_iter =
   8-element Vector{Int64}:
    3
    2
    2
    2
    2
    3
    2
    3
@ DFTK ~/work/DFTK.jl/DFTK.jl/src/scf/self_consistent_field.jl:76
 13   -7.926861681873      -14.75       -7.46    2.4    217ms
 14   -7.926861681873   +  -15.05       -7.75    1.0    152ms
 15   -7.926861681873   +    -Inf       -9.32    1.1    287ms

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

Reading a system using AtomsIO

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.921731818733                   -0.69    5.9    251ms
  2   -7.926163308194       -2.35       -1.22    1.0    181ms
  3   -7.926838429207       -3.17       -2.37    1.9    221ms
  4   -7.926861210646       -4.64       -3.02    2.1    218ms
  5   -7.926861637625       -6.37       -3.34    1.9    718ms
  6   -7.926861664105       -7.58       -3.68    1.5    151ms
  7   -7.926861679602       -7.81       -4.18    1.1    147ms
  8   -7.926861681807       -8.66       -5.20    1.8    156ms
  9   -7.926861681863      -10.25       -5.22    3.2    190ms
 10   -7.926861681872      -11.02       -6.51    1.0    143ms
 11   -7.926861681873      -12.40       -6.78    3.0    602ms
 12   -7.926861681873      -14.75       -7.22    1.2    150ms
 13   -7.926861681873      -14.75       -7.91    1.2    148ms
┌ Warning: Eigensolver not converged
  n_iter =
   8-element Vector{Int64}:
    2
    2
    3
    2
    2
    3
    2
    3
@ DFTK ~/work/DFTK.jl/DFTK.jl/src/scf/self_consistent_field.jl:76
 14   -7.926861681873   +  -15.05       -8.55    2.4    177ms

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:
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.921727148106                   -0.69    5.9    266ms
  2   -7.926166107149       -2.35       -1.22    1.0    159ms
  3   -7.926842576541       -3.17       -2.37    1.9    182ms
  4   -7.926864635414       -4.66       -3.03    2.1    206ms
  5   -7.926865057874       -6.37       -3.38    1.9    172ms
  6   -7.926865079219       -7.67       -3.75    1.5    155ms
  7   -7.926865090563       -7.95       -4.16    1.1    148ms

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₂, 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