6.4 Finite size correction and extrapolation formulae

The basic idea behind finite size correction formulae is to write the energy for the infinite system as

$$E_{\infty} = E_{N} + (E_{\infty} - E_{N})\;,$$ (6.2)

where the subscript denotes the system size. A highly accurate many-body calculation for the $N$-particle system is then performed to obtain an approximation for $E_{N}$, and the correction term in brackets is approximated using a much less expensive scheme which can be applied to very large systems:

$$E_{\infty} \simeq E_{N} + (E^{\prime}_{\infty} - E^{\prime}_{N})\;,$$ (6.3)

where the prime indicates that a less expensive scheme is used. In an extrapolation procedure the energy $E_{N}$ is calculated for a range of system sizes and fitted to a chosen functional form which contains parameters. Correction and extrapolation procedures can be combined to give an expression

$$E_{\infty} \simeq E_{N} + (E^{\prime}_{\infty} - E^{\prime}_{N}) + F(N)\;,$$ (6.4)

where the prime again indicates that a less expensive scheme is used and $F(N)$ is an extrapolation function. Clearly, the optimal form of $F(N)$ depends on the method used to calculate the correction term. The practical difference between correction and extrapolation is that correction requires a single calculation of $E_{N}$ using the accurate (and normally very costly) many-body technique, while extrapolation requires several such calculations for different values of $N$ and a subsequent fit to the chosen functional form $F(N)$. The extrapolation procedure is costly because it involves more calculations and is prone to inaccuracy because one has to perform a fit with only a few data points. In designing a correction/extrapolation procedure one therefore tries to make the extrapolation term as small as possible.

Candidate methods for evaluating the correction term include Kohn-Sham DFT and HF theory. The most convenient methods are the local density approximation (LDA) to DFT, or extensions such as the Generalized Gradient Approximation. These methods are very widely applied in periodic boundary conditions calculations and are computationally inexpensive, while retaining a realistic description of the system. Within an independent particle theory, such as Kohn-Sham DFT, calculations for periodic systems are normally performed by solving within the primitive unit cell.^6.1 To obtain the correct result for the infinite system it is necessary to integrate over k-space, and the integral is normally approximated by a sum over a finite set of k-points. A determinant formed from the occupied orbitals at a single k-point in the first Brillouin zone (BZ) of the primitive unit cell is a many-body wave function for a simulation cell of the size of the primitive unit cell. Adding a second k-point doubles the size of the determinant and is equivalent to doubling the size of the simulation cell, etc.

An important subtlety arises when finite size corrections derived from an independent particle theory such as LDA-DFT are used to correct the results of a true many-body method such as QMC. Suppose the many-body calculation is performed using the Ewald interaction, so that the Coulomb energy is that of an infinite periodic array of copies of the simulation cell. Given that the solid we are trying to model is crystalline, and hence the charge density is truly periodic, the Ewald interaction gives a good description of the classical Coulomb or Hartree energy. However, because the electronic positions are mirrored exactly in every copy of the simulation cell, the electronic correlations are also forced to be periodic, and the exchange-correlation (XC) energy corresponds to a system with a periodic XC hole. This unphysical approximation is particularly inaccurate when the simulation cell is small. In Kohn-Sham DFT calculations the XC energy is evaluated using a standard functional such as the LDA of Perdew and Zunger [88]. This was obtained by fitting to the results of DMC calculations for jellium, [13] which were extrapolated to the infinite-system-size limit prior to performing the fit.

The consequences of this difference may be illustrated by considering many-body and LDA calculations for jellium, in which the charge density is uniform. The LDA gives the XC energy for the infinite system irrespective of the size of the simulation cell, but the XC energy obtained from a many-body simulation using the Hamiltonian of equation 6.1 with the Ewald interaction gives an XC energy which depends on the size and shape of the simulation cell. Consequently, evaluating the correction term in equation 6.4 within the LDA does not give a good approximation to the finite size correction in the many-body simulation and a significant extrapolation term remains.

6.4.1 Previous work

Several methods for computing finite size corrections in many-body calculations occur in the literature.

The issue of finite size corrections to both the kinetic and interaction energies has been addressed by Ceperley and coworkers. [89,90,91,33] Their approach is to add separate extrapolation terms for the kinetic and interaction energies. In their work on hydrogen solids, Ceperley and Alder [90] performed DMC calculations for a number of different system sizes and fitted to the formula

$$E^{\rm DMC}_{\infty} \simeq E^{\rm DMC}_N + a(T_{\infty} - T_{N}) + \frac{b}{N}\;\;,$$ (6.5)

where $a$ and $b$ are parameters, and $T$ is the kinetic energy of the non-interacting electron gas. The $b/N$ term accounts for the finite size effects arising from the interaction energy and the difference of the parameter $a$ from unity accounts for the difference between the kinetic energies of the interacting and non-interacting systems. Normally $a>1$ because the interacting kinetic energy is larger than the non-interacting kinetic energy. Engel et al. [92] used the following formula for inhomogeneous systems:

$$E^{\rm QMC}_{\infty} \simeq E^{\rm QMC}_N + a(E^{\rm LDA}_{\infty} - E^{\rm LDA}_{N}) + \frac{b}{N}\;\;,$$ (6.6)

which reduces to equation 6.5 for a homogeneous system.