Your own crystal

In this page we compute a crystal that is not in the repository, and for concreteness we take diamond. We have to choose four things before a number from a plane-wave code has a meaning:

1. The structure, which is the Bravais lattice, the lattice constant and the atomic positions.
2. The pseudopotential of each element, which also fixes the functional.
3. The plane-wave cutoff.
4. The k-mesh, which says how finely the Brillouin zone is sampled.

The first two are physics. The last two are numerical parameters, so we have to converge them for the quantity we want.

The input file

We start with the input file of diamond:

docs/figures/own-crystal_convergence/diamond.in

! Diamond: two carbon atoms on the fcc lattice, a = 3.567 A = 6.74 bohr.
! LDA, ultrasoft. ecutrho is eight times ecutwfc, as an ultrasoft dataset needs.
! nbnd = 8 gives four occupied and four empty bands.
! The input of the figure "own-crystal_convergence".
 &control
    calculation = 'scf'
 /
 &system
    ibrav = 2, celldm(1) = 6.74, nat = 2, ntyp = 1,
    ecutwfc = 35.0, ecutrho = 280.0, nbnd = 8
 /
 &electrons
    conv_thr = 1.0d-8
 /
ATOMIC_SPECIES
 C  12.011  C.pz-rrkjus.UPF
ATOMIC_POSITIONS crystal
 C 0.00 0.00 0.00
 C 0.25 0.25 0.25
K_POINTS (automatic)
 6 6 6 0 0 0

It is a normal pw.x input file, and it means the same thing as in Quantum ESPRESSO. ibrav = 2 is the fcc lattice, celldm(1) is the cubic lattice constant in bohr, and the two atoms are given in crystal coordinates. An input file written for pw.x runs here without changes, and an input file written here runs in pw.x. We compute its total energy with diamond.get_scf():

from defuqumat import Calculator

diamond = Calculator.from_file("docs/figures/own-crystal_convergence/diamond.in",
                               pseudo_dir="tests/data/pseudo")
print(diamond.get_scf().total_energy)   # in Ry

The pseudopotentials are UPF files, and the code reads their names from the ATOMIC_SPECIES card and looks for them in pseudo_dir. If you do not give pseudo_dir, the code uses the pseudo_dir of the &control namelist, and then the folder of the input file.

Where to get a pseudopotential

The repository has only the pseudopotentials that its tests need, in tests/data/pseudo/. For other elements you can take them from three common sources:

• The table of Quantum ESPRESSO, at pseudopotentials.quantum-espresso.org.
• The PseudoDojo set.
• The SSSP library, which also gives a recommended cutoff for each element.

The code reads norm-conserving, ultrasoft and PAW files. The header of the file says which functional was used to generate it, and we use this functional unless you set input_dft. Notebook 05 shows what changes between LDA and PBE.

The functionals we have are LDA (PZ), PBE, revPBE and PBEsol. We write down only the energy of the functional, and the code takes the potential as its gradient. Two more functionals work the other way around. input_dft = 'tb09' (Tran and Blaha) and input_dft = 'bj06' (Becke and Johnson) are meta-GGA potentials with no energy behind them. They are used for band gaps, and for silicon tb09 gives a gap of 1.13 eV, where LDA gives 0.49 eV.

The two convergence tests

Let us now converge the two numerical parameters. A plane-wave basis has one parameter, the cutoff ecutwfc, and the integral over the Brillouin zone has one parameter, the k-mesh. Neither is a property of the crystal, so we cannot look them up, and we measure them instead.

The input is text, so a convergence test is just a loop over a string replacement:

text = open("docs/figures/own-crystal_convergence/diamond.in").read()

for ecut in (20, 30, 40, 50, 60):
    new = text.replace("ecutwfc = 35.0, ecutrho = 280.0",
                       f"ecutwfc = {ecut}, ecutrho = {8 * ecut}")
    calc = Calculator.from_text(new, pseudo_dir="tests/data/pseudo")
    print(ecut, calc.get_scf().total_energy)

In the left panel we show the total energy against ecutwfc, measured from the value at 60 Ry, on the \(6\times6\times6\) mesh of the file. In the right panel we show the total energy against the mesh \(n\times n\times n\), measured from the \(10\times10\times10\) mesh, at the 35 Ry of the file. The dashed line is \(10^{-4}\) Ry, which is roughly the accuracy that an energy difference between two structures needs.

We see that the cutoff converges from above, because a larger basis can only lower a variational energy. For this carbon dataset it converges slowly, and at 35 Ry the total energy is still \(10^{-3}\) Ry above the value at 60 Ry. The k-mesh converges fast, because diamond is an insulator. A metal has a Fermi surface, so it needs many more k-points for the same accuracy.

What you have to converge is the quantity that you will report, not only the total energy. A band gap or a lattice constant is often converged at a lower cutoff than the total energy, because the errors cancel in a difference.

A band gap from the SCF mesh

The smallest gap among the k-points of the SCF mesh is not the band gap. The conduction band minimum of diamond is at about 85% of the way from \(\Gamma\) to X, and this point is on none of the meshes. So we compute the gap from a band structure along a k-path, as in Band structure.

Code of this figure (docs/figures/own-crystal_convergence/)

compute.py

"""The total energy of diamond against the cutoff and against the k-mesh."""
import numpy as np
from defuqumat import Calculator
from common import PSEUDO, save

