4.3 Wavefunctions

4.3.1 Form of trial wavefunction

The trial wavefunctions used in this thesis all take the form of a sum of products of Slater determinants multiplied by a correlating Jastrow factor:

$$\Psi({\bf R})=\sum_{n} a_{n} D^{\uparrow}_{n}({\bf r}^{\upar... ...ldots,{\bf r}^{^\downarrow}_{N}) \exp \left[ J \right] \;\;\;.$$ (4.2)

The general trial wavefunction, $\Psi$, consists of a sum of products of up- and down-spin determinants. Separate determinants are used for the $N^\uparrow$ up-spin and $N^\downarrow$ down-spin electrons, consisting of sets of single-particle orbitals:

$$D_{n}^{\uparrow}=\left\vert \begin{array}{cccc} \phi_1({\bf... ...row}({\bf r}_{N^\uparrow})\\ \end{array} \right\vert \;\;\; .$$ (4.3)

The use of separate determinants for differing spins results in a trial function that is not antisymmetric with respect to interchange of opposite spin electrons, but gives the correct expectation values for spin-independent operators. [23] The Jastrow factor, $J$, in the most general case, is a function of all electron and atom positions. Jastrow functions are discussed in the next section.

The determinantal part of the wavefunction is usually obtained from a self-consistent method such as LDA-DFT or HF. [37,38] For most of the calculations in this thesis a single product of up and down-spin determinants of orbitals, $\phi$, is used to form the wavefunctions. In principle, configuration-interaction or multi-configuration calculations can be used to provide multiple determinants, including the coefficients $a_{n}$,^4.3 but many applications of QMC have demonstrated that a high accuracy may usually be obtained using only a single determinant.

The orbitals, $\phi$, are expanded in a computationally convenient basis-set. In principle the orbitals could be expanded in a basis that is optimal for QMC. The key criterion would be the speed of evaluation of the value and first and second derivatives of the orbitals at arbitrary points in space. However, to be obtainable from within a self-consistent scheme the basis must be convenient for solving the single-particle equations. Consequently, radial numerical grids, Gaussian functions [39] and plane-waves [40] are typically used for atomic, molecular and continuum systems respectively.^4.4

A further consideration is that the orbitals in the Slater determinants must be purely real functions. In supercell (periodic boundary condition) calculations, this is achieved by choosing sets of bloch orbitals such that simple linear combinations of the orbitals yield real functions.

4.3.2 Form of Jastrow factor

The Jastrow factor [22] is the primary method for including electron correlation beyond the mean-field level in QMC simulations. In VMC it is the only method. Although a Jastrow factor is not in principle required for DMC calculations, a HF wavefunction is usually insufficiently accurate for use as a guiding wavefunction. A good quality Jastrow factor should incorporate the physics of electron correlation in a compact and rapidly computable form that is convenient for optimisation. The Jastrow factor is therefore one of the keys to efficient QMC simulations.

The number of forms of Jastrow factor that have enjoyed repeated use in electronic structure applications is limited. In this section, three forms are presented, two of which are used in subsequent calculations. A thorough study of the efficiency of different forms would be difficult to produce due to the computational cost involved. Alexander and Coldwell [41] reviewed 118 trial wavefunction forms for the He, Li and Be atoms, but did not consider other species, or molecular and continuum systems.

Although the quality of a Jastrow function can be judged by the VMC energies obtained for candidate systems, compromises must inevitably be made. Benchmark accuracy Jastrow factors for molecular and continuum systems typically consist of very large numbers of terms and parameters. The increase in accuracy per additional term (of a given kind) inevitably reduces. Consequently, a wavefunction obtaining $80\%$ of the correlation energy and 30 parameters may be preferable to one obtaining $90\%$ of the correlation energy with 200 parameters. For use as a guiding DMC wavefunction, the simpler wavefunction would be preferable due to the lower computational cost. For VMC calculations, the preference would depend on the specific application and desired final accuracy.

4.3.2.1 Mitás' approach

Mitás et al. have successfully optimised wavefunctions for atomic systems, [42] molecular carbon [43,44] and silicon [45] systems, using a Jastrow function of the following form:

$$J=\sum_{I} \sum_{i}^{N} \sum_{j\neq i}^{N} u(r_{iI},r_{jI},r{ij}) \;\;\;,$$ (4.4)

$$u(r_{iI},r_{jI},r_{ij})=-\frac{c}{\gamma}e^{-\gamma r_{ij}} ... ...r_{iI})a_l(r_{jI})+a_k(r_{jI})a_l(r_{iI})]b_m(r_{ij}) \;\;\; ,$$ (4.5)

