6 Self Consistency


6.1 What are the units for quantity XYZ?

Unless otherwise specified, all PWscf input and output quantities are in atomic “Rydberg” units, i.e. energies in Ry, lengths in Bohr radii ($a_0$), etc.. Note that CP uses instead atomic “Hartree” units: energies are in Ha, not in Ry. The charge density is in $a_0^3$ units, i.e. if integrated on the cell, it yields the total number of electrons.

6.2 Self-consistency is slow or does not converge at all

In most cases: your input data is bad, or else your system is metallic and you are treating it as an insulator. If this is not the case: reduce mixing_beta to 0.3 ÷ 0.1 or smaller, try the mixing_mode value that is more appropriate for your problem.

6.3 What is the difference between total and absolute magnetization?

The total magnetization is the integral of the magnetization in the cell:

$$ M_{Tot}= \int_{cell}{ \left( n_{up}-n_{dw} \right) d^3r} $$

The absolute magnetization is the integral of the absolute value of the magnetization in the cell:

$$ M_{Abs} = \int_{cell}{ \lvert n_{up}-n_{dw} \rvert d^3 r} $$

In a simple ferromagnetic material they should be equal (except possibly for an overall sign). In simple antiferromagnets (like FeO, NiO) $M_{Tot}$ is zero and $M_{Abs}$ is twice the magnetization of each of the two atoms.

6.4 How can I calculate magnetic moments for each atom?

There is no ‘right’ way of defining the local magnetic moment around an atom in a multi-atom system. However an approximate way to define it is via the projected density of states on the atomic orbitals (code projwfc.x). This code generates many files with the density of states projected on each atomic wavefunction of each atom and a BIG amount of data on the standard output, the last few lines of which contain the decomposition of Löwdin charges on angular momentum and spin component of each atom.

6.5 What is the order of Ylm components in projected DOS / projection of atomic wavefunctions?

The order of m-components for each l in the output is:
$ 1, cos(\phi), sin(\phi), cos(2\phi), sin(2\phi), \dotsc ,cos(l\phi), sin(l\phi) $ See input data documentation for projwfc.x.

6.6 Why is the sum of partial Löwdin charges not equal to the total charge?

“Löwdin charges (as well as other conventional atomic charges) do not satisfy any sum rule. You can easily convince yourself that this is the case because the atomic orbitals that are used to calculate them are arbitrary to some extent. If you like, you can think that the missing charge is “delocalized” or “bonding” charge, but this would be another way of naming the conventional (to some extent) character of Löwdin charge.” (Stefano Baroni, Sept. 2008).

See also the definition of “spilling parameter”: Sanchez-Portal et al., Sol. State Commun. 95, 685 (1995). The spilling parameter measures the ability of the basis provided by the pseudo-atomic wavefunctions to represent the eigenstates of the crystal, by measuring how much of the subspace of the Hamiltonian eigenstates falls outside the subspace spanned by the atomic basis.

6.7 I cannot find the Fermi energy, where is it?

It is printed in the output. If not, the information on Gaussian smearing, needed to calculate a sensible Fermi energy, was not provided in input. In this case, pw.x prints instead the highest occupied and lowest unoccupied levels. If not, the number of bands to be calculated was not provided in input and pw.x calculates occupied bands only.

6.8 What is the reference level for Kohn-Sham energies? Why do I get positive values for Kohn-Sham levels?

The reference level is an ill-defined quantity in calculations in solids with periodic boundary conditions. Absolute values of Kohn-Sham eigenvalues are meaningless.

6.9 Why do I get a strange value of the Fermi energy?

“The value of the Fermi energy (as well as of any energy, for that matter) depends upon the reference level. What you are referring to is probably the “Fermi energy referred to the vacuum level” (i.e. the work function). In order to obtain that, you need to know what the vacuum level is, which cannot be said from a bulk calculation only”> (Stefano Baroni, Sept. 2008).

6.10 Why I don’t get zero pressure/stress at equilibrium?

If you make a calculation with fixed cell parameters, you will never get exactly zero pressure/stress, unless you use the cell that yields perfect equilibrium for your pseudopotentials, cutoffs, k-points, etc.. Such cell will anyway be slightly different from the experimental one. Note however that pressures/stresses in the order of a few KBar correspond to very small differences in terms of lattice parameters.

6.11 Why do I get different results from vc-relax and from scf on the same structure?

First of all, you should verify that the structure is really the same (hint: compare Ewald energies). The PW basis set used in a variable-cell calculations is determined by the cutoff and by the initial cell geometry. If you make a calculation with the final geometry at the same cutoff, you get slightly different results because the PW basis set is different. The difference should be small, though, unless you are using a too low cutoff for your system. Since v.4.3.1, a final scf step is performed at the end of the vc-relax run, with the PW basis set computed for the final geometry, to check for this.
Please also note that if you use a modified kinetic energy functional for variable-cell calculations, this affects the calculated pressure/stress.

6.12 Why do I get negative starting charge?

Self-consistency requires an initial guess for the charge density in order to bootstrap the iterative algorithm. This first guess is usually built from a superposition of atomic charges, constructed from pseudopotential data.

More often than not, this charges are a slightly too hard to be expanded very accurately in PWs, hence some aliasing error will be introduced. Especially if the unit cell is big and mostly empty, some local low negative charge density will be produced.

”This is NOT harmful at all, the negative charge density is handled properly by the code and will disappear during the self-consistent cycles”, but if it is very high (let’s say more than 0.001*number of electrons) it may be a symptom that your charge density cutoff is too low.” (L. Paulatto – November 2008)

6.13 How do I calculate the work function?

Work function = (average potential in the vacuum) – (Fermi Energy). The former is estimated in a supercell with the slab geometry, by looking at the average of the electrostatic potential (typically without the XC part). There is an example showing how to calculated the work function. 

6.14 Why are my forces zero / why BFGS does not relax my crystal?

In infinite systems there is a fundamental difference between displacements that do not change the unit cell (phonon modes at q=0) and those that do (elastic modes). The former are driven by forces, the latter by stresses. Equilibrium is reached when both forces and stresses are zero. In symmetric crystals, forces are often zero by symmetry. In that case, the BFGS algorithm, or any other algorithm that ignores stresses and keeps the unit cell constant, will do nothing.

6.15 My graphene sheet has a gap!

Congratulations! Get a suit for the Nobel price ceremony. No, wait: first check if your input is correct, FAQ 3.8.