Direct evaluation of the force constant matrix in quantum Monte Carlo

TL;DR

This paper introduces a formalism for directly evaluating force constant matrices in Quantum Monte Carlo, enabling highly accurate molecular relaxation and vibrational frequency calculations that outperform many existing methods.

Contribution

It presents a novel approach to compute force constants directly within Quantum Monte Carlo, improving the accuracy of molecular geometries and vibrational frequencies.

Findings

Bond lengths within 0.007 Å of experimental values.

Vibrational frequencies within 0.1% for H2 and HCl.

Vibrational frequencies within 1.1% for CO2 and CH4.

Abstract

We develop a formalism to directly evaluate the matrix of force constants within a Quantum Monte Carlo calculation. We utilize the matrix of force constants to accurately relax the positions of atoms in molecules and determine their vibrational modes, using a combination of Variational and Diffusion Monte Carlo. The computed bond lengths differ by less than 0.007{\AA} from the experimental results for all four tested molecules. For hydrogen and hydrogen chloride, we obtain fundamental vibrational frequencies within 0.1% of experimental results and ~10 times more accurate than leading computational methods. For carbon dioxide and methane, the vibrational frequency obtained is on average within 1.1% of the experimental result, which is at least 3 times closer than results using Restricted Hartree-Fock and Density Functional Theory with a Perdew-Burke-Ernzerhof (PBE) functional and comparable or better than Density Functional Theory with a semi-empirical functional.

Figures

Tables

Table 1. Table 1: Vibrational frequencies, evaluated at the computed equilibrium bond length, for the hydrogen molecule in units of cm − 1 superscript cm 1 \text{cm}^{-1} , where ω RHF subscript 𝜔 RHF \omega_{\text{RHF}} denotes the vibrational frequency obtained from a RHF calculation, ω PBE subscript 𝜔 PBE \omega_{\text{PBE}} and ω B3LYP subscript 𝜔 B3LYP \omega_{\text{B3LYP}} from DFT calculations with a PBE /B3LYP functional, ω EC subscript 𝜔 EC \omega_{\text{EC}} from the curvature of the DMC energy, ω FG subscript 𝜔 FG \omega_{\text{FG}} from the gradient of the DMC force, ω DMC subscript 𝜔 DMC \omega_{\text{DMC}} from the DMC matrix of force constants, and ω exp subscript 𝜔 exp \omega_{\text{exp}} from experiment. All values are presented at zero temperature and post anharmonic corrections. These quantities, as well as their associated errors, are discussed in Appendix  B .

mode	${\omega }_{\text{RHF}}$	${\omega }_{\text{PBE}}$	${\omega }_{\text{B3LYP}}$	${\omega }_{\text{EC}}$	${\omega }_{\text{FG}}$	${\omega }_{\text{DMC}}$	${\omega }_{\text{exp}}$Dickenson et al. (2013)
stretch	$4379$	$4116$	$4384$	$4170\pm 10$	$4180\pm 8$	$4166\pm 4$	$4161.1663\pm 0.0002$

Table 2. Table 3: Computed equilibrium bond lengths for the hydrogen, hydrogen chloride, carbon dioxide, and methane molecules, in units of Å, where x 0 RHF superscript subscript 𝑥 0 RHF x_{0}^{\text{RHF}} denotes the equilibrium bond length obtained from a RHF calculation, x 0 PBE superscript subscript 𝑥 0 PBE x_{0}^{\text{PBE}} and x 0 B3LYP superscript subscript 𝑥 0 B3LYP x_{0}^{\text{B3LYP}} from DFT calculations with a PBE/B3LYP functional, x 0 DMC superscript subscript 𝑥 0 DMC x_{0}^{\text{DMC}} from a DMC calculation, and x 0 exp superscript subscript 𝑥 0 exp x_{0}^{\text{exp}} from experiment. The details of the atomic relaxation calculation in DMC are discussed in Appendix  A .

molecule	$x_0^{\text{RHF}}$	$x_0^{\text{PBE}}$	$x_0^{\text{B3LYP}}$	$x_0^{\text{DMC}}$	$x_0^{\text{exp}}$ Sutton (1965)
${\text{H}}_2$	$0.736$	$0.753$	$0.745$	$0.7420\pm 0.0007$	$0.74130$
HCl	$1.260$	$1.286$	$1.278$	$1.280\pm 0.006$	$1.275$
${\text{CO}}_2$	$1.145$	$1.182$	$1.171$	$1.167\pm 0.003$	$1.1598$
${\text{CH}}_4$	$1.089$	$1.104$	$1.098$	$1.097\pm 0.002$	$1.093$

Equations

$$\hat{H}=$$

$$E=\frac{\int {\Psi }^{\ast }\hat{H}\Psi \mathrm{d}\mathbf{r}}{\int |\Psi {|}^2\mathrm{d}\mathbf{r}},$$

$$\frac{\text{d}E}{\text{d}\mathbf{R}_{I}}=\frac{\int\Psi^{*}\frac{\text{d}\hat{H}}{\text{d}\mathbf{R}_{I}}\Psi\mathrm{d}\mathbf{r}}{\int|\Psi|^{2}\mathrm{d}\mathbf{r}}+\Bigg{[}\frac{\int\frac{\text{d}\Psi^{*}}{\text{d}\mathbf{R}_{I}}(\hat{H}-E)\Psi\mathrm{d}\mathbf{r}}{\int|\Psi|^{2}\mathrm{d}\mathbf{r}}+\text{c.c.}\Bigg{]}\,,$$

$$\frac{{\text{d}}^{2}E}{\text{d}{\mathbf{R}}_I\text{d}{\mathbf{R}}_J}=\frac{\int {\Psi }^{\ast }\frac{{\text{d}}^2\hat{H}}{\text{d}{\mathbf{R}}_I\text{d}{\mathbf{R}}_J}\Psi \mathrm{d}\mathbf{r}}{\int |\Psi {|}^2\mathrm{d}\mathbf{r}}$$

$$\displaystyle\;\;\;\;\;\;\;\;\;+\frac{\int\Psi^{*}\Big{[}\frac{\text{d}\Psi}{\text{d}\mathbf{R}_{I}}\left(\Psi^{-1}\frac{\text{d}\hat{H}}{\text{d}\mathbf{R}_{J}}\Psi-\frac{\text{d}E}{\text{d}\mathbf{R}_{J}}\right)+(I\leftrightarrow J)\Big{]}\mathrm{d}\mathbf{r}}{2\int|\Psi|^{2}\mathrm{d}\mathbf{r}}$$

$$\displaystyle\;\;\;\;\;\;\;\;\;+\frac{\int\Big{[}\frac{\text{d}}{\text{d}\mathbf{R}_{J}}\left[\frac{\text{d}\Psi}{\text{d}\mathbf{R}_{I}}\left(\hat{H}-E\right)\Psi\right]+(I\leftrightarrow J)\Big{]}\mathrm{d}\mathbf{r}}{2\int|\Psi|^{2}\mathrm{d}\mathbf{r}}+\text{c.c.}.$$

