# Analytic nuclear forces and molecular properties from full configuration interaction quantum Monte Carlo

## Abstract

Unbiased stochastic sampling of the one- and two-body reduced density matrices is achieved in full configuration interaction quantum Monte Carlo with the introduction of a second, “replica” ensemble of walkers, whose population evolves in imaginary time independently from the first, and which entails only modest additional computational overheads. The matrices obtained from this approach are shown to be representative of full configuration-interaction quality, and hence provide a realistic opportunity to achieve high-quality results for a range of properties whose operators do not necessarily commute with the hamiltonian. A density-matrix formulated quasi-variational energy estimator having been already proposed and investigated, the present work extends the scope of the theory to take in studies of analytic nuclear forces, molecular dipole moments and polarisabilities, with extensive comparison to exact results where possible. These new results confirm the suitability of the sampling technique and, where sufficiently large basis sets are available, achieve close agreement with experimental values, expanding the scope of the method to new areas of investigation.

## 1Introduction

The full configuration interaction quantum Monte Carlo method (FCIQMC) and its initiator adaptation (-FCIQMC) are projector QMC techniques, capable of providing near-exact, systematically improvable descriptions of correlated wavefunctions expressed as linear combinations of Slater determinants.[1] This convergence is achieved by stochastically sampling the exponentially large (though finite) Hilbert spaces of configuration interaction theory *via* a population dynamics performed on an ensemble of signed walkers. Annihilation processes provide a means of combating the ill effects of the fermion sign problem which plagues projector approaches,[3] exploiting the sparsity of the wavefunction induced by a judicious choice of orbital basis. The approach requires substantially less computational effort than an iterative diagonalisation technique, and has thus found considerable success in studies of atomic and molecular systems,[6] model systems such as the homogeneous electron gas and the Hubbard model,[13] and solid-state systems.[16]

The principal focus of many of these studies has been to derive properties based upon total energies, for which an unbiased projected estimator is readily available, and which have included excitation and dissociation energies,[8] electron affinities,[7] ionisation potentials,[6] and equations of state.[16] Despite their success, however, the extension to include the calculation of a greater range of properties — expectation values of operators which do not necessarily commute with the hamiltonian — remains highly desirable. This focus has been the subject of considerable interest for the QMC community in general, and has posed a considerable challenge for decades.[17] These properties, which include static correlation functions and entropy estimators as well as the forces, multipole moments, and polarisabilities considered here, may be deduced from the effect of a perturbation from the corresponding operator, , upon the hamiltonian,

with the perturbation strength, such that the expectation value, , is given by the derivative of the energy with respect to , evaluated at :

In accordance with the Hellmann–Feynman theorem,[30] applicable to converged (normalised) -FCIQMC wavefunctions by analogy with deterministic and strictly variational FCI, this expression reduces to

or equivalently to the trace of with the appropriate rank of reduced density matrix.[31] It is worth noting that unconverged -FCIQMC wavefunctions need not rigorously obey the Hellmann–Feynman theorem, and so in this work we ensure that we are working in the large walker limit, such that systematic errors in the sampled distribution due to insufficient walker numbers have been minimized to the FCI-limit.

The effective stochastic acquisition of these reduced density matrices, therefore, has the capacity to broaden the scope of -FCIQMC significantly, and motivates the present work. We begin with a brief overview of the -FCIQMC algorithm, including its extension to non-integer walker weights,[32] before recapitulating some of the details of the “replica” density-matrix sampling technique.[33] Building upon that previous work, our discussion turns to consider the calculation of nuclear forces, molecular dipole moments, and atomic dipole polarisabilities, and in so doing confirms the high quality of the sampled one- and two-body reduced density matrices which is now achievable.

## 2Methodology

### 2.1-Fciqmc

Initiator full configuration interaction quantum Monte Carlo provides stochastic integration of the -electron, imaginary-time Schrödinger equation, yielding wavefunctions expressed as a linear combination of the set of Slater determinants, , formed from the underlying one-particle (most often Hartree–Fock) basis:

The coefficients of this wavefunction expansion are obtained by iterative application of the equations

representing the evolution of the coefficients over a timestep in imaginary time. This evolution is achieved by subjecting an ensemble of signed walkers to a three-step population dynamics algorithm of “spawning”, “death”, and “annihilation” steps, the walker populations, , becoming proportional to the coefficients. The full details of this approach have been expounded in previous papers,[1], and what follows should be regarded only as a brief summary.

