Total energy optimisation

The total energy is a functional of the charge density: $E = E[n(\mathbf{r})]$. The charge density itself is expanded in terms of the basis $\{ B(\mathbf{r}) \}$ and depends upon the density kernel elements, $\{ K^{\alpha\beta} \}$, and the NGWF expansion coefficients, $\{ C_{KLM,\alpha} \}$. Provided that it is $N$-representable, this dependence should be variational, i.e., the ground state energy, $E_{\mathrm{min}}$, is given by

$$E_{\mathrm{min}} = \min_{\{K^{\alpha\beta} \},\{ C_{KLM,\alpha} \}} E( \{K^{\alpha\beta}\}, \{C_{KLM,\alpha}\} ).$$ (39)

In this work we are concerned principally with the optimisation of the elements of the density kernel and we will consider the NGWF coefficients, $\{ C_{KLM,\alpha} \}$, as being fixed. We use the pseudo-atomic orbitals (PAOs) of Sankey et al. [20] as our NGWFs.

This minimisation must be performed under the constraints of constant electron number,

$$N_{\mathrm{e}} = \int_{V} d\mathbf{r} \: n (\mathbf{r}) = 2K^{\alpha\beta}S_{\beta\alpha} = 2\mathrm{Tr}[\mathbf{KS}],$$ (40)

and density-matrix idempotency,

$$\rho (\mathbf{r},\mathbf{r}') = \int_{V} d\mathbf{r}'' \: \r... ...lpha\beta} = K^{\alpha\gamma}S_{\gamma\delta} K^{\delta\beta},$$ (41)

where the overlap matrix, $S_{\alpha\beta}$ is given by

$$S_{\alpha\beta} = \int_{V} d\mathbf{r} \: \phi_{\alpha}(\mathbf{r})\phi_{\beta}(\mathbf{r}).$$ (42)

In order to avoid explicitly imposing the idempotency constraint (41), we use the method suggested by Li, Nunes and Vanderbilt [21] and independently by Daw [22], and generalised to the case of non-orthogonal functions by Nunes and Vanderbilt [23]. Our implementation follows the simplified version of Millam and Scuseria [24]. We define the following function of the density kernel $\mathbf{K}$:

$$L(\mathbf{K}) = E(\tilde{\mathbf{K}}) - \mu (2\mathrm{Tr}[\mathbf{KS}] - N_{\mathrm{e}}),$$ (43)

where $\tilde{\mathbf{K}}$ is the McWeeny purified density kernel [25],

$$\tilde{\mathbf{K}} = 3\mathbf{KSK} - 2\mathbf{KSKSK}.$$ (44)

The contravariant, tensor-corrected gradient [15,16] that is used in the steepest descent or conjugate gradient iterative minimisation is given by

$$\nabla L = \mathbf{S}^{-1}\frac{\partial L}{\partial \mathbf{K}} \mathbf{S}^{-1}$$

$$= 6\mathbf{KHS}^{-1} + 6\mathbf{S}^{-1}\mathbf{HK} - 4\mathbf{S}^{-1}\mathbf{HKSK}$$

$$- \mbox{ } 4\mathbf{KHK} -4\mathbf{KSKHS}^{-1} - 2\mu \mathbf{S}^{-1}.$$ (45)

The value of $\mu$ is set at each step such that $\mathrm{Tr}[\mathbf{S} \nabla L] = 0$. This ensures that the total electron number remains unchanged, thus we simply require that our initial guess for the density kernel gives the correct number of electrons.