$${\Psi }_{\text{T}}={\mathrm{e}}^J{D}_{\uparrow }{D}_{\downarrow },$$

$$E_{\text{VMC}}=\frac{\int |{\Psi }_{\text{T}}(\mathbf{r}){|}^2{E}_{\text{L}}(\mathbf{r})\mathrm{d}\mathbf{r}}{\int |{\Psi }_{\text{T}}(\mathbf{r}){|}^2\mathrm{d}\mathbf{r}},$$

$$E_{\text{DMC}}=\frac{\int \Phi (\mathbf{r}){\Psi }_{\text{T}}(\mathbf{r})E_{\text{L}}(\mathbf{r})\mathrm{d}\mathbf{r}}{\int \Phi (\mathbf{r}){\Psi }_{\text{T}}(\mathbf{r})\mathrm{d}\mathbf{r}}.$$

$$O=2O_{\text{DMC}}-O_{\text{VMC}}+\mathcal{O}\left[(\Phi -{\Psi }_{\text{T}}{)}^2\right].$$

$$\Delta \mathbf{R}=-2{\mathbf{M}}^{-1}\cdot {\nabla }_{\mathbf{R}}E-\frac{{\sum }_I{m}_I{\mathbf{R}}_I}{{\sum }_I{m}_I}-\mathbf{R}\times \mathbf{\theta },$$

$$\Delta \mathbf{R}=\Delta {\mathbf{R}}_{\text{e}}+\Delta {\mathbf{R}}_{\text{t}}+\Delta {\mathbf{R}}_{\text{r}},$$

$$E=E_0+\sum _{I=1}^{N_{\text{n}}}\frac{\mathrm{d}E}{\mathrm{d}{\mathbf{R}}_I}\Delta {\mathbf{R}}_I+\frac{1}{2}\sum _{IJ}\frac{{\mathrm{d}}^{2}E}{\mathrm{d}{\mathbf{R}}_I\mathrm{d}{\mathbf{R}}_J}\Delta {\mathbf{R}}_I\Delta {\mathbf{R}}_J,$$

$$\left[\frac{1}{2}\Delta {\mathbf{R}}^{⊺}\mathbf{M}+{\nabla }_{\mathbf{R}}E\right]\cdot \Delta \mathbf{R}=0,$$

$$\Delta {\mathbf{R}}_{\text{e}}=-2{\mathbf{M}}^{-1}{\nabla }_{\mathbf{R}}E,$$

$$\Delta {\mathbf{R}}_{\text{t}}=-\frac{{\sum }_I{m}_I{\mathbf{R}}_I}{{\sum }_I{m}_I}.$$

$$\sum _I{m}_I{\mathbf{R}}_I\times ({\mathbf{b}}_I-{\mathbf{R}}_I\times \mathbf{\theta })=0,$$

$$\sum _I{m}_I{\mathbf{R}}_I\times {\mathbf{b}}_I=\sum _I\left[m_I{\mathbf{R}}_I({\mathbf{R}}_I\cdot \mathbf{\theta })-m_I\mathbf{\theta }({\mathbf{R}}_I\cdot {\mathbf{R}}_I)\right],$$

$$\underset{\mathbf{a}}{\underbrace{\sum _I{m}_I\left(\begin{array}{c}{R}_{\mathrm{y}I}{b}_{\mathrm{z}I}-R_{\mathrm{z}I}{b}_{\mathrm{y}I}\\ R_{\mathrm{z}I}{b}_{\mathrm{x}I}-R_{\mathrm{x}I}{b}_{\mathrm{z}I}\\ R_{\mathrm{x}I}{b}_{\mathrm{y}I}-R_{\mathrm{y}I}{b}_{\mathrm{x}I}\end{array}\right)}}=\underset{\mathbf{B}}{\underbrace{\sum _I{m}_I\left(\begin{array}{ccc}-R_{\mathrm{y}I}^2-R_{\mathrm{z}I}^2& R_{\mathrm{y}I}{R}_{\mathrm{x}I}& R_{\mathrm{z}I}{R}_{\mathrm{x}I}\\ R_{\mathrm{x}I}{R}_{\mathrm{y}I}& -R_{\mathrm{x}I}^2-R_{\mathrm{z}I}^2& R_{\mathrm{z}I}{R}_{\mathrm{y}I}\\ R_{\mathrm{x}I}{R}_{\mathrm{z}I}& R_{\mathrm{y}I}{R}_{\mathrm{z}I}& -R_{\mathrm{x}I}^2-R_{\mathrm{y}I}^2\end{array}\right)}}\left(\begin{array}{c}{\theta }_{\mathrm{x}}\\ {\theta }_{\mathrm{y}}\\ {\theta }_{\mathrm{z}}\end{array}\right).$$

$$\Delta {\mathbf{R}}_{\text{r},I}=-{\mathbf{R}}_I\times \mathbf{\theta },$$

$$\hat{H}=\frac{{\hat{p}}^2}{2\mu }+\hat{V}$$

$$\hat{V}=V(x)=D[1-e^{-\alpha x}{]}^2,$$

$$E_n=\hslash {\omega }_0\left[\left(n+\frac{1}{2}\right)-x_{\text{e}}{\left(n+\frac{1}{2}\right)}^2\right],$$

$$\hat{H}=-\frac{{\hslash }^2}{2\mu }\frac{{\partial }^2}{\partial x^2}+\frac{{\hslash }^2{\alpha }^2}{2\mu }{\left(N+\frac{1}{2}\right)}^2(e^{-2x}-2e^{-x})$$

$$⟨0⟩\hat{x}0=\frac{\ln (2N+1)-\psi (2N)}{\alpha },$$

$$⟨0⟩\hat{x}0=\frac{3}{2}\sqrt{\frac{\hslash x_{\text{e}}}{2\mu {\omega }_0}}$$

$$V(x)=\frac{1}{2}\mu {\omega }_0^2{x}^2+\sqrt{\frac{{\mu }^3{x}_{\text{e}}{\omega }_0^5}{2\hslash }}{x}^3+\dots$$

$$x_{\text{e}}=\frac{\hslash }{18\sqrt{\mu }}{\left({\frac{{\text{d}}^{3}V}{\text{d}{x}^3}|}_{x_0}\right)}^2{\left({\frac{{\text{d}}^{2}V}{\text{d}{x}^2}|}_{x_0}\right)}^{-5/2}.$$

