Temperature and metallic systems

In this example we consider the modeling of a magnesium lattice as a simple example for a metallic system. For our treatment we will use the PBE exchange-correlation functional. First we import required packages and setup the lattice. Again notice that DFTK uses the convention that lattice vectors are specified column by column.

using DFTK
using Plots
using PseudoPotentialData
using Unitful
using UnitfulAtomic

a = 3.01794  # Bohr
b = 5.22722  # Bohr
c = 9.77362  # Bohr
lattice = [[-a -a  0]; [-b  b  0]; [0   0 -c]]

pseudopotentials = PseudoFamily("dojo.nc.sr.pbe.v0_4_1.standard.upf")
Mg = ElementPsp(:Mg, pseudopotentials)
atoms     = [Mg, Mg]
positions = [[2/3, 1/3, 1/4], [1/3, 2/3, 3/4]];

Next we build the PBE model and discretize it. Since magnesium is a metal we apply a small smearing temperature to ease convergence using the Fermi-Dirac smearing scheme.

kgrid = KgridSpacing(0.9 / u"Å")    # Minimal spacing of k-points,
#                                      in units of wavevectors (inverse Bohrs)
temperature = 0.01                    # Smearing temperature in Hartree
Ecut = 10                             # Kinetic energy cutoff in Hartree
smearing = DFTK.Smearing.FermiDirac() # Smearing method
#                                      also supported: Gaussian,
#                                      MarzariVanderbilt,
#                                      and MethfesselPaxton(order)

model = model_DFT(lattice, atoms, positions;
                  functionals=[:gga_x_pbe, :gga_c_pbe], temperature, smearing)
basis = PlaneWaveBasis(model; Ecut, kgrid);

Note, that in this example both the Ecut as well as the minimal $k$-point spacing 0.9 / u"Å" far too large to give a converged result. In the online documentation we have used these small values to obtain a fast execution time.

basis_default = PlaneWaveBasis(model)

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                 : 42.0 Ha
    fft_size             : (40, 40, 60), 96000 total points
    kgrid                : MonkhorstPack([9, 9, 5])
    num.   red. kpoints  : 405
    num. irred. kpoints  : 36

    Discretized Model(gga_x_pbe+gga_c_pbe, 3D):
        lattice (in Bohr)    : [-3.01794  , -3.01794  , 0         ]
                               [-5.22722  , 5.22722   , 0         ]
                               [0         , 0         , -9.77362  ]
        unit cell volume     : 308.37 Bohr³
    
        atoms                : Mg₂
        pseudopot. family    : PseudoFamily("dojo.nc.sr.pbe.v0_4_1.standard.upf")
    
        num. electrons       : 20
        spin polarization    : none
        temperature          : 0.01 Ha
        smearing             : DFTK.Smearing.FermiDirac()
    
        terms                : Kinetic()
                               AtomicLocal()
                               AtomicNonlocal()
                               Ewald(nothing)
                               PspCorrection()
                               Hartree()
                               Xc(gga_x_pbe, gga_c_pbe)
                               Entropy()

As can be seen the default discretisation selects the much finer discretisation parameters Ecut=42 and kgrid=[9, 9, 5].

Finally we run the SCF. Two magnesium atoms in our pseudopotential model result in four valence electrons being explicitly treated. Nevertheless this SCF will solve for eight bands by default in order to capture partial occupations beyond the Fermi level due to the employed smearing scheme. In this example we use a damping of 0.8. The default LdosMixing should be suitable to converge metallic systems like the one we model here. For the sake of demonstration we still switch to Kerker mixing here.

scfres = self_consistent_field(basis, damping=0.8, mixing=KerkerMixing());

n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime
---   ---------------   ---------   ---------   ----   ------
  1   -107.0208276793                    0.39    6.0   93.4ms
  2   -107.6655331400       -0.19       -0.47    2.7   90.3ms
  3   -107.6687243799       -2.50       -1.59    3.2   66.2ms
  4   -107.6687503021       -4.59       -2.77    1.0   54.4ms
  5   -107.6687513525       -5.98       -3.43    3.8   70.6ms
  6   -107.6687514223       -7.16       -3.96    3.0   67.7ms
  7   -107.6687514268       -8.35       -4.80    3.2   68.8ms
  8   -107.6687514272       -9.39       -5.06    3.0   68.0ms
  9   -107.6687514272      -10.86       -5.75    1.3   57.1ms
 10   -107.6687514272      -11.96       -6.72    3.0   64.1ms

scfres.occupation[1]

17-element Vector{Float64}:
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0
 1.9999999999999098
 1.9999272442014087
 1.91216798236583
 6.51142956540349e-17
 6.510827455231163e-17
 7.050065578943727e-18
 7.377378846608309e-19
 7.376804344778055e-19
 1.785400767465959e-19

scfres.energies

Energy breakdown (in Ha):
    Kinetic             48.3202481
    AtomicLocal         -101.1884682
    AtomicNonlocal      -19.3973009
    Ewald               -53.8605542
    PspCorrection       1.1973309 
    Hartree             30.7290239
    Xc                  -13.4577019
    Entropy             -0.0113289

    total               -107.668751427184

The fact that magnesium is a metal is confirmed by plotting the density of states around the Fermi level. To get better plots, we decrease the k spacing a bit for this step, i.e. we use a finer k-point mesh with more points.

bands = compute_bands(scfres, KgridSpacing(0.7 / u"Å"))
plot_dos(bands)