Creating and modelling metallic supercells

In this section we will be concerned with modelling supercells of aluminium. When dealing with periodic problems there is no unique definition of the lattice: Clearly any duplication of the lattice along an axis is also a valid repetitive unit to describe exactly the same system. This is exactly what a supercell is: An $n$-fold repetition along one (or multiple) axes of the original lattice.

The following code achieves this for aluminium:

using AtomsBuilder
using DFTK
using LinearAlgebra
using Unitful
using UnitfulAtomic
using PseudoPotentialData

function aluminium_setup(repeat=1; Ecut=7.0, kgrid=[2, 2, 2])
    # Use AtomsBuilder to setup aluminium cubic unit cell (4 Al atoms)
    # with provided lattice constant, see [AtomsBase integration](@ref) for details.
    unit_cell = bulk(:Al; a=7.65339u"bohr", cubic=true)
    supercell = unit_cell * (repeat, 1, 1)  # Make a supercell

    # Select standard pseudodojo pseudopotentials, construct an LDA model, discretize
    # Note: We disable symmetries explicitly here. Otherwise the problem sizes
    #       we are able to run on the CI are too simple to observe the numerical
    #       instabilities we want to trigger here.
    pseudopotentials = PseudoFamily("dojo.nc.sr.lda.v0_4_1.standard.upf")
    model = model_DFT(supercell; pseudopotentials, functionals=LDA(),
                      temperature=1e-3, symmetries=false)
    PlaneWaveBasis(model; Ecut, kgrid)
end;

As expected we obtain the unit cell for repeat=1:

aluminium_setup(1)
PlaneWaveBasis discretization:
    architecture         : DFTK.CPU()
    num. mpi processes   : 1
    num. julia threads   : 1
    num. DFTK  threads   : 1
    num. blas  threads   : 2
    num. fft   threads   : 1

    Ecut                 : 7.0 Ha
    fft_size             : (24, 24, 24), 13824 total points
    kgrid                : MonkhorstPack([2, 2, 2])
    num.   red. kpoints  : 8
    num. irred. kpoints  : 8

    Discretized Model(lda_x+lda_c_pw, 3D):
        lattice (in Bohr)    : [7.65339   , 0         , 0         ]
                               [0         , 7.65339   , 0         ]
                               [0         , 0         , 7.65339   ]
        unit cell volume     : 448.29 Bohr³
    
        atoms                : Al₄
        pseudopot. family    : PseudoFamily("dojo.nc.sr.lda.v0_4_1.standard.upf")
    
        num. electrons       : 12
        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()

and 5-fold as large supercell with repeat=5:

aluminium_setup(5)
PlaneWaveBasis discretization:
    architecture         : DFTK.CPU()
    num. mpi processes   : 1
    num. julia threads   : 1
    num. DFTK  threads   : 1
    num. blas  threads   : 2
    num. fft   threads   : 1

    Ecut                 : 7.0 Ha
    fft_size             : (96, 24, 24), 55296 total points
    kgrid                : MonkhorstPack([2, 2, 2])
    num.   red. kpoints  : 8
    num. irred. kpoints  : 8

    Discretized Model(lda_x+lda_c_pw, 3D):
        lattice (in Bohr)    : [38.267    , 0         , 0         ]
                               [0         , 7.65339   , 0         ]
                               [0         , 0         , 7.65339   ]
        unit cell volume     : 2241.5 Bohr³
    
        atoms                : Al₂₀
        pseudopot. family    : PseudoFamily("dojo.nc.sr.lda.v0_4_1.standard.upf")
    
        num. electrons       : 60
        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()

As we will see in this notebook the modelling of a system generally becomes harder if the system becomes larger.

  • This sounds like a trivial statement as per se the cost per SCF step increases as the system (and thus $N$) gets larger.
  • But there is more to it: If one is not careful also the number of SCF iterations increases as the system gets larger.
  • The aim of a proper computational treatment of such supercells is therefore to ensure that the number of SCF iterations remains constant when the system size increases.