$$\mathbf{M}=\left(\begin{array}{cc}\frac{{\mathrm{d}}^{2}E}{\mathrm{d}{R}_1^2}& \frac{{\mathrm{d}}^{2}E}{\mathrm{d}{R}_1\mathrm{d}{R}_2}\\ \frac{{\mathrm{d}}^{2}E}{\mathrm{d}{R}_2\mathrm{d}{R}_1}& \frac{{\mathrm{d}}^{2}E}{\mathrm{d}{R}_2^2}\end{array}\right).$$

$$\begin{gathered} {\omega }^2=\frac{1}{2}\left(\frac{{\mathrm{d}}^{2}E}{\mathrm{d}{R}_1^2}+\frac{{\mathrm{d}}^{2}E}{\mathrm{d}{R}_2^2}\right) \\ \pm \sqrt{\frac{1}{4}\left(\frac{{\mathrm{d}}^{2}E}{\mathrm{d}{R}_1^2}-\frac{{\mathrm{d}}^{2}E}{\mathrm{d}{R}_2^2}\right)+{\left(\frac{{\mathrm{d}}^{2}E}{\mathrm{d}{R}_1\mathrm{d}{R}_2}\right)}^2}, \end{gathered}$$

$${\mathbf{M}}_{\text{KM}}=\frac{1}{2}\left(\begin{array}{cc}\frac{{\mathrm{d}}^{2}E}{\mathrm{d}(R_1+R_2)}& 0\\ 0& \frac{{\mathrm{d}}^{2}E}{\mathrm{d}(R_1-R_2)}\end{array}\right),$$

Full text

I Introduction

Quantum Monte Carlo (QMC) is a leading class of approaches used to establish and study the electronic ground state of molecules and solids. Specifically, Diffusion Monte Carlo (DMC) is widely used to project out the exact electronic ground state wave function of a system, subject only to the fixed node approximation, fully accounting for correlation effects such as van der Waals interactions Wheatley et al. (2004); Wheatley and Lillestolen (2007). Although DMC is an ideal tool for studying the electronic wave function of the system Kolorenč and Mitas (2011), the determination of the wave function of the atoms – their expected positions and energy landscape – remains a challenge for the method. Several approaches have been put forward to calculate the force acting on the atoms with DMC Lee et al. (2007); Badinski et al. (2010, 2008); Moroni et al. (2014); Filippi and Umrigar (2002); Assaraf and Caffarel (2000); Attaccalite and Sorella (2008); Filippi et al. (2016), but a more comprehensive characterization of the atomic wave function requires the second derivative of the energy – the matrix of force constants – to both efficiently relax atomic positions and calculate vibrational modes.

We propose a method to directly calculate the matrix of force constants, $\text{d}^{2}E/\text{d}\mathbf{R}_{I}\text{d}\mathbf{R}_{J}$, where $\mathbf{R}_{I}$ is the position of the $I$th, and $\mathbf{R}_{J}$ the $J$th, atom in the system. The energy, $E=\braket{\hat{H}}$, is calculated in the Born-Oppenheimer approximation of Hamiltonian $\hat{H}$; that is with the electrons always in their ground state for the respective atomic configuration. The calculation is implemented in QMC through a new quantum mechanical expectation value, $\text{d}^{2}\langle\hat{H}\rangle/\text{d}\mathbf{R}_{I}\text{d}\mathbf{R}_{J}$, meaning that it can be evaluated with one configuration of the atoms to recover the entire matrix of force constants. The matrix of force constants allows us to efficiently relax atomic positions and determine the vibrational modes.

We start by introducing the formalism and the QMC methods in Sec. II. We subsequently outline the applications and implementation of the matrix of force constants in Sec. III, followed by a series of case studies in Sec. IV. We begin with atomic and diatomic hydrogen, and then move on to hydrogen chloride, carbon dioxide, and methane. For each molecule we derive the matrix of force constants, relax the positions of the atoms, and determine the vibrational modes. We critically evaluate the results with respect to existing computational methods: Restricted Hartree-Fock (RHF) Hartree (1928); Slater (1928); Fock (1930) and Density Functional Theory (DFT) Hohenberg and Kohn (1964); Kohn and Sham (1965). Finally, in Sec. V we summarize the results and discuss future opportunities for the new formalism.

II Formalism

In this section, we present the matrix of force constants. We then outline the numerics by discussing how the electronic orbitals are generated and the details of the QMC algorithms.

II.1 Matrix of force constants

We consider many-body quantum systems comprised of $N_{\text{n}}$ nuclei and $N_{\text{e}}$ electrons. The three-dimensional position vectors are denoted as $\mathbf{R}_{I}$ for nuclei and $\mathbf{r}_{i}$ for electrons, with $I=1,\dots,N_{\text{n}}$ and $i=1,\dots,N_{\text{e}}$. These are used to construct the corresponding multi-dimensional vectors in configuration phase space: $\mathbf{R}\equiv(\mathbf{R}_{1},\dots,\mathbf{R}_{N_{\text{n}}})$ and $\mathbf{r}\equiv(\mathbf{r}_{1},\dots\mathbf{r}_{N_{\text{e}}})$.

We use the non-relativistic Hamiltonian Foulkes et al. (2001)

[TABLE]

which is comprised of the electron kinetic energy, as well as the electron-electron, electron-ion, and ion-ion interactions111We work with Hartree atomic units $\hbar=\text{m}_{\text{e}}=\text{e}=1$.. $V_{I}$ and $Z_{I}$ represent the electron-ion pseudopotential and full nuclear charge, of atom $I$, respectively.

We use a Hartree-Fock average effective Trail-Needs pseudopotential Drummond et al. (2016), which has been specifically optimized for DMC calculations, to screen the effects of the core electrons and nucleus on the valence electrons.

In an electron position basis, the expectation value of the energy Badinski et al. (2010) may be written as

[TABLE]

where the many-body wave function, $\Psi$, and Hamiltonian, $\hat{H}$, are both functions of nucleus configuration, $\mathbf{R}$, and electron configuration, $\mathbf{r}$.

The force acting on ion I is defined as the negative total derivative of the energy with respect to the nuclear coordinates. Taking the first derivative of the energy with respect to atom position Badinski et al. (2010) yields

[TABLE]

which is decomposed into Hellmann-Feynman and Pulay terms, respectively. When the wave functions are exact eigenstates of the Hamiltonian, such that $(\hat{H}-E)\Psi=0$, the Pulay term vanishes. However in practice, the wave functions are not exact in Variational Monte Carlo (VMC) or DMC, so the Pulay term needs to be included to obtain the total force.

In this paper, we derive the matrix of force constants from the second derivative of the energy that takes the form of

[TABLE]

This comprises one component of the matrix of force constants, so we must cycle over all atom pairs $\{I,J\}$ to determine the entire matrix. The second derivative of the Hamiltonian with respect to atom position does not commute with the Hamiltonian, hence we approximate the pure expectation value of the force constants in the DMC procedure, as discussed in Sec. II.3. The first two terms of the matrix of force constants stem from the Hellmann-Feynman force, whereas the third is due to the Pulay force.