where

$$a_k(r) = \left(\frac{\alpha_kr}{1+\alpha_kr}\right)^2 \;,\;\alpha_k=\frac{\alpha_0}{2^{k-1}}\;,\;k>0$$ (4.6)

$$b_m(r) = \left(\frac{\beta_mr}{1+\beta_mr}\right)^2 \;,\;\beta_m=\frac{\beta_0}{2^{m-1}}\;,\;m>0 \;\;\; ,$$ (4.7)

and $r_{iI}$ is the separation between the $i^{th}$ of $N$ electrons and the $I^{th}$ ion, and $a_0(r)=b_0(r)=1$. This form can satisfy the electron-electron cusp condition and has no electron-nucleus cusp. This is desired as the orbitals in the determinant satisfy, or can be made to satisfy, the electron-nuclear cusp.

This functional form has yielded wavefunctions of sufficient accuracy to perform DMC calculations in several systems, but has not been systematically tested for accuracy and convergence of VMC energies.

Although successfully used in several calculations, the functional form has two undesirable features. The correlation function, $u$, is long ranged implying that the sum over electrons and ions is not readily cutoff and will increase in cost with system size. For certain values of $\gamma$, the long range of the exponential term requires the use of Ewald summation (section 4.6.1) in supercell calculations. As this sum involves all pairs of electrons, this operation would become costly.

Additionally, the function is not linear in the $\alpha_k$ parameters, which is desirable to avoid full recalculation of the function as the parameters are changed during optimisation. This increases the computational cost of the optimisation (see section 4.7.1).

4.3.2.2 Williamson's approach

Williamson et al. [46,47] developed a homogeneous Jastrow function designed for solid-state periodic boundary condition applications, which avoids potentially costly sums. This was used in conjunction with a one-body ``$\chi$'' term that had been used in Fahy's early bulk calculations. [26]

$$J=\sum_{i=1}^{N}\chi({\bf r}_{i}) -\sum_{i=1}^{N}\sum_{j\neq i}^{N}u({r}_{ij}) \;\;\;,$$ (4.8)

where ${\bf r}_{i}$ is the absolute position of electron $i$, and ${r}_{ij}$ is a minimum-image [17] electron-electron distance. The one-body term consists of a plane-wave expansion

$$\chi({\bf r}_{i}) = \sum_{{\bf G}}\chi({\bf G})e^{i{\bf G}.{\bf r}_{i}} \;\;\;,$$ (4.9)

where reciprocal lattice vectors are denoted by ${\bf G}$. The number of independent coefficients $\chi({\bf G})$ is reduced by imposing symmetry.

The electron-electron dependent part of the Jastrow function consists of cusp-satisfying part and a separate power expansion, both of which are constrained to be short ranged. For the cusp-satisfying part, the following form was chosen:

$$u_0(r)=\frac{A}{r}\left(1-\exp (-\frac{r}{F})\right)\exp\left(-\frac{r^2}{L_{0}^{2}}\right) \;\;\;,$$ (4.10)

where $F$ is chosen so that the cusp condition is obeyed and $L_0$ is chosen so that $u_0$ is effectively zero at the Wigner-Seitz cell edge in supercell calculations. A is a variational parameter, and $L_0$ is fixed at approximately $30\%$ of the Wigner-Seitz cell radius. A power expansion in Chebyshev polynomials, over a pre-determined range $L$, was used in a second term:

$$\begin{array}{llll} f(r) & = & B(\frac{L}2+r)(L-r)^2\;\;\;\;... ...M} c_l T_l(\overline{r}) \\ & = & 0 & r>L \end{array}\;\;\; ,$$ (4.11)

where $B$ and $c_{l}$ are variational coefficients and $T_l$ denotes the $l$th Chebyshev polynomial. The expansion range, $(0,L)$, of the correlation function is rescaled to the orthogonality interval of the Chebyshev polynomials, $[-1,1]$,

$$\overline{r}=\frac{2r-L}{L} \;\;\;.$$ (4.12)

The parameters $B$, $c_{l}$ are different for spin parallel and anti-parallel electron pairs, to model the additional Pauli repulsion experienced by parallel spin electrons. These functions satisfy the electron-electron cusp condition and are continuous in value and first derivative. The expansion range, $L$, is chosen to be less than or equal to the Wigner-Seitz radius in supercell calculations avoiding the need to perform an Ewald sum. A further advantage of the Chebyshev expansion is that unlike the terms in a simple power expansion, the terms carry equal weight over their range and consequently exhibit greater stability during optimisation.

