Energy cutoff smearing

A technique that has been employed in the literature to ensure smooth energy bands for finite Ecut values is energy cutoff smearing.

As recalled in the Problems and plane-wave discretization section, the energy of periodic systems is computed by solving eigenvalue problems of the form

In [1]:
#$$
\begin{equation} H_ku_k=\varepsilon_k u_k \end{equation}

$$

for each $k$-point in the first Brillouin zone of the system. Each of these eigenvalue problem is discretized with a plane-wave basis $\mathcal{B}_k^{E_c}=\{x\mapsto e^{iG\cdot x} \;\;|\;G\in\mathcal{R}^*,\;\; |k+G|^2\leq 2E_c\}$ whose size highly depends on the choice of $k$-point, cutoff energy $\rm E_c$ or cell size. As a result, energy bands computed along a $k$-path in the Brillouin zone or with respect to the system's unit cell volume - in the case of geometry optimization for example - display big irregularities when $E_c$ is taken too small.

Here is for example the variation of the ground state energy of face cubic centered (FCC) silicon with respect to its lattice parameter, around the experimental lattice constant.

In [2]:
using DFTK
using Statistics

a0 = 10.26 # Experimental lattice constant of silicon in bohr
a_list = range(a0 - 1/2, a0 + 1/2; length=20)

Ecut    = 5         # very low Ecut to display big irregularities
kgrid   = (2, 2, 2) # very sparse k-grid to fasten convergence
n_bands = 8         # Standard number of bands for silicon

function compute_ground_state_energy(a; Ecut, kgrid, kinetic_blowup, kwargs...)
    lattice = a / 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]
    model = model_PBE(lattice, atoms, positions; kinetic_blowup)
    basis = PlaneWaveBasis(model; Ecut, kgrid)
    self_consistent_field(basis; kwargs...).energies.total
end

callback = info->nothing # set SCF to non verbose
E0_naive = compute_ground_state_energy.(a_list; kinetic_blowup=BlowupIdentity(),
                                        Ecut, kgrid, n_bands, callback);

To be compared with the same computation for a high Ecut=100. The naive approximation of the energy is shifted for the legibility of the plot.

In [3]:
E0_ref = [-7.839775223322127, -7.843031658146996, -7.845961005280923,
          -7.848576991754026, -7.850892888614151, -7.852921532056932,
          -7.854675317792186, -7.85616622262217,  -7.85740584131599,
          -7.858405359984107, -7.859175611288143, -7.859727053496513,
          -7.860069804791132, -7.860213631865354, -7.8601679947736915,
          -7.859942011410533, -7.859544518721661, -7.858984032385052,
          -7.858268793303855, -7.857406769423708]

using Plots
shift = mean(abs.(E0_naive .- E0_ref))
p = plot(a_list, E0_naive .- shift, label="Ecut=5", xlabel="lattice parameter a (bohr)",
         ylabel="Ground state energy (Ha)", color=1)
plot!(p, a_list, E0_ref, label="Ecut=100", color=2)
Out[3]:

The problem of non-smoothness of the approximated energy is typically avoided by taking a large enough Ecut, at the cost of a high computation time. Another method consist in introducing a modified kinetic term defined through the data of a blow-up function, a method which is also referred to as "energy cutoff smearing". DFTK features energy cutoff smearing using the CHV blow-up function introduced in [REF of the paper to be submitted], that is mathematically ensured to provide $C^2$ regularity of the energy bands.

Let us lauch the computation again with the modified kinetic term.

In [4]:
E0_modified = compute_ground_state_energy.(a_list; kinetic_blowup=BlowupCHV(),
                                           Ecut, kgrid, n_bands, callback, );

!!! note "Abinit energy cutoff smearing option" For the sake of completeness, DFTK also provides the blow-up function proposed in the Abinit 1 quantum chemistry code. This function depends on a parameter Ecutsm fixed by the user. Note that for the right choice of Ecutsm, the Abinit blow-up function corresponds to the CHV one with a choice of coefficients ensuring $C^1$ regularity. The Abinit options is chosen with kinetic_blowup=BlowupAbinit(Ecutsm)


  1. Abinit software suite user guide

We can know compare the approximation of the energy as well as the estimated lattice constant for each strategy.

In [5]:
estimate_a0(E0_values) = a_list[findmin(E0_values)[2]]
a0_naive, a0_ref, a0_modified = estimate_a0.([E0_naive, E0_ref, E0_modified])

shift = mean(abs.(E0_modified .- E0_ref))  # again, shift for legibility of the plot

plot!(p, a_list, E0_modified .- shift, label="Ecut=5 + BlowupCHV", color=3)
vline!(p, [a0], label="experimental a0", linestyle=:dash, linecolor=:black)
vline!(p, [a0_naive], label="a0 Ecut=5", linestyle=:dash, color=1)
vline!(p, [a0_ref], label="a0 Ecut=100", linestyle=:dash, color=2)
vline!(p, [a0_modified], label="a0 Ecut=5 + BlowupCHV", linestyle=:dash, color=3)
Out[5]:

The smoothed curve obtained with the modified kinetic term allow to clearly designate a minimal value of the energy with respect to the lattice parameter $a$, even with the low Ecut=5 Ha.

In [6]:
println("Error of approximation of the reference a0 with modified kinetic term:"*
        " $(round((a0_modified - a0_ref)*100/a0_ref, digits=5))%")
Error of approximation of the reference a0 with modified kinetic term: 0.50393%