To calculate the entire matrix of force constants using the Monte Carlo algorithm, we need to compute the ion-ion and electron-ion components at a cost of at most $O(N_{\textrm{n}}^{3}N_{\textrm{e}})+O(N_{\textrm{n}}^{2}N_{\textrm{e}}^{2})$. This leads to an overall dominant scaling of $O(N_{\textrm{n}}^{4})$ assuming $O(N_{\textrm{n}})\propto O(N_{\textrm{e}})$.

Having evaluated the matrix of force constants and implemented the formalism we can then use it to study atomic relaxation and vibrational modes.

II.2 Variational Monte Carlo

For the fermionic many-body trial wave function in the VMC method McMillan (1965), we take a Slater-Jastrow wave function of the form Foulkes et al. (2001):

[TABLE]

where $D_{\uparrow}(D_{\downarrow})$ denotes the Slater determinant of the molecular spin-up(down) orbitals. Here, the usual Hartree-Fock ansatz, $\Psi_{\text{HF}}=D_{\uparrow}D_{\downarrow}$, which encodes Pauli exclusion through the anti-symmetry of the Slater determinant, is multiplied by a Jastrow factor, $\text{e}^{J}$, which is an optimizable function used to impose further constraints on $\Psi_{\text{T}}$.

Initially, we compute the VMC energy, which is the expectation value of the Hamiltonian operator with respect to the trial wave function Badinski et al. (2010):

[TABLE]

where $E_{\text{L}}=\Psi_{\text{T}}^{-1}(\mathbf{r})\hat{H}\Psi_{\text{T}}(\mathbf{r})$ is the local energy, $\mathrm{d}\mathbf{r}$ is the infinitesimal hypervolume element in electron configuration phase space, and the integrals are performed using Monte Carlo222Specifically, the Metropolis algorithm is used to generate a set of configurations distributed according to the square modulus of a trial wave function over which the local energy is averaged. in the CASINO program Needs et al. (2010).

Single-particle orbitals for the different molecular structures were calculated using the CRYSTAL program Dovesi et al. (2014). The RHF and DFT calculations with two exchange-correlation functionals the Perdew-Burke-Ernzerhof (PBE) Perdew et al. (1996) containing no exact orbital exchange and the B3LYP Becke (1993); Stephens et al. (1994); Vosko et al. (1980) hybrid functional containing a fixed amount of exact exchange were performed with triple-$\zeta$-valence Gaussian basis sets, as well as polarization and diffuse basis functions Weigend and Ahlrichs (2005). The exact exchange-correlation functional is unknown and the choice of functional depends heavily on the system and the property of interest. PBE as a general functional was chosen for its greater predictive power across all simulations and properties Ernzerhof and Scuseria (1999), though it is less likely to achieve the accuracy of a semi-empirical functional such as B3LYP, which was chosen in addition for its good agreement with post-DFT methods within its range of applicability on molecules Mardirossian and Head-Gordon (2017); Paier et al. (2007). We use a Jastrow factor in its most general form comprising of an electron-electron term, an electron-nucleus term, and an electron-electron-nucleus term. The wave function parameters are optimized by using the variance minimization method Drummond and Needs (2005) first, followed by the energy minimization method Umrigar et al. (1988, 2007); Brown et al. (2007).

II.3 Diffusion Monte Carlo

DMC evolves a wave function, $\Phi$, according to the imaginary-time Schrödinger equation, in order to project out the lowest energy eigenstate, $\Phi_{0}$, with the same nodal surface as the trial wave function Foulkes et al. (2001).

The efficiency of the DMC algorithm is improved by importance sampling Ceperley and Alder (1980). By multiplying the wave function, $\Phi$, by a trial wave function, $\Psi_{\text{T}}$ from VMC, we may solve the Schrödinger equation for the mixed distribution $f(\mathbf{r},\tau)=\Phi(\mathbf{r},\tau)\Psi_{\text{T}}(\mathbf{r})$, where $\tau$ denotes imaginary time. We tested, and found negligible error, with time-steps of $\tau=0.01~{}\text{a.u.}$, and so this is used throughout Lee et al. (2011). The fixed-node approximation Anderson (1975, 1976) is introduced to overcome the fermion sign problem by constraining the nodal surface of $\Phi_{0}$ to match that of $\Psi_{\text{T}}$ 333The fixed-node approximation is the only uncontrolled approximation in a DMC simulation of all-electron systems..

The expectation value of the energy in the DMC method Badinski et al. (2010) is given by

[TABLE]

This is an unbiased estimator, up to the approximations made, since $E_{\text{DMC}}$ does not depend on the trial wave function used. However, the mixed expectation value of an operator that does not commute with the Hamiltonian is biased. In these cases, we approximate the pure expectation value of an operator $\hat{O}$ with the extrapolation formula Ceperley and Kalos (1979)

[TABLE]

Alternatively, the future-walking method may be used, for example, to obtain an exact pure estimator Barnett et al. (1991). Although the extrapolation formula improves the results, this procedure depends on an almost complete error cancellation and is strongly dependent on the quality of the wave function employed. We run the simulations for longer to systematically reduce the statistical error associated with variational techniques.

II.4 Contributions of the Hellmann-Feynman and Pulay terms

Both the Hellmann-Feynman and Pulay terms contribute to the force and matrix of force constants, and both Pulay terms are zero at the exact electronic ground state. However, the Pulay contribution to the matrix of force constants contains a second derivative of the wave function with respect to atom position, and so is more susceptible to steep gradients due to an incorrect trial wave function. Therefore, when using the electronic structure methods, it is useful to determine the relative contributions of the Hellmann-Feynman and Pulay terms so that we can gauge the importance of refining the electronic wave function.

The interatomic force in a hydrogen molecule and methane molecule is shown in Fig. 1. Different bond lengths within the vicinity of the equilibrium were chosen and forces were evaluated using the methods described in Secs. II.2 and II.3. The addition of the Pulay force to the Hellmann-Feynman force shifts the equilibrium bond length by 2% in both examples, therefore the Pulay force is crucial for finding the correct equilibrium geometry.

We now turn to consider the calculation of the matrix of force constants – the gradient of the force. We first note from Fig. 1 that the Pulay force is remarkably constant with respect to bond length across all molecules tested in this paper, regardless of the molecular geometry. This means that the gradient of the force is negligible, and therefore does not significantly contribute to the matrix of force constants. We find that, when directly evaluated, the value of the Pulay term in the matrix of force constants is smaller than its standard error, as well as the standard error of the contribution from the Hellmann-Feynman term. Furthermore, based on the analysis of variance at the $\alpha=0.05$ level, the gradient of the Pulay force does not significantly deviate from zero. This means that the change in vibrational frequency due to the Pulay force gradient is just 1% of that from the Hellmann-Feynman gradient for both hydrogen and methane. This conclusion is also backed up by independent studies: taking a numerical derivative of the results for $\text{H}_{2}$ and LiH reported by Casalegno et al. Casalegno et al. (2003), $\text{CO}_{2}$ reported by Lee et al. Lee et al. (2007), as well as adenine-thymine reported by Ruiz-Serrano et al. Ruiz-Serrano et al. (2012), confirms the small contribution of the Pulay term to the matrix of force constants.