Typically initialised with a single walker placed upon the Hartree–Fock determinant, , a simulation using integer walkers proceeds with a coupled determinant, , being randomly selected for each walker on parent determinant, , with a probability . The determinant selected, the parent walker then attempts to spawn a child on to it, with a probability

If the attempt is successful, the sign of the spawned walker matches that of its parent if and is inverted if . The initiator adaptation, -FCIQMC, modifies the spawning step by introducing a parameter, , which specifies a lower population threshold under which a parent determinant is prevented from spawning on to unoccupied determinants. Each walker next attempts to die, with a probability given by

in which is a population control parameter — known as the “shift” — which tends to the ground-state energy in the long- limit.

These two steps are themselves sufficient to describe Equation 2 fully, but are insufficient to provide convergence to a fermionic wavefunction. Instead, a third step — “annihilation” — is required in order to suppress the deleterious effects of the fermion sign problem.[3] After each iteration, walkers of opposite sign on the same determinant annihilate, and in so doing ensure that each determinant is populated by walkers of only one sign for the next iteration. The success of these processes relies on the sparsity of the wavefunction induced by the underlying basis — typically chosen to be Hartree–Fock orbitals — which confines it to a generally small region of the Hilbert space. In so doing, it ensures that annihilation events are numerous enough to maintain the sign structure of the sampled wavefunction accurately.

Although the walkers of ()-FCIQMC were initially conceived as an ensemble of discrete particles, there is some merit in instead positing a set of non-integer walkers.[32] Such an approach reduces the amount of random number generation required, reduces the instantaneous fluctuations in the populations on a given determinant, and hence the fluctuations in the energy estimators in imaginary time.

This formulation is achieved by applying the spawning, death, and annihilation steps introduced earlier *continuously*, rather than discretely. Thus, instead of spawning a walker of signed integer weight from a determinant to a coupled determinant with a probability , a walker of weight is spawned with probability . Likewise, the death step is remodelled such that it simply involves reducing the population on a determinant by . Annihilation is achieved by taking the signed sum of walkers on each determinant on a given iteration as the residual population for the next iteration. For -FCIQMC calculations, the parameter is recast as a continuous variable rather than an integer.

The continuous nature of the spawned walkers in this approach does not, however, imply that the number of spawning *events* becomes continuous. As in the integer formulation, where there are exactly spawning attempts from determinant with a population on each iteration, there are a discrete number of attempts per determinant per iteration. For practical purposes, a continuous spawning threshold, , is introduced such that if , walkers are spawned with a probability . This implementation is designed to alleviate the significant cost of low-weighted spawnings compared to their effect on the overall wavefunction, as well as ensuring that the wavefunction remains compact and expressible instantaneously by a number of walkers far smaller than the size of the space.

Whilst the death step requires no extra modification of this kind, some additional considerations must be addressed for annihilation. In order that determinants can become completely depopulated, and we are not forced to store large numbers whose populations are very close to, but not exactly, zero, a minimum occupation threshold, , is imposed upon them. If, after annihilation, the population on a determinant , its population is set either to with probability , or else to with probability .

As a final practical means of alleviating the computational burden of this approach, it is possible to treat only a subspace of the full Hilbert space with non-integer walkers, continuing to describe the remainder in a discretised fashion. In order to preserve the benefits of the non-integer approach on the fluctuations of the energy estimators, the truncation is specified by an excitation level, , with only -fold and lower excitations from the reference included in the non-integer subspace. A typical choice of parameters , ( is used here), and entails only a modest increase in the computational cost of the calculation over the integer implementation, while retaining many of the benefits of the full non-integer approach.

-FCIQMC provides two essentially independent energy estimators, which, taken together, provide a useful confirmation of the validity of the obtained result. The first, to which we have already alluded, is the shift, . This is initially held constant (typically at zero) to facilitate an exponential growth in the number of walkers, before being allowed to vary dynamically to keep the population constant. At convergence, this variation fluctuates around the energy of system, and thus provides an energy on the basis of the growth rate of the entire ensemble of walkers. A projected energy estimator, of the form

