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 of the axes of the original lattice.
The following code achieves this for aluminium:
using DFTK
using LinearAlgebra
using ASEconvert
function aluminium_setup(repeat=1; Ecut=7.0, kgrid=[2, 2, 2])
a = 7.65339
lattice = a * Matrix(I, 3, 3)
Al = ElementPsp(:Al, psp=load_psp("hgh/lda/al-q3"))
atoms = [Al, Al, Al, Al]
positions = [[0.0, 0.0, 0.0], [0.0, 0.5, 0.5], [0.5, 0.0, 0.5], [0.5, 0.5, 0.0]]
unit_cell = periodic_system(lattice, atoms, positions)
# Make supercell in ASE:
# We convert our lattice to the conventions used in ASE, make the supercell
# and then convert back ...
supercell_ase = convert_ase(unit_cell) * pytuple((repeat, 1, 1))
supercell = pyconvert(AbstractSystem, supercell_ase)
# Unfortunately right now the conversion to ASE drops the pseudopotential information,
# so we need to reattach it:
supercell = attach_psp(supercell, Al="hgh/lda/al-q3")
# Construct an LDA model and discretise
# 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.
model = model_LDA(supercell; temperature=1e-3, symmetries=false)
PlaneWaveBasis(model; Ecut, kgrid)
end;
As part of the code we are using a routine inside the ASE, the atomistic simulation environment for creating the supercell and make use of the two-way interoperability of DFTK and ASE. For more details on this aspect see the documentation on Input and output formats.
Write an example supercell structure to a file to plot it:
setup = aluminium_setup(5)
convert_ase(periodic_system(setup.model)).write("al_supercell.png")
┌ Warning: Skipping atomic property pseudopotential, which is not supported in ASE. └ @ ASEconvert ~/.julia/packages/ASEconvert/neSFl/src/ASEconvert.jl:123 ┌ Warning: Skipping atomic property pseudopotential, which is not supported in ASE. └ @ ASEconvert ~/.julia/packages/ASEconvert/neSFl/src/ASEconvert.jl:123
Python None
As we will see in this notebook the modelling of a system generally becomes harder if the system becomes larger.
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:
[^HL2021]: 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);
┌ Warning: Skipping atomic property pseudopotential, which is not supported in ASE. └ @ ASEconvert ~/.julia/packages/ASEconvert/neSFl/src/ASEconvert.jl:123 n Energy log10(ΔE) log10(Δρ) Diag Δtime --- --------------- --------- --------- ---- ------ 1 -8.298524216091 -0.85 5.5 2 -8.300217858134 -2.77 -1.25 1.0 130ms 3 -8.300438491804 -3.66 -1.88 1.8 154ms 4 -8.300461243153 -4.64 -2.73 2.6 165ms 5 -8.300464407522 -5.50 -3.23 2.6 180ms 6 -8.300464523519 -6.94 -3.36 5.0 217ms 7 -8.300464579853 -7.25 -3.51 1.6 200ms 8 -8.300464614218 -7.46 -3.65 5.2 220ms 9 -8.300464635949 -7.66 -3.84 1.6 161ms 10 -8.300464640061 -8.39 -3.96 2.0 163ms 11 -8.300464643685 -8.44 -4.28 2.0 232ms
self_consistent_field(aluminium_setup(2); tol=1e-4);
┌ Warning: Skipping atomic property pseudopotential, which is not supported in ASE. └ @ ASEconvert ~/.julia/packages/ASEconvert/neSFl/src/ASEconvert.jl:123 n Energy log10(ΔE) log10(Δρ) Diag Δtime --- --------------- --------- --------- ---- ------ 1 -16.67370293141 -0.70 5.9 2 -16.67794951680 -2.37 -1.15 1.0 334ms 3 -16.67916905052 -2.91 -1.87 4.2 432ms 4 -16.67927692484 -3.97 -2.69 2.4 379ms 5 -16.67928540256 -5.07 -3.02 4.5 529ms 6 -16.67928619756 -6.10 -3.46 2.9 365ms 7 -16.67928621800 -7.69 -3.97 2.4 353ms 8 -16.67928622157 -8.45 -4.55 2.9 416ms
self_consistent_field(aluminium_setup(4); tol=1e-4);
┌ Warning: Skipping atomic property pseudopotential, which is not supported in ASE. └ @ ASEconvert ~/.julia/packages/ASEconvert/neSFl/src/ASEconvert.jl:123 n Energy log10(ΔE) log10(Δρ) Diag Δtime --- --------------- --------- --------- ---- ------ 1 -33.32726864166 -0.56 7.9 2 -33.33448338055 -2.14 -1.00 1.5 1.27s 3 -33.33595848559 -2.83 -1.72 5.2 1.69s 4 -33.33611213616 -3.81 -2.62 5.6 1.63s 5 -33.33675603604 -3.19 -2.19 5.0 1.90s 6 -33.33679359508 -4.43 -2.24 1.0 1.16s 7 -33.33694305616 -3.83 -3.39 2.0 1.37s 8 -33.33694374177 -6.16 -3.64 5.8 1.95s 9 -33.33694377401 -7.49 -3.75 1.5 1.15s 10 -33.33694377382 + -9.72 -4.27 1.1 1.11s
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());
┌ Warning: Skipping atomic property pseudopotential, which is not supported in ASE. └ @ ASEconvert ~/.julia/packages/ASEconvert/neSFl/src/ASEconvert.jl:123 n Energy log10(ΔE) log10(Δρ) Diag Δtime --- --------------- --------- --------- ---- ------ 1 -8.298664190119 -0.85 5.1 2 -8.300278883579 -2.79 -1.59 1.2 109ms 3 -8.300375578874 -4.01 -2.30 2.9 141ms 4 -8.300323660104 + -4.28 -2.18 3.6 170ms 5 -8.300464003153 -3.85 -3.47 1.0 116ms 6 -8.300464521231 -6.29 -3.70 7.9 294ms 7 -8.300464632305 -6.95 -4.08 1.9 179ms
self_consistent_field(aluminium_setup(4); tol=1e-4, mixing=SimpleMixing());
┌ Warning: Skipping atomic property pseudopotential, which is not supported in ASE. └ @ ASEconvert ~/.julia/packages/ASEconvert/neSFl/src/ASEconvert.jl:123 n Energy log10(ΔE) log10(Δρ) Diag Δtime --- --------------- --------- --------- ---- ------ 1 -33.32454347069 -0.56 6.5 2 -33.27675355090 + -1.32 -1.25 1.6 1.03s 3 +12.47587617986 + 1.66 -0.23 8.5 2.96s 4 -33.31822551715 1.66 -1.68 6.4 2.86s 5 -33.08167393688 + -0.63 -1.22 4.2 1.89s 6 -30.82774946587 + 0.35 -0.87 3.6 1.69s 7 -33.31784035850 0.40 -1.88 4.2 1.67s 8 -33.33254948581 -1.83 -2.00 2.1 1.26s 9 -33.33486989943 -2.63 -2.12 2.4 1.38s 10 -33.33672375444 -2.73 -2.74 1.9 1.24s 11 -33.33679364530 -4.16 -2.85 3.4 1.47s 12 -33.33688397935 -4.04 -3.03 1.6 1.14s 13 -33.33693892172 -4.26 -3.32 2.8 1.35s 14 -33.33694020020 -5.89 -3.56 2.1 1.21s 15 -33.33694024850 -7.32 -3.73 3.4 1.51s 16 -33.33694349164 -5.49 -4.22 2.8 1.39s
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.