The effect of the Pulay term in the force is significant for force analysis and, when suitably formulated, can reduce the statistical noise in the expectation value and improve the convergence of the optimization algorithm Ruiz-Serrano et al. (2012). However, as the Pulay force is almost constant with interatomic bond length, its contribution to the matrix of force constants is negligible. Therefore, we expect the Pulay contribution to the matrix of force constants to be insensitive to the quality of the trial wave function. Another corollary is that for our zero-variance scheme Assaraf and Caffarel (1999, 2003), the expected $-5/2$ power law tail associated with infinite variance that could arise from Pulay terms Trail (2008); Badinski et al. (2010) will only make a limited contribution to the tail of the total probability distribution and as shown in Fig. 4 we have not been hindered by this problem in our practical application.

III Applications of Force Constants

III.1 Atomic relaxation

The primary requirement for a versatile geometry model, is the ability to minimize the energy of an arbitrary configuration of atoms Wagner and Grossman (2010). For quantum mechanical simulations, this is often performed using VMC, due to the algorithm’s efficiency. In this paper, we relax the positions of the atoms first with VMC using the additional information provided by the matrix of force constants. The wave function from VMC is then optimized in DMC and we perform the same iteration steps using DMC to confirm convergence and further reduce the error Needs et al. (2010); Al-Saidi (2008); Barnett and Whaley (1993); Drummond et al. (2006).

Requiring that the energy of the system is constant up to quadratic order in atomic displacement, and explicitly correcting for global translation and rotation, as well as anharmonicity, we find that the atomic displacement, $\Delta\mathbf{R}$, is given by

[TABLE]

where $\mathbf{M}\equiv\text{d}^{2}E/\text{d}\mathbf{R}_{I}\text{d}\mathbf{R}_{J}$ is the matrix of force constants; $\nabla_{\mathbf{R}}E\equiv\text{d}E/\text{d}\mathbf{R}_{I}$ is the multi-particle gradient of the energy with respect to atom position; $\boldsymbol{\theta}$ is the three-dimensional global angular displacement vector for the configuration $\mathbf{R}$. On each step, we displace the atoms by $\Delta\mathbf{R}$ in order to compare with other methods in determining the minimum energy of the system. This yields the interatomic bond length at the minimum of the total potential. The details of the calculation are outlined in Appendix A.

Though $\langle\hat{\mathbf{R}}\rangle$ minimizes the total potential after the atomic relaxation procedure, if the potential is not symmetric then the expected separation of the atoms does not coincide with the minimum. We capture the lowest-order difference with the addition of an anharmonic correction term $\Delta\mathbf{R}_{\text{a}}$, outlined in Appendix A.4.

III.2 Vibrational modes

One main motivation for incorporating the matrix of force constants into the QMC procedure is the ability to calculate vibrational modes and frequencies directly. In this paper, we use a variety of methods to calculate the vibrational frequency for a cross comparison.

Up to a simple mapping, the eigenvectors of the matrix of force constants correspond to the vibrational modes of the system, and the eigenvalues correspond to the vibrational frequencies. We can, therefore, use a complete diagonalization of the matrix of force constants to estimate the eigenmodes and frequencies. To benchmark the results, we also calculate the frequencies from a numerical second derivative of the energy with respect to bond length – referred to as the energy curvature (EC) method – and from a numerical derivative of the force – force gradient (FG) method.

A discussion of all of these methods, including the analysis of statistical uncertainty and the anharmonic correction, is detailed in Appendix B.

IV Case Studies

In this section, we evaluate the effectiveness of the matrix of force constants formalism for a selection of molecules. We first confirm the theory with the simplest possible molecules, before testing the generalizability of the formalism on molecules containing more atoms.

IV.1 Hydrogen atom and molecule

We begin by analyzing the simplest physical system: the hydrogen atom. By performing a DMC calculation, we verify that the hydrogen atom obeys Newton’s laws since it has a net force of $(3.68\pm 5.17)\times 10^{-3}~{}E_{\text{h}}\text{\AA}^{-1}$ acting on it, which is zero within standard error. Furthermore, the hydrogen atom has a computed eigenfrequency of $0~{}\text{cm}^{-1}$. This system behaves as expected and confirms the translational invariance.

From this, it is natural to increment the complexity by adding another hydrogen atom to form an $\mathrm{H}_{2}$ molecule. This is the simplest physical example that allows us to verify the eigenfrequencies from our DMC method, which has no nodes and gives an exact wave function, by comparing them against both experimental results in the literature, and RHF/DFT predictions from the CRYSTAL program.

The energy, force, and diagonal elements of the matrix of force constants for the hydrogen atom is shown as a function of bond length in Fig. 2. We verify that the energy is at a minimum and the force is zero at the correct equilibrium bond length of $74.13~{}\text{pm}$ Sutton (1965), within standard error. Furthermore, in Fig. 2b, we show that in the vicinity of the equilibrium bond length for the hydrogen molecule, the energy curvature, the force gradient, and the direct computation, all agree within error bounds. The entire matrix of force constants is sparse. If we are only interested in the vibrational mode, it can be reduced to be a $2\times 2$ matrix, with the off-diagonal force constants to be minus the diagonal elements within error, as required by symmetry 444The diagonal elements of the matrix of force constants are the same for the hydrogen molecule, and hence no particular atom is specified when referring to $\mathrm{d}^{2}E/\mathrm{d}R^{2}$ in Figs. 2 and 4.. Note that here the diagonal elements of the matrix of force constants are not constant across the range of bond lengths shown, as can also be seen in the slight curvature of the force in Fig. 2a. This is due to the anharmonicity of the potential in a diatomic molecule. We may use the gradient of the force constant to calculate the anharmonic constant, and correct for the anharmonicity, as discussed in Appendix A.4.

For the hydrogen molecule, it is possible to extract the matrix of force constants efficiently from the force, or energy, because the computational cost of obtaining the numerical derivatives is low. However, for more complicated molecules, where structural optimization is influential, using the matrix of force constants would be beneficial, as it provides both the movement direction and amplitude towards the minimum energy configuration. In these cases, the equivalent information would take considerably longer to extrapolate from either energy or force, if possible.

