Are CHARGE=2 and INIWAV=1 contradictory in the VASP setting?

Question

When initializing the charge density in VASP, I set the parameters as follows:

Startparameter for this run:
ISTART =0         #job   :0-new  1-continue 2-samecut   charge:0-initial
ICHARG =2         #charge:0-initial orbitals 1-from CHGCAR 2-ISTART=0
INIWAV =1         #0-lowe 1-random: only ISTART=0 use! default-1

As shown on the wiki page of ICHARG and INIWAV,

ICHARG=2    #Take superposition of atomic charge densities
INIWAV=1    #Fill wavefunction arrays with random numbers.

ICHARG=2 looks like it tells VASP to store the atomic charge density from POTCAR files, and add them together to form the initial charge density.

INIWAV=1 seems to assume the wavefunctions in the unit cell are plane waves, and sets random occupation number for each wavefunction.

Aren't these two options contradictory?

Answer 1

If you solve the Kohn-Sham equations using density mixing (or any non-self-consistent method) then the density and wavefunctions are treated independently. This means that you are free to initialise the density and wavefunctions independently, and in this case it appears that the density is initialised from a sum of (pseudo-)atomic charge densities, and the wavefunctions are initialised randomly. This is quite a common initialisation.

You might wonder how the wavefunctions and densities can be treated independently and still optimised to give the ground state. The usual approach is the "self-consistent field" method (SCF) where the wavefunctions are optimised by solving the Kohn-Sham equations (iteratively, in the case of a plane-wave-based code like VASP) for a fixed input density,

$\mathrm{H}_{k}[\rho_{in}]\psi_{bk}(r)=\epsilon_{bk}\psi_{bk}(r).$

Once the optimised wavefunctions have been calculated, an "output density" is computed from them,

$\rho_{out}(r) = \sum_{bk} f_{bk}|\psi_{bk}(r)|^2$,

where $b$ and $k$ are band and k-point indices, respectively, and $f_{bk}$ are the band occupancies.

If $\rho_{out}=\rho_{in}$, then $\{\psi_{bk}\}$ are the correct ground state Kohn-Sham states; however, this is not the case in general, and a new trial input density is contructed from $\rho_{out}$ and $\rho_{in}$ according to a different optimisation method (for example, a Kerker-Pulay mixing method). This constitutes one "SCF cycle".

The updated density is then used as the input density for the next SCF cycle, and the wavefunctions optimised for this new density. As the SCF cycles progress, the goal is for the input and output densities to converge towards each other to within some specified precision, at which point it is assumed that $\{\psi_{bk}\}$ are sufficiently accurate to be treated as the ground state Kohn-Sham states. At this point $\rho_{out}\approx \rho_{in}$ and the solution is said to be (approximately) "self-consistent".