In applications to bulk carbon and silicon, wavefunctions obtaining $80-85\%$ of the fixed-node DMC correlation energy were obtained. The high accuracy of these wavefunctions indicates that the limited range (form) of the correlation functions does not severely limit the accuracies obtained.

4.3.2.3 Atom centred approach

For the calculations in this thesis, a development of the above approach has been used. The plane-wave $\chi$ function was replaced by atom-centred functions, expanded over a short range in a complete polynomial. By extending the idea of short ranged power expansions to atom centred functions, partially transferable correlation functions result. This is distinct from the plane-wave $\chi$ function which is not easily transfered to other systems. The high accuracy of the previous functions were further increased by adding parameterised terms dependent on electron-electron-ion distances. These are present in Mitás' function and are important in all-electron calculations and systems with tightly bound electrons.

A general power expansion, in $r_i$ and $r_{ij}$, about an atom over a range $L$ must satisfy similar conditions to Williamson's function, above: the function must go to zero at $L$ and be continuous in value and first derivative. The atom centred power expansion in electron-ion distance, $r_i$ and electron-electron distance $r_{ij}$ contains 5 independent series: [48,49]

$$S_{1} = (r_{i}-L)^{2}r_{i}^{2}\sum_{l\ge 0}\alpha_{l}T_{l}(\bar{r}_{i})+B(r_{i}-L)^{2}\left(\frac{L}{2}+r_{i}\right)$$ (4.13)

$$S_{2} = \left(r_{ij}-2L\right)^{2}r_{ij}^{2}\sum_{n\ge0}\beta_{n}T_{n}(\bar{r}_{ij})+B^{\prime}\left(r_{ij}-2L\right)^{2}\left(\frac{2L}{2}+r_{ij}\right)$$

$$S_{3} = (r_{i}-L)^{2}(r_{j}-L)^{2}r_{i}^{2}r_{j}^{2}\sum_{l\ge0}\sum_{m\ge 0}\gamma_{lm}T_{l}(\bar{r}_{i})T_{m}(\bar{r}_{j})$$

$$S_{4} = (r_{i}-L)^{2}\left(r_{ij}-2L\right)^{2}r_{i}^{2}r_{ij}^{2}\sum_{l\ge0}\sum_{n\ge 0}\epsilon_{ln}T_{l}(\bar{r}_{i})T_{n}(\bar{r}_{ij})$$

$$S_{5} = (r_{i}-L)^{2}(r_{j}-L)^{2}r_{i}^{2}r_{j}^{2}r_{ij}^{2}\sum_{l\ge ... ...e 0}\omega_{lmn}T_{l}(\bar{r}_{i})T_{m}(\bar{r}_{j})T_{n}(\bar{r}_{ij}) \;\;\;,$$

where it is implicitly assumed that terms are zero for $r_{i}>L$ or $r_{ij}>2L$ and

$$\overline{r}_{ij}=\frac{r_{ij}-L}{L} \;\;\;.$$ (4.14)

The overall Jastrow factor consists of sets of the above functions supplemented by the electron-electron cusp-satisfying $u_0$ term, equation 4.10, and the homogeneous $r_{ij}$-dependent $f(r_{ij})$ term, equation 4.11. This is shown schematically in figure 4.1.

Figure 4.1: Example Jastrow factor parameterisation for $\mathrm{CO}_3$ ( $\mathrm{C}_{2v}$ symmetry). Overlapping atom-centred functions are placed on the carbon (filled circle) and oxygen (empty circle) atoms and describe electron-ion, electron-electron and electron-electron-ion correlations within a specified range (shading). The functions for the oxygen atoms are identical by symmetry. An additional longer-ranged electron-electron term (not shown) describes $r_{ij}$-dependent correlations near atoms and also in the valence regions.

For the main applications in this thesis, only the terms $S_{1}$ were used as sufficient accuracy was obtained and the full series had yet to be implemented. In tests on the ground state of an all-electron oxygen atom, using exact HF orbitals, $89\%$ of the experimental correlation energy was obtained using all the terms, approximately $20\%$ more than a wavefunction using only the $S_{1}$ terms. This accuracy is similar to that obtained in pseudopotential calculations of oxygen atoms, demonstrating the importance of the additional correlation terms when treating core electrons.