Self Consistent Field Iteration

After the Kohn-Sham equations are discretized by planewaves, they become a set of nonlinear algebraic eigenvalue equations

\[H(\rho(X)) X = X\Lambda\]

in which the matrix to be diagonalized is a function of the eigenvectors \(X\) to be computed. The desire eigenvalues to be computed (on the diagonal of \(\Lambda\)) are the \(n_e\) smallest eigenvalues of \(H\), where \(n_e\) is proportional to the number of electrons. Due to this nonlinear coupling between the matrix operator and its eigenvectors, the Kohn-Sham equations are much more difficult to solve than a standard linear eigenvalue problem.

Currently, the most widely used numerical method for solving the this type of problem is the Self Consistent Field (SCF) iteration. An alternative approach which is designed to minimize the Kohn-Sham total energy directly will be reviewed later. Both methods have been implemented in KSSOLV to allow users to compute the solution to the Kohn-Sham equations associated with various molecules and solids.

To use the SCF algorithm in KSSOLV, we can simply type:

>> [mol, info] = scf(mol);

where mol is a previously constructed Molecule object. The function call returns a Molecule object which may contain updated atomic forces. The total energy of the molecule obtained at the end of the SCF procedure, as well as the eigenvalues and eigenvectors of the Kohn-Sham Hamiltonian are stored in the info structure.

Given an initial guess of the eigenvector matrix \(X\) and electron density \(\rho\)

[mol, H, X, Hprec, nocc] = scfinit(mol);

Construct initial Hamiltonian in scfinit function associated with mol .:

H = Ham(mol);

H is a Ham object that contains information about kinetic energy and total potential(excluding the non-local ionic pseudopotential). Initial guess of eigenvector matrix \(X\)

X = genX0(mol);

X is a Wavefun object that contains information about grids and factors of wave functions in Fourier space. Generates the preconditioner for the Hamiltonian, which is in the format of Wavefun.:

Hprec = genprec(H);

And nocc is the number of occupy state. Calculate Ewald and Ealphat and then begin to SCF iteration.:

Eewald  = getEewald(mol);
Ealphat = getEalphat(mol);

the SCF procedure can be described by the following loop

iter = 1
while (not converged and iter < maxiter)
   Compute the ne smallest eigenvalues and corresponding eigenvector of H
   Update the eletron density, potential and H
   iter = iter + 1
end while

In SCF interation, KSSOLV compute \(X^{i+1}\) given by \(H^{i}X^{i+1}=X^{i+1}\Lambda^{i+1}\):

[X,ev] = KSeigs(mol, H, X, Hprec, options);

Update density function \(\rho\) by getting square modulus of the wave functions in real space \(\rho^{i+1}=diag(X^{i+1}*X^{i+1})\) and update H, using the \(\rho\) coming from the computation of the new eigen functions.:

rho = getcharge(mol, X, nocc);

Mix density or mix potential later:

[rho, dfmat, dvmat, cdfmat] = potmixing(mol, rhoin, rho, iterscf, mixtype, betamix, dfmat, dvmat, cdfmat, mixdim, brank);

Kinetic energy and some additional energy terms are calculated by:

Ekin = (2/nspin)*sum(ev(1:nocc));

Ionic and external potential energy was included in Ekin along with incorrect Ecoul and Exc. Need to correct them later. Compute an energy correction term:

Ecor = getEcor(mol, rho, vtot, vion, vext);

Compute Hartree and exchange correlation energy and potential using the new charge density and update the total potential:

[vhart, vxc, uxc2, rho, uxcsr] = getVxc(mol,rho);

Update total potential:

vtot = getVtot(vion,vext,vhart,vxc);

Mix potential if user didn’t choose density mixing:

[vtot, dfmat, dvmat, cdfmat] = potmixing(mol, vtotin, vtot,iterscf, mixtype, betamix, dfmat, dvmat, cdfmat, mixdim, brank);

Calculate the coulomb and exchange correlation potential energy based on the new potential:

Ecoul = getEcoul(mol, abs(rho), vhart);
Exc   = getExc(mol, abs(rho), uxc2);

Calculate total energy in this SCF iteration:

Etot = Eewald + Ealphat + Ekin + Ecor + Ecoul + Exc;

Check convergence and break loop either SCF converge or reach the max number of iteration steps. After SCF iteration, calculate the atom force:

mol.xyzforce = getFtot(mol, H, X, rho);

However, this simple procedure rarely converges because it corresponds to a simple fixed point iteration for a set of nonlinear equations satisfied by the electron density

\[\rho = f_\mu (\rho), \sum_i \rho(i) = n_e,\]

where the function \(f_\mu(\rho)\) can be written as the diagonal part of a step function applied to \(H(\rho)\). The step function assume the value of 1 on the interval \((-\infty,\mu)\) and 0 on the interval \([\mu,\infty]\). The value of \(\mu\) is called the chemical potential. It is not uniquely defined by the above equation if there is a gap between the \(n_e\)-th and the \(n_e+1\)-st eigenvalues of the converged \(H(\rho)\).

Because \(f_\mu(\rho)\) is generally not a global contraction, there is no guarantee that a simple fixed point iteration that starts from any initial guess of \(\rho\) can converge. A more sophisticated iterative algorithm is required to solve the nonlinear system successfully. We will discuss this alogrithm in the next section.