SCF Loop in Quantum Espresso

Summarized from Shobhana Narasimhan's presentation: PWscf: The SCF Loop and Some Relevant Input Parameters

Goal: Solve the Kohn-Sham equations:

$\left[\frac{1}{2}\nabla ^2 + V_{nuc}(r)+V_{H}[{n(r)}]+V_{xc}[n(r)] \right] \psi_i(r) = \varepsilon_i\psi_i(r) $

Since Hamiltonian, $H$, depends on solution $n(r)$, a self-consistent solution is required.

                                     Self-consistent Iterative Solution

=======================================================================
Q: Why using plane waves for periodic systems?
The most naive answer: as plane wave is periodic, it is nature to choose plane waves as the basis sets for periodic systems.

The wavefunction is a linear composition of plain waves,
$\psi_k(r)=\frac{1}{V}\sum_G c_{k,G}e^{i(k+G)r}$, where G is reciprocal lattice vector, k is the k points in BZ.
We can truncate the expansion at some value of |k+G|, as the contribution from higher Fourier components (large |k+G|) is usually small.
Related parameter: ecutwfc

$\frac{\hbar^2|k+G|^2}{2m}\le E_{cut}$

=======================================================================

0. Define the System

Related parameter: ibrav, celldm(i), nat, ntyp, 
Field: Atomic_Positions

1. Obtaining $V_{nuc}$

Given that (1) core wavefunctions sharply peaked near nucleus and (2) valence wavefunctions have lots of wiggles near nucleus, plane-wave basis sets with very high energy cutoff (|k+G|) are required to correctly describe the total wavefunctions. The problem is solved (computational expense is reduced) via using pseudopotentials. 
Related parameter: ecutwfc, ecutrho
Note:

"ecutwfc and ecutrho are the kinetic energy cutoff of the plane wave basis used to describe the wavefunctions and the charge density, respectively. The choice of the cutoffs depends only on the pseudopotentials (PPs) you are using and on which quantities you are interested in (for instance, magnetic properties are harder to converge wrt total energy). Once you've chosen your PPs you can test convergence of the energy wrt to cutoffs in the cheapest system (usually in the bulk). The threshold of convergence is somehow your choice: personally I don't go below 1mRy since this is the expected order of magnitude of the errors due to PPs."
For normconserving PPs you don't need to specify ecutrho, since the default (4*ecutwf) is good. For ultrasoft PPs you will need higher ecutrho due to the augmentation terms in the charge. "

2. Initial Guess for $n(r)$

Related parameter: startingwfc
a. superpositions of atomic densities, atomic;
b. converged n(r) from previous calculation, file;
c. random numbers, random.

3. $V_H$ and $V_{xc}$

Exchange-Correlation functions
LDA: pz
GGA: pw91, pbe

4. Potential & Hamiltonian

Kohn-Sham equations in plane wave basis

5. Diagonalization

Exact diagonaization can give $N_{PW}$ eigenvalues and eigenfunctions,  but the process is highly expensive as 
a. $N_{PW}$ is usually large
b. $N_{PW}^2$ elements to sore 
c. $N_{PW}^3$ operations are required to do the diagonalization
d. We usually only care limited number of bands 

Iterative diagonalization: recast diagonalization as a minimization process
Conjugate Gradient and Davidson
Related parameter: diagonalization, nbnd

6. New charge density: $n(r)=\sum_i|\psi_i(r)|^2$

Sum over all the k points in Brillouin Zone

k-point meshes:
1. Special points [Chadi & Cohen]
2. Monkhorst-Pack

Field: K_POINTS
Related parameter : nk1, nk2, nk2, k1, k2, k3

7. Test for Convergence 

a. compare new and old wavefunctions, charge densities, and total energies
b. compare with energy estimated with Harris-Foulkes functional
Related parameter: conv_thr (1.0e-8), electron_maxstep

8. Mixing to generate new $n(r)$

Take some combinations of input and output densities
 Related parameter: mixing_mode, mixing_beta

9. Output Quantities

$E_{tot}$: sum over eigenvalues, correct for double-counting of Hartree & XC terms, add ion-ion interactions
Force: Hellmann-Feynman theoreom