For achieving the latter DFTK by default employs the LdosMixing preconditioner [HL2021] during the SCF iterations. This mixing approach is completely parameter free, but still automatically adapts to the treated system in order to efficiently prevent charge sloshing. As a result, modelling aluminium slabs indeed takes roughly the same number of SCF iterations irrespective of the supercell size:

M. F. Herbst and A. Levitt. Black-box inhomogeneous preconditioning for self-consistent field iterations in density functional theory. J. Phys. Cond. Matt 33 085503 (2021). ArXiv:2009.01665

self_consistent_field(aluminium_setup(1); tol=1e-4);
n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime
---   ---------------   ---------   ---------   ----   ------
  1   -9.355337359065                   -1.10    5.8    186ms
  2   -9.356790306988       -2.84       -1.43    1.0   83.6ms
  3   -9.357059455208       -3.57       -2.80    4.0    121ms
  4   -9.357110755747       -4.29       -3.04    6.4    192ms
  5   -9.357111071416       -6.50       -3.17    1.0   77.2ms
  6   -9.357111281526       -6.68       -3.31    1.0   75.0ms
  7   -9.357111409833       -6.89       -3.47    2.2    105ms
  8   -9.357111465200       -7.26       -3.63    1.4   87.7ms
  9   -9.357111497203       -7.49       -3.89    1.0   79.1ms
 10   -9.357111505457       -8.08       -4.14    2.2   92.3ms
self_consistent_field(aluminium_setup(2); tol=1e-4);
n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime
---   ---------------   ---------   ---------   ----   ------
  1   -18.74731195672                   -0.97    6.5    403ms
  2   -18.74990534443       -2.59       -1.36    1.0    205ms
  3   -18.79220120375       -1.37       -1.74    5.8    351ms
┌ Warning: Eigensolver not converged
  n_iter =
   8-element Vector{Int64}:
     1
    11
     1
     1
     1
     1
     1
     4
@ DFTK ~/work/DFTK.jl/DFTK.jl/src/scf/self_consistent_field.jl:76
  4   -18.79246436625       -3.58       -1.76    2.6    284ms
  5   -18.79255758853       -4.03       -1.84    2.0    249ms
  6   -18.79245095600   +   -3.97       -2.30    2.4    249ms
  7   -18.79238683060   +   -4.19       -2.36    2.6    274ms
  8   -18.79250287266       -3.94       -2.53    1.4    234ms
  9   -18.79256074465       -4.24       -2.69    1.0    224ms
 10   -18.79256068865   +   -7.25       -2.75    1.4    216ms
 11   -18.79255892515   +   -5.75       -2.76    1.0    209ms
 12   -18.79257315679       -4.85       -2.86    1.0    200ms
 13   -18.79258685210       -4.86       -2.97    1.0    208ms
 14   -18.79259041543       -5.45       -3.01    1.0    197ms
 15   -18.79259768654       -5.14       -3.14    1.0    208ms
 16   -18.79260229750       -5.34       -3.29    3.0    249ms
 17   -18.79260629352       -5.40       -3.89    1.9    252ms
 18   -18.79260571160   +   -6.24       -3.73    3.5    326ms
 19   -18.79260348534   +   -5.65       -3.39    1.9    248ms
 20   -18.79260638498       -5.54       -4.52    2.5    286ms
self_consistent_field(aluminium_setup(4); tol=1e-4);
n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime
---   ---------------   ---------   ---------   ----   ------
  1   -37.55052136156                   -0.84    9.9    1.62s
  2   -37.55827994778       -2.11       -1.22    2.2    811ms
┌ Warning: Eigensolver not converged
  n_iter =
   8-element Vector{Int64}:
    14
     9
    15
    18
    24
    25
    15
    12