on the other hand, depends only upon the populations of the determinants coupled to the reference state, . Whilst the error in this projection is formally first-order in the wavefunction error, its non-variationality tends to mean that it converges rather faster to the exact, infinite-walker limit than does a variational estimator, owing to favourable cancellation of errors. The projected energy is thus typically preferred when the wavefunction is dominated by the Hartree–Fock determinant, but a projection on to a multi-reference trial wavefunction or the variational estimator provided by the density matrices (which is second-order in the wavefunction error) are often more useful in more strongly-correlated cases.[29] Once the ensemble has equilibrated, the simulation is allowed to evolve in imaginary time until the statistical errors in both and have been satisfactorily reduced, upon which a Flyvbjerg–Petersen blocking analysis is performed to estimate the error in the obtained result.[37]

### 2.2Stochastic density-matrix sampling

In terms of the wavefunction ansatz of -FCIQMC (Eq. Equation 1) and the creation and annihilation operators, the one- and two-body reduced density matrices, and , may be formulated in terms of the wavefunction expansion and the conventional creation and annihilation operators, and , as

and,

respectively,[38] and an important recent development of the theory allows these objects to be sampled in an efficient, stochastically unbiased fashion.[40]

The diagonal elements of these objects, of the form

may be calculated straightforwardly, as each determinant, , contributes to each of the matrix elements involving its occupied orbitals. The corresponding explicit generation of all the required determinant pairs for the off-diagonal elements is not practical, but the observation that the relevant pairs are at most double excitations of one another allows both and to be sampled *via* the spawning steps.[40] Thus, the existing computational effort required for the communication of the spawning event need only be slightly accentuated (by the need now to convey both the amplitude and the identity of the parent determinant to the child) to allow the contributions to the off-diagonal matrix elements from determinant pairs to be calculated on the fly.

As these off-diagonal contributions are only added upon a *successful* spawning event, it is necessary that they be rescaled according to the probability of such an event taking place. That is, a contribution will instead be accumulated as

with the probability that at least one spawning attempt from to is successful on a given iteration. Depending on whether integer or non-integer walkers are considered, this probability is given by

with the instantaneous walker population residing on and the probability that no walker is spawned between and in a single attempt. For an integer spawning event, this probability is

but this must be modified in the case of continuous spawning to

where is the continuous spawning threshold, if used.

A naïve implementation of the above sampling is satisfactory for the accumulation of approximate density matrices, but is beset by a number of shortcomings which should be considered.[29] As contributions to the off-diagonal matrix elements are only added upon a successful spawning attempt, problems can arise when the spawning events are discretised. In this case, the probability that such an event occurs is proportional to the coupling hamiltonian matrix element, , and pairs of determinants which are connected by large matrix elements are correspondingly sampled more often than pairs which are only weakly coupled. Thus, if two highly-weighted determinants contained in the stochastically-sampled, integer walker space are connected by a small hamiltonian element, their contribution to the density matrices may be severely under-represented, or neglected entirely.

This problem is most notably in evidence in the case of single excitations of the reference determinant, for which the coupling matrix elements are strictly zero according the Brillouin’s theorem. This is countered in the present implementation by accounting for these contributions to the density matrices explicitly, and hence removing the dependence upon a successful spawning event. Other contributions, however, whose sampling will still be proportional to the reduced hamiltonian,[41] defined in terms of the one- and two-electron integrals, and , as

will give rise to an biasing error in density matrices in the long- limit for determinant pairs where , but whose amplitudes are both significantly non-zero. However, modifications to the algorithm to treat the bias remaining beyond that already defined by Brillouin’s theorem explicitly — such as introducing additional events to spawn walkers proportionally to the inverse of the hamiltonian element — have been shown to be of little additional benefit due to the negligible nature of this bias in numerical studies to date.[29]

In a separate difficulty, it has been shown previously that a straightforward implementation of the above sampling gave rise to a convergence of the density matrices with increasing which was rather slower than that of, say, the projected energy. This behaviour stems not simply from undersampling, but rather from a bias in the statistical sampling technique itself. In particular, appropriate contributions to the matrix elements are approximated by

ignoring the potentially significant covariance, , between the two amplitudes and introducing a bias, whether or not the averaged walker populations are themselves unbiased. It is, to that end, unsurprising that this problem is at its greatest for diagonal elements, for which the “two” amplitudes are perfectly correlated, and — the error being of a single sign — there is no possibility of error cancellation.

