Regularization through partial occupation¶
For systems with \(n_e\) valence electron pairs, the Kohn-Sham problem is ill-posed if the multiplicity of the \(n_e\) th eigenvalue of of the Kohn-Sham Hamiltonian \(H(\rho)\) is larger than 1. The ill-posedness of the problem can be seen from the ambiguity of definition of electron density
where \(\psi_i\) is the \(i\) th eigenfunction associated with the \(i\) th eigenvalue \(\varepsilon_i\), ordered as \(\varepsilon_1 \leq \varepsilon_2 \leq \cdots \leq \varepsilon_n\).
The degeneracy or near degeneracy of eigenvalues at the Fermi level \(\mu\) is a feature of metallic systems that exhibit no gap or a tiny gap between the valence (occupied) and conducting (unoccupied) states.
Computationally, the ill-poseness of such problems makes it difficult for the standard SCF iteration to converge. One often experience the so-called “charge sloshing” problem where the largest magnitude of \(\rho\) oscillates among a few spatial locations.
To overcome this problem, a regularization technique is used to modify the fixed point equation to be solved as
where \(f_\beta\) is the Fermi-Dirac distribution function defined as
The parameter \(\beta\) is a regularization parameter that is related to the temperature \(T\) in the original definition of the Fermi-Dirac distribution via
where \(k_B\) is the Boltzman constant.
At zero temperature, \(\beta = \infty\), and \(f_\beta\) becomes a step function that is 1 for \(t \leq \mu\) and 0 for \(t > \mu\). This corresponds to the standard Kohn-Sham model in which all eigenvalues below \(\mu\) are occupied. If there is a gap between occupied and unoccupied states, \(f_\beta(\varepsilon_i) = 1\) for \(i = 1,2,...,n_e\), and \(f_\beta(\varepsilon_i) = 0\) for \(i > n_e\).
When \(T > 0\) (often known as finite temperature in physics literature), \(0 \leq f_{\beta}(\varepsilon_i) \leq 1\), i.e., electronic states can be partially occupied. In particular, it is possible that \(f_{\beta}(\varepsilon_i) > 0\) for \(i > n_e\). However, as \(i \rightarrow n\), \(f_{\beta}(\varepsilon_i) \rightarrow 0\).
The Fermi-level \(\mu\) is chosen to ensure \(\int \rho(r) dr = 2n_e\).
Note that Fermi-Dirac distribution is not the only way to regularize the Kohn-Sham problem. The regularization does change the problem to be solved. However, since the original Kohn-Sham model is an approximation, such a modification is justified as long as the modified model is effective in applications such as geometry optimization and ab initio molecular dynamics.
The regularization introduces additional computation to be performed in the SCF iteration:
- We need to compute more than \(n_e\) eigenpairs in the diagonalization procedure in order to evaluate the electron density as
for some \(k>n_e\).
- We need to determine \(\mu\) by solving \(\int \rho(r) dr = 2n_e\). Because the left hand side is a monotonic function with respect to \(\mu\), \(\mu\) can be determined by a bisection procedure.
In KSSOLV2.0, the regularized Kohn-Sham equation can be solved by the function scf4m. The interaface for this function is identical to that of scf. To perform a finite temperature calculation for a regularized model,
- we need to define the
temperatureattribute (in Kelvin) in theMoleculeorCrystalinput object. For example,
>> cry = Crystal('supercell',C, 'atomlist',atomlist, 'xyzlist' ,xyzlist, ...
'ecutwfc',30, 'name','Cu', 'autokpts', autokpts, 'temperature', 300 );
- In addition, we need to initialize the
Wavefuninput object to have more than \(n_e\) wavefunctions.
Note that we keep both scf.m and scf4m.m so that the difference between these two solvers can be easily compared, even though it is possible to use scf4m.m to perform a zero temperature calculation.