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)`

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))%")
```