This problem is rather more serious than the previous concerns over discretised spawning, but one for which a rather simple solution exists. Unbiased density matrices can be calculated with the introduction of a second, uncorrelated walker ensemble, to which the stochastic spawning, death, and annihilation steps are applied independently, and whose statistics are acquired separately, from the first.[29] This adaptation, known as replica sampling, achieves the unbiasing by ensuring that all the products of determinant amplitudes are calculated using populations from both simulations, and has previously found application in the stochastic sampling of the -electron density matrix known as density matrix quantum Monte Carlo,[42] and the recently-introduced Krylov-projected quantum Monte Carlo.[34] That is, for example, a successful spawning event from to in replica , occurring with a probability , gives rise to a contribution of:

This approach bears some conceptual similarity with the bilinear sampling algorithm in Green’s function Monte Carlo, introduced by Zhang and Kalos, in that both seek a means of finding expectation values of operators which do not commute with the hamiltonian, *via* two sets of independent walker distributions.[33] The main difference, though, is that the bilinear approach transforms the Schrödinger equation such that there are two related wavefunctions to sample, while in the present work the walker ensembles are independent samples of the same underlying object. In providing a stochastically unbiased route to the density matrices, the replica sampling technique thus provides the first realistic opportunity to achieve high-accuracy *ab initio* results for the sizeable suite of properties that can be derived therefrom.

## 3Nuclear forces

The force acting on a nucleus in a molecule or cluster is defined as the negative gradient of the molecular energy with respect to the nuclear coordinates:

In Eq. the symbol denotes the nuclear force vector, the energy of a molecule at a fixed geometry in the electronic ground state, and refers to the vector of nuclear coordinates in the centre of mass frame of reference. A comprehensive review of techniques and explicit expressions to compute derivatives of the electronic energy with respect to nuclear coordinates is available in the literature.[43] The following discussion is thus limited to the basic concepts for the calculation of nuclear forces using all-electron FCI wavefunctions as obtained as a statistical average using the -FCIQMC method, once the calculation has been converged with respect to the number of walkers. In the present work, we have adjusted the total number of walkers to achieve a walker population of 50000 at the reference (i. e. highest populated) determinant. Preceding work confirmed that, at such a population levels, noise arising from small stochastic populations of random determinants is sufficiently suppressed and the wavefunction converged.

The first derivative of the electronic energy of a CI wavefunction generally depends on the derivatives of the atomic orbitals (AOs) and the molecular-orbital (MO) and CI coefficients. All these terms depend upon the nuclear coordinates, and the computation of nuclear forces requires knowledge of the first derivatives with respect to all considered degrees of freedom. However, electronic wave functions obtained from -FCIQMC optimizations are variational with respect to the CI coefficients and a component of the nuclear force vector can be expressed in terms of the reduced density matrices as

in which the terms represent the one-electron integrals from the hamiltonian, and is the contribution from the fixed nuclei. Moreover, all-electron FCI wavefunctions considered in this work are also invariant under variation of the MO coefficients. The nuclear forces can thus be expressed solely in terms of the one- and two-electron density matrices and the skeleton derivative integrals of the basis functions:

where

is an element of the lagrangian and is an element of the overlap matrix. In particular, neither the computation of derivatives of the CI hamiltonian matrix nor the solution of the coupled-perturbed Hartree–Fock equation for the derivatives of the MO coefficients are required. We have implemented an interface to `MOLPRO`

to compute the integrals and nuclear forces from the -FCIQMC density matrices.[46]

As a first benchmark, we have applied the -FCIQMC methodology to compute the nuclear forces at several points along the dissociation curve of molecular nitrogen, as the electronic wavefunction changes from single- to strong multi-reference character.

Figure 1 (top) compares the potential energy computed with -FCIQMC and the FCI program in `MOLPRO`

using a small 6-31G basis set to allow for comparison to exact (FCI) results. The accuracy of the -FCIQMC methodology for the computation of total energies was already evaluated [2], and we generally find excellent agreement between the -FCIQMC and FCI data set.

In Figure 1 (bottom), the nuclear forces for the same geometries are illustrated. Comparison with FCI results obtained from numerical gradients provides a direct measure of the quality of the reduced density matrices computed from the replica algorithm based on -FCIQMC, and, once again, the data shows excellent agreement between the analytic -FCIQMC forces and the FCI results for all geometries.