text = open("diamond.in").read()       # 35 Ry and a 6x6x6 mesh


def energy(text):
    return Calculator.from_text(text, pseudo_dir=PSEUDO).get_scf().total_energy


cutoffs = [20, 25, 30, 35, 40, 50, 60]     # ecutwfc in Ry, with ecutrho = 8 ecutwfc
e_cutoff = [energy(text.replace("ecutwfc = 35.0, ecutrho = 280.0",
                                f"ecutwfc = {c:.1f}, ecutrho = {8 * c:.1f}"))
            for c in cutoffs]

meshes = [2, 4, 6, 8, 10]                  # n x n x n at 35 Ry
e_mesh = [energy(text.replace(" 6 6 6 0 0 0", f" {n} {n} {n} 0 0 0"))
          for n in meshes]

save(cutoffs=np.array(cutoffs, float), e_cutoff=np.array(e_cutoff),
     meshes=np.array(meshes), e_mesh=np.array(e_mesh))

plot.py

"""How far the energy is from its last value, against the cutoff and the k-mesh."""
import numpy as np
import style

d = np.load("data.npz")


def draw():
    fig, (a, b) = style.figure(2)
    for ax, x, e, label in ((a, d["cutoffs"], d["e_cutoff"], "ecutwfc (Ry)"),
                            (b, d["meshes"], d["e_mesh"], r"k-mesh $n\times n\times n$")):
        ax.semilogy(x[:-1], np.abs(e[:-1] - e[-1]), "o-")
        ax.axhline(1e-4, **style.reference(), label=r"$10^{-4}$ Ry")
        ax.set_xlabel(label)
        ax.legend(loc="upper right")
    a.set_ylabel(r"$|E - E_{\rm last}|$ (Ry)")
    b.set_xticks(d["meshes"][:-1])
    return fig


style.main(draw)

Changing the crystal from Python

We can also change the crystal without writing a new input file. Three methods give a new calculator with a changed crystal, and the old calculator does not change. When this makes sense, each method passes the converged density to the new calculator as a starting guess.

• calc.with_positions(positions): the same cell with the atoms moved. The positions are cartesian, in bohr, with shape (atoms, 3). The density of the old geometry is the starting guess, so the new SCF needs fewer iterations.
• calc.with_cell(at, positions=None): a deformed cell. at holds the three lattice vectors in bohr, as rows. The density is not used as a starting guess, because a new cell has a new grid.
• calc.with_kpoints(kpoints): the same crystal on another k-set, for example the full mesh next to its irreducible wedge. The two must give the same energy, so this is a check of the symmetry reduction.

Before you trust a new crystal

Before we trust the numbers of a new crystal, we check four things.

• The functional that the code reports. It comes from the header of the pseudopotential file, unless input_dft sets it.
• The symmetry that the code found. An atom typed one thousandth of a bohr off its site removes symmetry operations, and then the run uses more k-points and takes longer, but the result is the same.
• The vacuum, for a slab or a molecule. It must be thick enough that the periodic images do not see each other, and this is a convergence test of its own.
• The number of bands. An insulator needs the occupied bands and the ones you want to look at, and a metal needs bands well above the Fermi level.

Limits

An input variable that changes the physics and is not implemented gives an error with its name, so it is never read and ignored. These are a charged cell (tot_charge), the corrections for an isolated system (assume_isolated), two chemical potentials (twochem, in &system or in &control), the occupations of an isolated atom (one_atom_occupations), the constant-cutoff kinetic functional (qcutz) and a fixed FFT grid (nr1 to nr3s). The four symmetry switches no_t_rev, force_symmorphic, use_all_frac and nosym_evc also give an error, because each one removes symmetry operations that the code would still use. A value equal to the default of pw.x gives no error, so an ordinary input never meets these. calculation = 'md' and calculation = 'vc-md' give an error, because we have no molecular dynamics.

The old DFT+U variables of &system give an error, as in pw.x 7.5. These are lda_plus_u, lda_plus_u_kind, Hubbard_U, Hubbard_U_back, Hubbard_J0, Hubbard_J, Hubbard_V, U_projection_type, Hubbard_parameters and backall, and you write the HUBBARD card instead. The test is on the value, so lda_plus_u = .false. and Hubbard_U = 0.0 give no error.

A lattice given twice gives an error, as in pw.x. ibrav different from 0 together with a CELL_PARAMETERS card gives an error, and to use the card you set ibrav = 0. celldm together with A gives an error. A CELL_PARAMETERS bohr or angstrom card together with celldm(1) or A gives an error, and you either write the card in alat units or remove celldm(1) or A.

Only UPF version 2 files are read. A UPF version 1 file gives an error, and you can convert it with upfconv.x of Quantum ESPRESSO. Pseudopotentials generated with different functionals give an error, unless you set input_dft, and then the code warns that the functional differs from the files.

A functional that is not implemented gives an error, and the code never replaces it with LDA. This covers the hybrid functionals (EXX) and the meta-GGAs with an energy (TPSS, SCAN, M06L). tb09 and bj06 give an error with an ultrasoft dataset, because it has no partial waves to build the kinetic energy density inside the spheres, and you use a PAW or norm-conserving one instead. They also give an error with a spin spiral and with DFT+U (a HUBBARD card). With tb09 or bj06 the code warns that the total energy is not the energy of any functional, and forces, stress, phonons and the responses give an error.

The notebook for this page is 25, your own crystal.