Equipped with reliable results for the matrix of force constants directly from DMC at each bond length, we may now exploit this information to efficiently relax the bond length of the molecule. The force tells us the direction to move the atoms, and the matrix of force constants additionally tells us how far to move them, on each step (Eq. 1). Owing to the anharmonicity of the potential, we must relax to the equilibrium bond length of $74.13~{}\text{pm}$ over several steps. The predicted bond length on the next atomic relaxation step as a function of initial bond length is shown in Fig. 3. We see that the PBE curve intersects the equilibrium line at $0.753~{}\text{\AA}$ and the B3LYP curve intersects at $0.745~{}\text{\AA}$ ; whereas our DMC calculation intersects at $0.7420\pm 0.0007~{}\text{\AA}$, in close agreement with the experimental value of $0.74130~{}\text{\AA}$ Sutton (1965)555Note that the turning points of the PBE, B3LYP and DMC curves in Fig. 3 do not correspond to fixed points, but rather the fixed points are given by the intersections of the curves with the equilibrium line.. Furthermore, we note that PBE and B3LYP curves have a similar shape as a result of sharing the same optimization algorithm. However, both are steeper than the DMC curve in the vicinity of equilibrium and therefore converge more slowly, due to the fact that the DFT algorithm uses an inaccurate estimate for the force constant. The number of iteration steps reduced are particularly apparent when there are multiple atoms in a molecule, however the lower number of steps does not necessarily indicate a more efficient algorithm, as the complexity of each step needs to be taken into consideration. In some cases where complex molecules cannot be relaxed sufficiently for a long time using DFT, our approach may be more efficient in giving the ground state geometry. In general, due to DFT’s inaccurate estimate of the force constant, we observe at least a slight improvement for all molecules.

An additional important check, before we proceed, is an analysis of the probability distribution of the matrix of force constants generated by DMC. Fig. 4 shows the histogram of a force constant value for the DMC run at the computed equilibrium bond length. From this we can see that the force constant heavy tails decay with the same $\sim|M-M_{0}|^{-4}$ power law as the energy and Hellmann-Feynman force Trail (2008); Badinski et al. (2010). This is as expected, since the effective remaining term is the Hellmann-Feynman term due to the quasi-constant Pulay contribution in proximity to the ground state, as seen in Sec. II.4. Reassuringly, the expected value of the distribution is also the modal value.

Now that the configuration is relaxed, we may analyze the fundamental vibrational modes. For the hydrogen molecule, we obtain six eigenmodes, as expected for a diatomic molecule. Three modes correspond to global translation, two correspond to global rotation, and one corresponds to a symmetric stretch. The symmetric stretch mode has the largest eigenfrequency. We extract the frequency using a selection of methods, outlined in Appendix B, for a cross-comparison. In this case, we obtain a fundamental vibrational frequency of $\omega_{\text{DMC}}=4166\pm 4~{}\text{cm}^{-1}$, compared to the experimental value of $\omega_{\text{exp}}=4161.1663\pm 0.0002~{}\text{cm}^{-1}$, which is $4.83~{}\text{cm}^{-1}$ away. This is a firm statistical confirmation of the accuracy of DMC compared to RHF, PBE and B3LYP results, which have deviations of $218~{}\text{cm}^{-1}$, $45~{}\text{cm}^{-1}$ and $223~{}\text{cm}^{-1}$ from experiment, respectively. All of our computational values for the vibrational frequency from DMC – matrix of force constants, force gradient, and energy curvature – agree with each other within standard error and show close agreement to experiment. The DMC procedure yields no translational or rotational modes, just as for the hydrogen atom. A summary of the results is shown in Table 1. Note that the less computationally expensive calculation of the energy was run for four times longer, compared to the force gradient and matrix of force constant methods, in order to give the error bars of the energy curvature eigenfrequency to a comparable value.

IV.2 Hydrogen chloride

Now that we have verified that the matrix of force constants agrees with numerical estimates, and that by exploiting this information it is possible to relax the molecule more efficiently, and outperform RHF and DFT estimates for the fundamental vibrational frequency for the hydrogen case, we move onto a more complex molecule: hydrogen chloride. We increment the complexity of our case study in order to verify that our formalism can cope with an asymmetric system with unequal masses.

The hydrogen chloride molecule again relaxes quickly to equilibrium, with a computed bond length of $128.0\pm 0.6~{}\text{pm}$, which agrees with the experimental value of $127.5~{}\text{pm}$ within standard error. Both atoms have the appropriate displacement to ensure that the center of mass is stationary. We obtain six eigenmodes for the system, including one symmetric stretch mode with eigenfrequency $\omega_{\text{DMC}}=2995\pm 8~{}\text{cm}^{-1}$. This result agrees with the experimental value of $\omega_{\text{exp}}=2990.946\pm 0.003~{}\text{cm}^{-1}$ within standard error, whereas RHF, PBE, B3LYP methods are $107~{}\text{cm}^{-1}$, $112~{}\text{cm}^{-1}$, $50~{}\text{cm}^{-1}$ away, respectively. A summary of the results is shown in Table 2a.

IV.3 Carbon dioxide

In the previous two examples, we found that the matrix of force constants can correctly calculate the modes of a diatomic molecule. Building on this, we increment the complexity to carbon dioxide: a three-atom system with several non-trivial vibrational modes, some of which are in orthogonal directions.

In this case, the $\mathrm{O=C=O}$ configuration is relaxed to an equilibrium $\mathrm{C=O}$ bond length of $116.7\pm 0.3~{}\text{pm}$ along one axis, which is within three standard deviations of the experimental value of $115.98~{}\text{pm}$. For carbon dioxide, we obtain nine vibrational modes: three of which correspond to global translations, two to global rotations, and four to vibrational modes. Of the vibrational modes, we obtain one symmetric stretch mode, one asymmetric stretch mode, and two bending modes along orthogonal axes.

The modes examined in this section show a consistent improvement over the RHF and PBE calculations, with DMC eigenfrequency deviations from experiment of $-1.80$% (symmetric), $-1.58$% (antisymmetric), and $-0.75$% (bending). The recovery of the non-trivial antisymmetric mode is our first example to break the underlying symmetry of the molecule, and the bending mode shows that our formalism can extend to atoms moving in orthogonal directions. On average, our DMC result is $22~{}\text{cm}^{-1}$ away from the experimental value, which is an improvement over the RHF results (on average $122~{}\text{cm}^{-1}$ away) and PBE results (on average $45~{}\text{cm}^{-1}$ away). We note that in this particular case the B3LYP results are on average $13~{}\text{cm}^{-1}$ away from the experimental results, which is why it is a popular choice for non metal-containing molecules Peverati and Truhlar (2011).