As second example for the calculation of analytic gradients and nuclear forces, we considered symmetric displacements of the atoms in a water molecule along the OH bonds. In a small 6-31G basis set, exact (FCI) diagonalisation of the hamiltonian matrix is still feasible and Figure 2 illustrates results from FCI reference and -FCIQMC calculations. The nuclear forces as shown in Figure 2 have been obtained from the Cartesian force vectors as the absolute force acting on either a hydrogen or the oxygen atom with the sign taken from the z-component of the force vector, which has been aligned with the symmetry principal axis. Although there is no computational advantage over direct diagonalisation methods for basis sets as small as the 6-31G basis, the replica algorithm implemented in -FCIQMC can be applied to much larger molecules and basis sets, providing essentially numerically exact nuclear forces. In order to demonstrate the scope of the -FCIQMC replica technology, we have also computed the all-electron forces within a cc-pVTZ basis set, evidently an infeasible task for current deterministic FCI algorithms, where the many-body basis now spans determinants. Figure 2 (dashed lines) illustrates the notably larger forces at intermediate stretching of the OH bonds if accurate cc-pVTZ basis set are combined with this level of theory in the calculations. This would have implications for dynamics calculations, as well as providing the basis for highly accurate geometry optimisations for systems with electronic ground states of strong multi-reference character.

## 4The dipole moment of CO

The interaction of an electronic system of charge with an external electric field, , in an external potential, , may be expressed as an expansion in terms of multipoles,

with the rank-1 dipole moment, the rank-2 quadrupole moment, and so on. It is the dipole moment itself with which we are presently concerned, and which may be calculated according to:

where, in the last line, the substitution (for electrons) has been made. Applying the Slater–Condon rules,[47] this expression can be recast in terms of the one-body reduced density matrix and one-electron molecular-orbital integrals for an arbitrary Cartesian component, , as,

to which the contribution from the (fixed) nuclei with charges and positions has been added. Thus, given the molecular-orbital integrals, , , and , which are readily available,[49] the one-body reduced density matrix obtained from -FCIQMC provides direct access to the dipole moment, and, more generally, to multipole moments of arbitrary rank.

As an interesting application of this approach, we consider the well-known problem of the dipole moment of CO at its equilibrium bond length, .[50] This system, with its subtle combination of and effects, is difficult to predict intuitively *a priori*, and Hartree–Fock theory notably suggests the polarity to be CO, while it is experimentally known to be CO.

We use the large aug-cc-pVZ-DK basis sets for this study and adopt the second-order Douglas–Kroll–Hess hamiltonian.[51] Although relativistic effects are small for these comparatively light atoms, the calculation of the dipole moment tends to be strongly basis-set dependent, and the use of a large set becomes correspondingly desirable. To that same end, it is desirable to be able to extrapolate finite-basis dipole moments to the complete-basis-set limit, as such extrapolations have previously been useful in -FCIQMC studies.[12] It has been shown that the asymptotic convergence of the correlation part of the dipole moment with the cardinality of the basis set, , is suitably described by the form

in much the same way as the correlation energy itself.[56] The complete-basis-set limit correlation contribution to the dipole moment, , may thus be derived from two consecutive finite-basis results, of cardinality and , according to

to which the Hartree–Fock contribution in a suitably large basis (aug-cc-pV5Z-DK is used here, for which ) may then be added to obtain the total dipole moment.

D | T | Q | (DT) | (TQ) | |

HF | -0.10135 | -0.10435 | -0.10369 | - | - |

MRCI | 0.07175 | 0.07203 | 0.07066 | 0.07419 | 0.06929 |

CCSD | 0.06829 | 0.05594 | 0.05087 | 0.05278 | 0.04681 |

-FCIQMC | 0.05893(3) | 0.05200(4) | 0.0474(4) | 0.05112 | 0.0437 |

Table 1 presents the results of this approach, with the extrapolations performed from the double- and triple- and the triple- and quadruple- basis sets, alongside the analogous results from coupled-cluster theory and multi-reference CI. The rapid convergence of with number of walkers in the -FCIQMC dynamic is also demonstrated in Figure 3.

