Matrix diagonalization¶
| Diagonalization method | eigmethod |
|---|---|
| LOBPCG | ‘lobpcg’ |
| Davidson | ‘davidson’ |
| Chebyshev filtering | ‘chebyfilt’ |
| Preconditioned Project Conjugate Gradient | ‘ppcg’ |
| Krylov Schur | ‘diagbyeigs’ |
To seek the fixed point of the Kohn-Sham map
we need to evaluate the matrix function \(f_\mu[\cdot]\), which is either a step function that jumps from 0 to 1 at \(\mu\), or a Fermi-Dirac distribution of the form
where :math`k_B` is the Boltzman constant and \(T\) is the temperature.
Such an evaluation can be performed by computing the eigenvalues and eigenvectors of the Kohn-Sham hamiltonian \(H\).
Because the evaluation of \(f_{\mu}\) does not need to be fully accurate when \(\rho\) or \(X\) is away from the ground state density or wavefunctions, and because a good initial guess from the previous SCF iteration is often available, an iterative diagonalization method is more appropriate.
KSSOLV2.0 provides several functions/methods for computing the lowest
\(k\) eigenpairs of a fixed Kohn-Sham Hamiltonian \(H\), where
\(k\) is specified by the user. It is typically set to either the
number of valence electron pairs or a number that is slightly larger
(when a finite temperature SCF calculation is performed.)
A user can choose a specific solver by setting the eigmethod
field of the options class to one of the choices listed in
the table above and passing options to the updateX() function.
Locally Optimal Block Conjugate Gradient lobpcg()¶
LOBPCG is the default eigensolver used in KSSOLV to compute the lowest \(k\) eigenpairs of a fixed Kohn-Sham Hamiltonian \(H\). The algorithm can be viewed as a way to solve the equivalent constrained minimization problem
Given a starting guess \(X^{(0)}\) to the minimizer, LOBPCG iterates using the following updating formula:
where \(C_1\), \(C_2\), and \(C_3\) are coefficient matrices chosen to minimize the trace of \(H\) within the space spanned by \(\{{X^{(i)}, W^{(i)}, P^{(i-1)}}\}\) while maintaining the orthonormality constraint \(\langle X^{(i+1)}, X^{(i+1)} \rangle = I\).
The matrix block \(W^{(i)}\) is the preconditioned gradient of the Lagrangian \(\frac{1}{2} X^* H X - (X^*X - I)\Theta\), where \(\Theta = X^* H X\), i.e.,
where \(K\) is an appropriately chosen preconditioner.
The matrix block \(P^{(i)}\) defines the \(i\) th search direction, which satisfies the recurrence formula
with \(P^{(1)}\) set to \(\emptyset\), and \(P^{(2)}\) set to \(W^{(1)}\).
The coefficient matrices \(C_1, C_2, C_3\) are obtained by solving a projected eigenvalue problem
where \(Y = \left(X^{(i)}, W^{(i)}, P^{(i-1)}\right)\), \(C = (C_1^*, C_2^*, C_3^*)^*\), and \(\Theta\) is a diagonal matrix containing the \(k\) smallest eigenvalues of the project Hamiltonian. The solution of the projected eigenvalue problem (3) and the udpate of the approximate eigenvector in eq:update are often collectively called a Rayleigh Ritz procedure.
The LOBPCG solver can be called as follows
>> [X,lambda] = lobpcg(H, X0, prec, tol, maxit, verbose);
where H is a Hamiltonian object, X0 is a wavefunction object that stores the initial guess to the desired eigenvectors, prec is a wavefunction object that stores a diagonal preconditioner (in reciprocal space), tol is a user specified convergence tolerance, maxit is the maximum number of iterations allowed, and verbose controls the amount of diagnostic output during the LOBPCG iteration.
Krylov Schur diagbyeigs()¶
KSSOLV2.0 can use MATLAB’s sparse eigensolver eigs() which implements
the implicitly restarted Lanczos (or Krylov Schur) method. This is not
the fastest solver because it does not take full advantage of approximations
to several eigenvectors produced from a previous SCF iteration. Also, because
the Hamiltonian is multiplied with one wavefunction (coefficients) at a time,
there is limited opportunity for parallelization or data movement optimization.
The Kyrlov Schur solver can be called as follows
>> [X,lambda] = diagbyeigs(H, neig, tol, maxit);
where neig is the number of eigenpairs to be computed, and maxit is the maximum number of implicit restart used in the Krylov Schur algorithm.
Chebyshev Filtering chebyfilt()¶
In Chebyshev filtering, we use a subspace iteration to compute dominant eigenvalues of \(T_d(H;\alpha,\lambda_{\mathrm{ub}})\), where \(T_d(t;\alpha,\lambda_{\mathrm{ub}})\) is a \(d\) th degree shifted and scaled Chebyshev polynomial that is bounded by -1 and 1 on the interval \([\alpha,\lambda_{\mathrm{ub}}]\), and \(\alpha\) is a carefully chosen parameter that is slightly above the largest eigenvalue of \(H\) that we need to compute (e.g., the largest occupied state), and \(\lambda_{\mathrm{ub}}\) is an upper bound of the spectrum of \(H\). The polynomial is chosen to amplify eigenvalues to the left of \(\alpha\).
Chebyshev filtering can be called as follows
>> [X,lambda] = chebyfilt(H, X0, degree, lb, ub);
where degree is the degree of the polynomial, lb is the parameter \(\alpha\) described above and ub is the upper bound of the spectrum \(\lambda_{\mathrm{ub}}\).
Davidson davidson()¶
To some extent, the Davidson’s algorithm is a generalization of the LOBPCG algorithm. In the i th iteration, the LOBPCG algorithm extracts approximations to the desired eigenpairs from a subspace spanned by \(\{ W^{(i)}, X^{(i)}, X^{(i-1)} \}\). In a Davidson’s algorithm, we can extract approximation from a larger subspace
However, as the dimension of the subspace becomes larger, the cost of solving the projected eigenvalue problem as well as maintaining an orthonormal basis of such as subspace can become much higher, and the algorithm can be come unstable due to the difficult in finding an orthonormal basis of the subspace.
PPCG ppcg()¶
The preconditioned projected conjugate gradient (PPCG) algorithm is a variant of the LOBPCG designed to reduce the cost of solving the projected problem (3) when the number of eigenpairs to be computed is relatively large. If \(X^{(i)}\) is sufficiently close to the solution of (1), we can partition \(X^{(i)}\) into several smaller blocks of vectors and apply LOBPCG to each block separately. The orthonormalization of all vectors in \(X^{(i)}\) and the Rayleigh Ritz procedure is applied periodically to prevent vectors in different blocks from converging to the same set of eigenvectors.