It is worthwhile to mention that the experimental results come with a larger error for carbon dioxide when compared to smaller molecules, as shown in Table 2b. The symmetric stretch mode is Raman active and infrared inactive, whereas for the other modes, the opposite is true Shimanouchi (1972). The Raman measurement typically comes with a larger uncertainty than infrared spectroscopy in this case, complicating the comparison to DMC. Additionally, for these larger molecules, as the number of modes increases, the chances of eigenfrequency interference is increased. Here we notice that the symmetric stretch mode eigenfrequency is quasi-degenerate with twice the bending mode eigenfrequency in Table 2b, which could also potentially contribute to the increased uncertainty of the symmetric stretch mode McCluskey and Stoker (2006). Finally, we note that the precise eigenfrequencies for arbitrarily large molecules have not been studied as extensively. In contrast, the eigenfrequency calculation for hydrogen especially, as well as for other common diatomic molecules, is often used as an experimental benchmark Dickenson et al. (2013). Together these factors motivate the need for improving the accuracy and precision of electronic structure calculations, such as QMC.

IV.4 Methane

For the last example, we extend our formalism to a three-dimensional molecule containing five atoms: methane, to demonstrate that the formalism can be applied to diverse configurations of atoms.

We find that the configuration relaxes to a $\mathrm{C-H}$ bond length of $109.7\pm 0.2~{}\text{pm}$, within two standard deviations of the experimental value of $109.3~{}\text{pm}$, in fewer iterations than existing methods. The equilibrium bond lengths for all case studies are summarized in Table 3. In this case, we obtain fifteen eigenmodes of the system: three corresponding to global translation, three corresponding to global rotation, and nine corresponding to non-trivial vibrational modes. Of the vibrational modes, we select two modes to examine in detail: the symmetric stretch mode and a scissor mode, as summarized in Table 2c and illustrated in Fig. 5.

An analysis of these modes yields a DMC deviation from experiment of $-1.47$% for the symmetric stretch mode and an expected agreement for the scissor mode, which is generally comparable to the results for carbon dioxide i.e. still of the order of $1$% from the experimentally measured values. The successful recovery of these modes demonstrates that the formalism holds in three dimensions, and the excellent agreement for the scissor mode demonstrates that we are able to capture a non-trivial symmetry for this molecule. The symmetric stretch DMC eigenfrequency is $43~{}\text{cm}^{-1}$ away from experiment, whereas the RHF, PBE and B3LYP results are $184~{}\text{cm}^{-1}$, $117~{}\text{cm}^{-1}$ and $157~{}\text{cm}^{-1}$ away, respectively.

V Conclusion

In this paper, we develop and implement a formalism to evaluate the matrix of force constants in QMC. We calculate vibrational frequencies and improved estimates for the atomic displacements on each relaxation step, as well as correcting for anharmonicity. We report statistically significant improvements over RHF and DFT methods in the vast majority of cases, both in terms of the vibrational frequency and the efficiency of the atomic relaxation, for the hydrogen, hydrogen chloride, carbon dioxide, and methane molecules.

The ability to calculate the matrix of force constants within DMC, in particular, makes us well-positioned to calculate vibrational modes where high accuracy is a necessity and relax atomic positions in complex systems with many degrees of freedom where the extrapolation from energy or force is difficult, if not impossible, to optimize the geometry. The approach applies to both molecules and periodic configurations. This will be especially beneficial in systems with heavy atoms that are challenging to analyze accurately with DFT, systems with significant anharmonic corrections, and also those with strong van der Waals interactions, such as layered materials and surfaces.

Acknowledgements.

The authors thank Andrew Fowler, Neil Drummond, Christian Carbogno, Ryan Hunt, Pablo López Ríos, John Trail, Thomas Gasenzer, and Lars Schonenberg for useful discussions. The QMC calculations in this paper were performed using CASINO. Y.Y.F.L. acknowledges support from the Agency for Science, Technology and Research; B.A. acknowledges support from the Engineering and Physical Science Research Council under grant EP/M506485/1; and G.J.C. acknowledges support from the Royal Society. There is open access to this paper and data available at https://www.openaccess.cam.ac.uk.

Appendix A Atomic relaxation calculation

In this section, we describe in detail how the configuration coordinates are adjusted on each step during the atomic relaxation process.

Let us define the atomic displacement on each Monte Carlo step as

[TABLE]

where $\Delta\mathbf{R}_{\text{e}}$ is the energy-minimizing term, $\Delta\mathbf{R}_{\text{t}}$ is the correction for global translations, $\Delta\mathbf{R}_{\text{r}}$ is the correction for global rotations. We adjust the atomic displacements from $\mathbf{R}$ to $\mathbf{R}+\Delta\mathbf{R}$ on each step, so as to minimize the total energy of the system. Once the equilibrium is reached, the anharmonic correction $\Delta\mathbf{R}_{\text{a}}$ is applied.

A.1 Minimizing the energy

Consider a system of $N_{\text{n}}$ atoms in three dimensions. Taylor expanding the total energy of the system as a function of atomic displacements, up to quadratic order, yields

[TABLE]

where $E_{0}$ is a constant. Demanding that the sum of the first- and second-order terms in the energy are zero at the minimum, gives

[TABLE]

which, excluding the trivial solution, implies

[TABLE]

where $\mathbf{M}$ is the matrix of force constants, and $\nabla_{\mathbf{R}}E$ is the multi-atom energy gradient with respect to the configuration atomic displacement vector, $\mathbf{R}$. This is the bare estimate for the atomic displacement correction, up to second order in the energy.

A.2 Correction for global translations

In order to ensure that the origin of our configuration is fixed and that we have no global translational mode, we explicitly subtract the center of mass motion of the configuration.

Given $N_{\text{n}}$ atoms, each with mass $m_{I}$, this implies that the global translation correction term is

[TABLE]

This term is particularly important for non-symmetric molecules, such as hydrogen chloride in Sec. IV.2.

A.3 Correction for global rotations

Similarly, to ensure that the bond length corrections do not result in a rotation of the configuration, or atomic pair rotations, we explicitly subtract global rotational modes.

The law of moments states that the total moment about the center of mass of any atomic pair, as well as the total moment about the origin of the configuration, is zero, which gives $\sum_{I}m_{I}\mathbf{b}_{I}=\mathbf{0}$ and $\sum_{I}m_{I}\mathbf{R}_{I}=\mathbf{0}$, where $\mathbf{b}$ is the half-bond length between an atomic pair. Together, these relations imply

[TABLE]

where $\mathbf{b}_{\mathrm{r},I}\equiv\mathbf{b}_{I}-\mathbf{R}_{I}\times\boldsymbol{\theta}$ is the corrected half bond length to be found. Hence, an expression for the angular displacement of the molecule $\boldsymbol{\theta}_{I}\equiv(\theta_{\mathrm{x}I},\theta_{\mathrm{y}I},\theta_{\mathrm{z}I})$ is needed. Using the vector triple product identity, we find that Eq. 2 reduces to

[TABLE]

which after rearrangement becomes

[TABLE]

This implies that the atomic displacement correction for global rotations, is

[TABLE]

where $\boldsymbol{\theta}=\mathbf{B}^{-1}\mathbf{a}$.