By comparison with the experimental dipole moment, variously given as and ,[61] it is apparent that -FCIQMC performs rather better than MRCI, and is comparable to CCSD. However, it can be seen that CCSD actually overestimates the dipole moment compared to -FCIQMC, which can be taken as close to exact in each of the finite basis sets, and this feature of CCSD allows for favourable cancellation of errors with the basis-set incompleteness, yielding the fortuitously accurate extrapolated result. The remaining disparity between these results and experiment should not be ascribed to an inadequacy of the -FCIQMC density matrices, but is rather largely attributable to basis-set incompleteness error. Indeed, the larger (TQ) extrapolations are rather more satisfactory than the corresponding (DT) results, highlighting the sensitivity of such approaches to the adequacy of the choice of basis. This effect has been previously observed in the context of ionisation potentials,[12] but is magnified in this instance by the stronger basis-set dependence of dipole moments than correlation energies. It is also worth noting that a small vibrational contribution to the dipole moment is expected,[65] but the results of this study support the view expounded in that work by Luis and coworkers, that an accurate treatment of electron correlation in a sufficiently large basis set is adequate for close agreement with experiment.

The quality of the -FCIQMC density matrices may be illustrated by considering the CO problem in a small cc-pVDZ basis, for which deterministic FCI results can be obtained. In this case, whose results are summarised in Table 2, -FCIQMC reproduces the FCI dipole moment to within , whilst the CCSD and MRCI results are in error by and , respectively. Also of note is that, whilst the quoted -FCIQMC result was obtained using walkers, it can be obtained just as well, and with apparently negligible initiator error, with only .

Abs. relative error (%) | ||
---|---|---|

HF | -0.0915 | 201.10 |

MRCI | 0.0973 | 7.61 |

CCSD | 0.0996 | 9.94 |

CCSDT | 0.0931 | 2.87 |

CCSDTQ | 0.0906 | 0.11 |

CCSDTQP | 0.0905 | - |

-FCIQMC | 0.09045(3) | 0.06 |

FCI | 0.0905 | - |

These results, therefore, bear out the supposed high quality of the sampled density matrices, and in demonstrating the compatibility of -FCIQMC with the Hellmann–Feynman theorem, suggest that future studies of energy derivatives and their associated properties may well prove fruitful.

## 5Atomic dipole polarisabilities

The previous section began by noting the dependence upon the permanent dipole moment of a system’s interaction with an applied electric field as given by . Of course, the application of such a field will, in reality, affect the distribution of charge, and hence the dipole moment itself. Expanding the dipole moment as a function of the field, therefore, we may write a given component, , as

where and represent elements of the polarisability and first hyperpolarisability tensors, respectively.[68] is the zero-field permanent dipole, which is always zero for atomic species. Whilst the effect of the induced dipole moment is generally less significant for polar systems, it is the leading-order term in the expansion of the dipole moment for atoms which has no static dipole. It is thus crucial in accounting for the dipole-dipole dispersion interactions which often bind such species, and indeed will be the first-order response not only to static, but also to dynamic fields.[69] The calculation of thus provides an interesting study in and of itself, as well as a probing test of the calculation of reduced density matrices with -FCIQMC. We here consider the noble-gas atoms, Ne, Ar, and Kr, as archetypal examples of the problem.

It is apparent from Equation 6 that the polarisability may be thought of as the derivative of the dipole with respect to the field,

evaluated at . As for many response properties, this may be calculated by solution of the coupled perturbed Hartree–Fock equations,[70] but for our purposes it is convenient to suppose that a particular component, , might be effectively calculated by a finite-difference approach,

in which is a small field applied in the direction, and is the th component of the dipole moment induced by so doing. The second equality holds for the spherically-symmetric atomic systems under consideration here since , and the errors resulting formula are second-order since it is now equivalent to a central-difference approximation.

Straightforward and appealing though this implementation is, it is useful before proceeding to have some notion of its performance relative to analytic gradient methods. In particular, analytic gradients are readily and rapidly available from MP2 theory,[49] and this approach thus provides a useful framework in which to assess the suitability of the finite-field method.

System | Analytic | Finite-field | Abs. relative error |

Ne | 2.438384 | 2.437962 | 1.73 |

Ar | 10.841398 | 10.842952 | 1.43 |

Kr | 16.674939 | 16.680012 | 3.04 |

Analytic | Finite-field | Abs. relative error | |

Ne | 2.620174 | 2.619806 | 1.40 |

Ar | 11.128300 | 11.131348 | 2.74 |

Kr | 16.792296 | 16.799500 | 4.29 |

The results of this comparison, with the finite-field polarisabilities performed in a field of strength 0.005 , are summarised in Table 3. The mean absolute percentage error inherent in the approach is found to be of the order of , demonstrating both its suitability for the purpose, and also that the field chosen is sufficiently small to establish the pseudo-linear dependence of the induced dipole upon the field. That this dependence is established without having to use a very small field is encouraging, since in the stochastic formulation provided by -FCIQMC, the stochastic error in the induced dipole must be divided by the field strength to obtain the equivalent error bounds in the polarisability. This behaviour is illustrated for Ne in the aug-cc-VTZ-DK basis in Figure 4, which highlights the balance which must be achieved between minimising second-order effects and maintaining a suitable level of stochastic error.