@ DFTK ~/work/DFTK.jl/DFTK.jl/src/scf/self_consistent_field.jl:76
  3   -37.56437284777       -2.22       -2.11   16.5    1.68s
  4   -37.56488045641       -3.29       -2.35    9.6    1.53s
  5   -37.56494469596       -4.19       -3.11    2.1    867ms
  6   -37.56494834906       -5.44       -3.37    8.6    2.97s
  7   -37.56494975204       -5.85       -3.87    1.8    783ms
  8   -37.56494985174       -7.00       -4.23    4.1    1.21s

When switching off explicitly the LdosMixing, by selecting mixing=SimpleMixing(), the performance of number of required SCF steps starts to increase as we increase the size of the modelled problem:

self_consistent_field(aluminium_setup(1); tol=1e-4, mixing=SimpleMixing());
n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime
---   ---------------   ---------   ---------   ----   ------
  1   -9.355188882949                   -1.10    6.2    158ms
  2   -9.356819948545       -2.79       -1.90    1.0   82.9ms
  3   -9.357089285669       -3.57       -2.63    4.2    126ms
  4   -9.357086993871   +   -5.64       -2.56    4.8    147ms
  5   -9.357111043587       -4.62       -3.47    1.0   69.0ms
  6   -9.357111477446       -6.36       -4.07    3.9    117ms
self_consistent_field(aluminium_setup(4); tol=1e-4, mixing=SimpleMixing());
n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime
---   ---------------   ---------   ---------   ----   ------
  1   -37.55050041320                   -0.84    8.9    1.60s
  2   -37.55358401490       -2.51       -1.61    1.4    668ms
┌ Warning: Detected large minimal occupation 0.0020691154172238894. SCF could be unstable. Try switching to adaptive band selection (`nbandsalg=AdaptiveBands(model)`) or request more converged bands than 29 (e.g. `nbandsalg=AdaptiveBands(model; n_bands_converge=32`)
@ DFTK ~/work/DFTK.jl/DFTK.jl/src/scf/self_consistent_field.jl:87
  3   -30.97344236171   +    0.82       -0.66    8.8    1.56s
  4   -37.52889757654        0.82       -1.71    7.2    1.55s
  5   -37.53657744639       -2.11       -1.77    1.1    702ms
  6   -37.53628815314   +   -3.54       -1.77    1.0    696ms
  7   -37.15545737574   +   -0.42       -1.26    4.9    1.09s
  8   -37.52985802489       -0.43       -1.75    5.0    1.14s
  9   -37.55577121387       -1.59       -1.90    2.1    841ms
 10   -37.56295089956       -2.14       -2.27    1.5    734ms
 11   -37.56410269439       -2.94       -2.38    3.2    974ms
 12   -37.56414701107       -4.35       -2.38    1.2    733ms
 13   -37.56401844705   +   -3.89       -2.39    1.0    695ms
 14   -37.56435868546       -3.47       -2.54    1.5    767ms
 15   -37.56463037975       -3.57       -2.52    3.0    945ms
 16   -37.56452862037   +   -3.99       -2.55    1.1    707ms
 17   -37.56482645980       -3.53       -2.65    1.0    700ms
 18   -37.56479263026   +   -4.47       -2.69    1.0    699ms
 19   -37.56487974351       -4.06       -2.82    1.0    697ms
 20   -37.56491900069       -4.41       -2.98    1.0    699ms
 21   -37.56493820401       -4.72       -3.20    2.4    871ms
 22   -37.56494568068       -5.13       -3.59    2.1    870ms
 23   -37.56494947407       -5.42       -4.00    3.2    1.06s
 24   -37.56494968308       -6.68       -4.24    4.0    1.07s

For completion let us note that the more traditional mixing=KerkerMixing() approach would also help in this particular setting to obtain a constant number of SCF iterations for an increasing system size (try it!). In contrast to LdosMixing, however, KerkerMixing is only suitable to model bulk metallic system (like the case we are considering here). When modelling metallic surfaces or mixtures of metals and insulators, KerkerMixing fails, while LdosMixing still works well. See the Modelling a gallium arsenide surface example or [HL2021] for details. Due to the general applicability of LdosMixing this method is the default mixing approach in DFTK.