A.4 Correction for anharmonicity

Up to this point in the analysis, we have assumed that the interaction between atomic pairs is harmonic. Although this is a valid approximation at short distances, at larger distances this approximation breaks down and so a correction term is necessary. One of the most well-studied models used to capture anharmonicity in the interaction between diatomic molecules is the Morse Hamiltonian, which we use as an approximation for our case studies. The Morse Hamiltonian is given by

[TABLE]

with a Morse potential

[TABLE]

where $D$ is the $x=x_{0}$ energy minimum depth relative to the dissociation limit at $x\to\infty$ and $\alpha$ determines the curvature of the potential Morse (1929).

The eigenvalues of the Morse Hamiltonian are

[TABLE]

where $\omega_{0}=\sqrt{2D\alpha^{2}/\mu}$ is the fundamental frequency, $x_{\text{e}}=\hbar\omega_{0}/4D$ is the anharmonic constant, and $n\in\mathbb{Z}^{+}$ is the principal quantum number.

Note that the minima of the harmonic and Morse potentials are the same. However, due to the dissociative limit of the Morse potential, the expectation value of position is shifted in the positive x direction in the Morse case. One of the main advantages of this model is that the majority of its properties can be expressed analytically.

By setting $D=\frac{\hbar^{2}\alpha^{2}}{2\mu}(N+1/2)^{2}$, the Morse Hamiltonian may be written as

[TABLE]

up to a constant term. The expectation value of position with respect to the ground state Morse wave function is then

[TABLE]

where $\psi$ is the digamma function de Lima and Hornos (2005). Expanding the expectation value of position gives

[TABLE]

up to leading order in $x_{\text{e}}$. This is the shift in the equilibrium bond length due to the anharmonicity of the Morse potential.

In order to evaluate this shift, an estimate for the anharmonic constant is needed. Expanding the Morse potential (Eq. 3) about the equilibrium displacement $x=x_{\text{0}}$ in powers of $x$, we find that

[TABLE]

up to a constant term. Comparing quadratic and cubic terms in $x$ with the general form of the Taylor expansion, and solving simultaneously, yields

[TABLE]

Conventionally, the third derivative of the energy is extracted from the curvature of the force, however now utilizing the new information available, we extract the anharmonic constant directly from the gradient of the force constant.

Appendix B Vibrational modes calculation

In this section, we describe in detail the methods used to determine the vibrational modes and frequencies of atomic configurations, as well as their associated statistical uncertainties.

B.1 Exisiting computational approaches

In order to calculate an estimate for the frequency using the RHF and DFT methods, we use the default scheme, PBE and B3LYP exchange-correlation functionals, respectively, within the CRYSTAL program Dovesi et al. (2014).

B.2 Matrix of force constants approach

The direct method to obtain the vibrational frequencies of a molecule is from an exact diagonalization of the matrix of force constants. Consider, for example, a diatomic molecule in one dimension, such as the hydrogen molecule discussed in Sec. IV.1. The matrix of force constants for this system may be written as

[TABLE]

By exactly diagonalizing the matrix, we obtain the eigenmodes, and eigenfrequencies of the system given by

[TABLE]

where the positive frequencies are physical. The errors are calculated using Monte Carlo, as discussed in Sec. B.5.

There are two possible disadvantages of this method for obtaining the vibrational frequencies of a configuration. First, since it is a complete diagonalization method, it uses all of the entries in the matrix of force constants. However, many of these entries are related by symmetries, and so these calculations are potentially redundant. Second, due to numerical inaccuracy, Eq. 5 may result in an overestimate of the frequencies if the diagonal terms in Eq. 4 are not equal.

Following from the previous example, by imposing the known modes of a diatomic molecule in one dimension, we may write the matrix of force constants as

[TABLE]

which now yields the eigenfrequencies

[TABLE]

Notice that $|\omega_{\text{KM}}|\leq|\omega|$ due to the absence of the diagonal terms in the square root of Eq. 5.

For a general system, we may input a set of known modes $\{\mathbf{x}\}$. These $3N_{\text{n}}$-dimensional row vectors act on the $3N_{\text{n}}\times 3N_{\text{n}}$ dynamical matrix, $\mathbf{D}$, to extract the corresponding eigenfrequency, such that

[TABLE]

with corresponding error

[TABLE]

where the hats denote normalization, $\boldsymbol{\Sigma}$ is the standard error matrix corresponding to $\mathbf{M}$, and the dynamical matrix, $\mathbf{D}$, is the matrix of force constants weighted by the atomic masses.

By imposing known modes on the system, we can reduce the potential for numerical error and speed up the matrix diagonalization. However, these advantages only hold if the correct eigenmodes are known a priori, and therefore we do not employ this scheme as standard for our DMC calculations.

B.3 Approaches based on derivatives of the force and energy

Further to the methods based on the matrix of force constants, we also consider traditional techniques, for comparison.

We obtain an estimate of the frequency ($\omega_{\text{FG}}$) from the gradient, $\kappa$, of the interatomic force against bond length graph. The error in the gradient of the slope is the asymptotic standard error from a linear regression, and this is propagated to the vibrational frequency in the usual way:

[TABLE]

Similarly, an additional estimate of the vibrational frequency ($\omega_{\text{EC}}$) is obtained by computing the second derivative of the energy at a series of displacements along the trajectory of an eigenmode. For this, we use a numerical central difference scheme. Since this result is based on a linear superposition of energy data points, the errors add in quadrature.

B.4 Correction for anharmonicity

All of the above methods for calculating the vibrational frequency rely on the harmonic potential approximation. However, there are certain cases where anharmonic vibration is dominant and a correction to these frequencies needs to be applied. As for atomic relaxation, we apply an approximate correction, due to a Morse potential, which for the fundamental vibrational frequency, is given as $\Delta\omega=-x_{\text{e}}/4$, where $x_{\text{e}}$ is the anharmonic constant.

B.5 Monte Carlo uncertainty

The matrix of force constants $\mathbf{M}$ comes with an associated standard error matrix, $\boldsymbol{\Sigma}$, from the reblocking method in CASINO Flyvbjerg and Petersen (1989). Calculating the errors in eigenvalues given the errors in the matrix elements is a non-trivial task and one which has been studied extensively in pure mathematics Alefeld and Spreuer (1986); Behnke (1988); Dongarra et al. (1983); Rump (1989); Yamamoto (1980, 1982). For the purposes of this paper, we calculate the eigenvalue errors using Monte Carlo.

For each Monte Carlo run we generate a dynamical matrix, whose matrix elements are normally distributed, with a mean equal to the original matrix elements and a standard deviation equal to the corresponding standard errors. We then perform many runs until the average eigenvalues converge to the true eigenvalues, and we use the standard errors of these Monte Carlo runs as the errors in the eigenvalues.