Secure in the knowledge of the suitability of the finite-field approach, we may proceed with an assessment of the performance of -FCIQMC compared with other methods. Specifically, as in the previous section, we calculate the polarisabilities using CCSD and MRCI for comparison,[49] though in this case the extrapolations are performed from results at the triple- and quadruple- basis sizes, reflecting the increased sensitivity to basis set incompleteness error which this quantity entails.

System | aug-cc-pVTZ-DK | aug-cc-pVQZ-DK | CBS | Experiment |

Ne | 2.42(1) | 2.596(2) | 2.65 | 2.57 |

Ar | 10.855(5) | 11.092(3) | 11.08 | 11.23 |

Kr | 16.81(4) | 16.86(6) | 16.82 | 16.73 |

As might have been expected, the error incurred by extrapolating is somewhat reduced upon including the larger quadruple- treatment, and the -FCIQMC results given in Table 4 bear correspondingly close agreement with experiment. The remaining errors — in the region of to — are nonetheless still likely to be artefacts of the basis sets, as the application of a field accentuates the importance of describing the intricacies of the more diffuse regions of electron density. Thus, although “augmented” basis sets are employed, there is likely still something to be gained from a more complete description of this behaviour. This suggestion is, once again, further strengthened by the fact that -FCIQMC is able to recover FCI-quality results for basis sets in which direct comparison is possible, reproducing the polarisability of Ne in a small cc-pVDZ basis to within , for instance.

The same results, computed using Hartree–Fock theory, CCSD,[49] and MRCI,[58] are listed in Table 5. The mean (absolute) error for the MRCI calculations is , whilst that for -FCIQMC, and coupled-cluster theory, is around . The comparability is unsurprising, given the ascription of much of the error to finite-basis effects. However, it is now necessary to investigate the impact of stronger correlation on this quantity in more challenging systems, where we expect more significant advantages to come from -FCIQMC.

System | HF | CCSD | MRCI | HF | CCSD | MRCI | CCSD | MRCI |

Ne | 2.20 | 2.42 | 2.41 | 2.33 | 2.59 | 2.58 | 2.64 | 2.64 |

Ar | 10.45 | 10.81 | 10.36 | 10.72 | 11.03 | 11.53 | 11.00 | 12.20 |

Kr | 16.21 | 16.78 | 17.15 | 16.36 | 16.81 | 17.22 | 16.74 | 17.18 |

## 6Conclusions

The results presented in this work serve to confirm the high quality of the stochastically-obtained reduced density matrices available *via* replica sampling in -FCIQMC, capable as they are of reproducing FCI-quality results for nuclear forces, dipole moments, and polarisabilities, and in some cases close agreement with experimental values. In so doing, they cement the place of the replica technique as an important extension to the theory, and widen its scope considerably.

In addition to the most obvious extension of an ability to compute a larger range of properties for a wider variety of systems, there remain a number of theoretical and technical challenges to be addressed in future studies. Perhaps the most pressing task is to extend this work to encompass results from open-shell systems, in which correlation effects are likely to be more important. Moreover, if comparisons to experimental results are to be further sought and achieved for dipole moment properties, there is some motivation to explore larger basis sets with multiple levels of augmentation,[73] which may be of particular use in better describing the more diffuse electron densities of finite-field calculations, and more generally in describing larger and heavier atoms of interest.

## 7Acknowledgements

The authors would like to thank Gerald Knizia for helpful discussions. The authors gratefully acknowledge Trinity College, Cambridge, and the Royal Society *via* a University Research Fellowship for funding. This work was also supported through a research fellowship of the Deutsche Forschungsgemeinschaft (D. O.) and by EPSRC under Grant No. EP/J003867/1. The calculations made use of the facilities of the Max Planck Society’s Rechenzentrum Garching.

### References

- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- (), bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- , ed. (, , )bibitemNoStop
- (, , )bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- in (, , )bibitemNoStop
- (), bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- (, , )bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- (, , )bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemStop
- bibitemNoStop
- , ed. (, , )bibitemNoStop
- , ed. (, , )bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop