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.
Both the Ecut and the kgrid keyword argument in PlaneWaveBasis are optional. If the user does not specify these values, DFTK will try to determine reasonable default values:
- Kgrid default:
kgrid=KgridSpacing(2π * 0.022 / u"bohr"), which usually gives reasonable results for a first calculation. - Ecut default: DFTK will consult the PseudoPotentialData.jl library for a recommended kinetic energy cutoff and use the maximal value over all atoms of the calculation. See the Pseudopotentials chapter for more details on using pseudopotentials with DFTK. For cases where no recommended values can be determined, DFTK will throw an error and expects the user to manually provide a value for
Ecut.
Therefore we could also construct a more reasonable basis as follows:
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
Estimated memory usage (per MPI process):
nonlocal projectors : 0.98MiB
single ψ : 874KiB
single ρ : 750KiB
peak memory (SCF) : 16.6MiB
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]. In production calculations it is often advisable to refine these values by Performing a convergence study.
The above detailed basis printing with parallelisation information, discretisation parameters, structure and model is the default output generated when returning a basis in an interactive session or a Pluto / Jupyter notebook. In order to obtain the same information in a script (e.g. to stdout) execute
show(stdout, "text/plain", basis)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.0209913727 0.39 5.5 96.0ms
2 -107.6655344435 -0.19 -0.47 4.2 89.8ms
3 -107.6687199919 -2.50 -1.59 2.7 76.0ms
4 -107.6687498424 -4.53 -2.76 1.0 52.9ms
5 -107.6687513383 -5.83 -3.38 4.0 77.4ms
6 -107.6687514218 -7.08 -3.87 2.5 64.7ms
7 -107.6687514266 -8.32 -4.57 4.3 85.4ms
8 -107.6687514272 -9.28 -5.00 2.7 68.4ms
9 -107.6687514272 -10.64 -6.11 2.0 61.2msscfres.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.9999272441691354
1.9121679887267309
6.511418155317659e-17
6.510816118632296e-17
7.050056190740388e-18
7.37737728241726e-19
7.376802725128064e-19
1.7854031377715975e-19scfres.energiesEnergy breakdown (in Ha):
Kinetic 48.3202528
AtomicLocal -101.1884735
AtomicNonlocal -19.3973028
Ewald -53.8605542
PspCorrection 1.1973309
Hartree 30.7290267
Xc -13.4577025
Entropy -0.0113289
total -107.668751427183The 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)