Atomic forces by quantum Monte Carlo: application to phonon dispersion calculation

We report the first successful application of the {\it ab initio} quantum Monte Carlo (QMC) framework to a phonon dispersion calculation. A full phonon dispersion of diamond is successfully calculated at the variational Monte Carlo (VMC) level, based on the frozen-phonon technique. The VMC-phonon dispersion is in good agreement with the experimental results, giving renormalized harmonic optical frequencies very close to the experimental values, by significantly improving upon density functional theory (DFT) in the generalized gradient approximation. Key to success for the QMC approach is the statistical error reduction in atomic force evaluation. We show that this can be achieved by using well conditioned atomic basis sets, by explicitly removing the basis-set redundancy, which reduces the statistical error of forces by up to two orders of magnitude. This leads to affordable and accurate QMC-phonons calculations, up to $10^{4}$ times more efficient than previous attempts, and paves the way to new applications, particularly in correlated materials, where phonons have been poorly reproduced so far.

The accurate description of phonons in a solid is one of the central research topics in the field of condensed matter physics and materials science for discussing phase stability (i.e., Gibbs-free energy), electron-phonon interaction, structural phase transitions of materials 1,2 . Ab initio phonon calculations based on Density Functional Theory (DFT) 3,4 have been successful for many compounds, but they often fail in strongly-correlated materials. For example, DFT calculations severely underestimate the highest frequency of the optical phonons of graphene at K point 5,6 , because the electronelectron correlation is not taken into account with sufficient accuracy. Another example is the elemental cerium, whose phonon dispersions measured by neutron scattering strongly mismatch with the calculated DFT-PBE ones 7 . Some effort has been made to include correlation in phonon calculations in the DFT+U framework 8 and also within the DMFT framework 9,10 . In both cases, this requires modeling correlations by an empirical parameter, though physically motivated (i.e., the Hubbard U). Indeed, a genuine ab initio framework applicable to strongly correlated materials without any empirical parameters remains, so far, a very important theoretical challenge.
The ab initio quantum Monte Carlo (QMC) framework, including variational quantum Monte Carlo (VMC) and the diffusion quantum Monte Carlo (DMC) schemes, is among the state-of-art numerical methods to obtain highly accurate many-body wave functions 11 , and cope with the electron correlation more rigorously than DFT. It has been successfully applied to challenging materials that DFT cannot tackle, such as cuprates 12 , iron-arsenides 13 , and graphene 12,14 . So far, unfortunately, almost all QMC applications are mainly based on energy and its first derivative (i.e., atomic force) calculations. Indeed, it is at present an open problem how to evaluate, with an affordable computational effort, second and higher-order derivatives, which are essential for computing various physical properties.
There are three routes to compute the second derivatives (i.e., ∂ 2 E ∂R α ∂R β ), which are needed for evaluating harmonic phonon properties, by ab initio calculations, i.e., ( i ) potential energy surface (PES) fitting, (ii) finite difference expression based on atomic force evaluations , and (iii) direct evaluation of second derivatives. All the above attempts have been successful for isolated molecular systems [15][16][17][18] . On the other hand, for solids, only the strategy ( i ) has been successful so far within the QMC framework. For example, Maezono et al. calculated Raman frequencies of diamond (phonons at q = Γ point) 19 . However, QMC-phonon calculations of solids have been limited to a single high-symmetry q-point. To the best of our knowledge, the full (q-resolved) phonon dispersion has been unaccessible so far at both VMC and DMC levels.
In this paper, we report the first successful phonon dispersion calculation of diamond at the VMC level by adopting the strategy (ii) , or the so-called frozen phonon technique 4 . The key point to success is to reduce the statistical error of atomic forces. We found that removing the nearly linear dependency of the basis set used for the trial wave function parametrization 20 dramatically lowers the statistical error of forces. Its decrease reaches the order of ∼ 10 −2 , which corresponds to a speed-up of ∼ 10 4 times in a VMC computation. This drastic reduction enables us to construct a dynamical matrix within an affordable computational time, and to eventually apply VMCphonon calculations to new interesting materials from first principles.
All Variational Monte Carlo (VMC) and lattice regularized diffusion Monte Carlo (LRDMC) 28,29 calculations in this study were performed by the TurboRVB 30 SISSA quantum Monte Carlo package. We employed the Jastrow Slaterdeterminant(JSD) ansatz, i.e, Ψ = Ψ SD × exp J, where Ψ SD and J are the Slater determinant and Jastrow terms, respectively. The Slater determinant part is expressed in terms of molecular orbitals φ k (r) = L i=1 c i,k ψ i (r) expanded in a periodized Gaussian basis set {ψ i }. The valence triple-zeta (VTZ) basis set accompanying an energy-consistent pseudo potential developed by Burkatzki et al. 31 was employed for the primitive Gaussian atomic orbitals. The coefficients of atomic orbitals (i.e., c i,k ) were obtained by a DFT calculation and were left unchanged during the VMC optimization. The Jastrow factor was composed of inhomogeneous one-, two-and threebody contributions (J = J 1 inh + J 2 + J 3 ). For details, see the Supplemental Information (SI). The variational parameters in the Jastrow terms were optimized by the stochastic reconfiguration 32 and/or the modified linear method 30,33 implemented in TurboRVB. LRDMC calculations were performed by the single-grid scheme 28 with a lattice space, a = 0.2 Bohr.
In this paper, we focus on diamond (Space group: Fm3m) as a proof of concept for the first example of a VMC-based phonon dispersion calculation. 2 × 2 × 2 conventional supercells (256 electrons / 64 carbon atoms in a simulation cell) were used for most calculations, while 3 × 3 × 3 conventional supercells (864 electrons / 216 carbon atoms in a simulation cell) were also used for several calculations to investigate the finite size errors. L-twist (i.e., k = π, π, π) was employed for alleviating the so-called one-body finite-size effects 19,34,35 . Dynamical matrices and the corresponding phonon dispersions were calculated based on the frozen-phonon technique implemented in the Phonopy module 4 , where a 0.15 Bohr displacement of carbon atoms was large enough to work with reasonable signal/noise ratio in QMC forces. This displacement underestimates harmonic frequencies only by ∼ 0.1 THz, as shown in Fig. S-3 of the supplemental information. Error bars in a phonon dispersion were estimated by the Jackknife method 36 . Phonon calculations based on DFT were performed using the Quantum Espresso package 37 with LDA-PZ 38 and GGA-PBE 39 exchange-correlation functionals at the experimental lattice parameter. The phonon dispersion of diamond has already been studied using DFT calculations by many groups so far, at the theoretical 40,41 and experimental 19 lattice parameters, that makes this system a very good testbed for any new methodological implementation of phonons calculations. As shown later, our DFT calculations are well consistent with the previous study. Fig. 1 shows the phonon dispersion obtained by our VMC calculations using the conventional 2 × 2 × 2 diamond supercell 46 . Observed experimental frequencies 21,23,26 are also plotted for comparison. A phonon dispersion obtained by the finite displacement method does not include anharmonic effects, which can decrease harmonic frequencies by up to ∼ 5-10% for the lightest elements 47 . Therefore, for comparison with the experimental results, anharmonic corrections were added to the VMC phonon dispersion in this study. The phonon dispersion before the correction is shown in Fig. S-1. ! " FIG The anharmonic renormalizations were estimated in this study using path integral molecular dynamics simulations 48 , giving −0.411 THz, −0.177 THz, and −0.280 THz for Γ, X, and L, respectively (see the SI for details). The value at Γ is consistent with a reported estimate of −17.4 cm −1 = −0.522 THz 24 . Notice that other possible sources of error were also considered. The phonon dispersion in Fig. 1 also includes the one-body finite size corrections that were estimated by DFT calculations (see Fig. S-2 in the SI for details). The two-body finite-size error is negligible because the average density and the volume do not change in phonon calculations. Tab. I shows a detailed comparison of the highest harmonic phonon frequencies at three q points, i.e., Γ = (0,0,0), X = (2π,0,0), and L = (π,π,π). Raman frequencies at Γ obtained by a direct fit of VMC energies, VMC forces, and LRDMC energies of the structures displaced along the corresponding eigenmode are also plotted in Fig. 1. 49 In Tab. I and hereinafter, (P) denotes the interpolated frequencies obtained by Phonopy, (F) denotes the phonon frequency at Γ obtained by force fitting, and (E) denotes the phonon frequency at Γ obtained by energy fitting. The corresponding Raman frequencies are consistent within the error bars, indicating that the Slater determinant obtained by DFT is almost optimal also in the presence of the Jastrow factor, because otherwise this consistency would not be possible 50 .
All VMC(P), VMC(F), and VMC(D) calculations give the harmonic phonon frequency at Γ very close to the experimental value, considering the renormalization of the anharmonic effect, i.e., the discrepancy is just ∼ 0.3 THz. On the other hand, both DFT calculations with LDA-PZ and GGA-PBE exchange-correlation functionals underestimate the highest frequency at Γ by ∼ 1.8 THz and ∼ 1.5 THz, respectively. Tab. I shows that our VMC phonon dispersion calculation also gives accurate harmonic frequencies at other q points, i.e., at X and L. Compared with the previous VMC study 19 , our VMC frequency at Γ point is closer to the experimental value, as reported in Tab. I. 51 The improvement at the VMC level certainly derives from the more flexible Jastrow factor employed in this study, while the one used in the previous VMC study was much simpler 19,52 .
It is intriguing that our LRDMC calculation gives 1.2(2) THz higher Raman frequency than the experimental one, as shown in Tab. I. The previous DMC study also overestimated the Raman frequency by 0.9(1) THz 19 . To discuss the origin of the discrepancy, we investigated the effect of the lattice-space error in our LRDMC calculations. An extrapolation (a → 0) with four lattice spaces (i.e., 0.20, 0.30, 0.40, and 0.50 Bohr) yields 41.89 (44) THz, suggesting that the lattice-space error is not the origin of the discrepancy. We also suspected that the experimental lattice parameter employed in the phonon calculation could be significantly different from the theoretical one, while this is also not the origin, as shown later. Therefore, the discrepancy should arise from the fixednode error and this should be alleviated by a nodal surface optimization 30 , which is however prohibitive in the 2 × 2 × 2 conventional supercells.
Since the equilibrium lattice parameter also affects phonon frequencies, we investigated the equation of states (EOSs) of diamond. Fig. 2 shows plots of internal energies vs. volumes fitted by the Vinet EOS 44 . Previous VMC and DMC results 19 , and the experimentally observed equilibrium lattice parameter are also plotted in Fig. 2 and summarized in Tab. II. In Fig. 2, the zero point energy (ZPE) contribution 19,53 is subtracted, to make the comparison possible with internal energies computed at T = 0 and on static lattice. Tab. II shows that our VMC calculation reproduces the previous VMC study, while our LRDMC calculation gives a slightly smaller lat- Reducing the statistical errors of atomic forces is key to a successful phonon calculation. We found that alleviating the linear dependency of a localized atomic basis set drastically decreases the statistical error. In general, a basis set optimized for molecular systems is not suitable for solid state calculations, because it is strongly linear dependent 54 . The quality of the basis set is systematically improved by a general and efficient scheme implemented in TurboRVB 20 . Indeed, the linear dependency of a localized atomic basis set (ψ µ ( r)) is characterized by the condition number, κ(S), where S µ,ν = ψ µ |ψ ν is the overlap matrix 55 . In TurboRVB, a redundant basis set is converted into a well-conditioned one, by disregarding small eigenvalues and the corresponding eigenvectors of the overlap matrix S 20,30 . Fig. 3 shows the plots of VMC energies, VMC forces, and their statistical errors vs. the inverse of the condition number (κ(S) −1 ). Fig. 3 (a) indicates that the statistical error of the energy is independent of the condition number, at variance with the statistical error of the force, which instead strongly depends on it (Fig. 3 (b)). Even with the space warp coordinate transformation (SWCT) 45 , the error bar in the forces amounts to ∼ 4.0 × 10 −2 (Ha/Bohr) when the atomic basis set is strongly linear dependent (i.e., κ(S) −1 = 10 −15 ), a condition that introduces also some bias in the forces because, as we have verified, they are no longer consistent with the finite difference energy derivatives for κ(S) −1 < 10 −12 (see Fig. 3b). On the other hand, the statistical error becomes much smaller ∼ 4.0 × 10 −4 (Ha/bohr) by removing the linear dependency (i.e., κ(S) −1 = 10 −7 ).
We analyze now in detail the reason of this behavior. Tur-boRVB evaluates atomic forces in a differential form (i.e, by the algorithmic differentiation) 45 : where R = (R 1 , . . . , R α , . . .) and ∆R α = (0, . . . , ∆R α , . . .). This equation suggests that the statistical error on forces depends on how much the wavefunction changes after an atom is displaced. In other words, to minimize the stochastic error, the overlap Ψ R+∆R α |Ψ R should be close to unity. To investigate the effect of the linear dependency on the overlap, we calculated Ψ R+∆R α |Ψ R with linear-dependent and linear-independent basis sets, using correlated sampling techniques 30 , where only one carbon atom in the 1 × 1 × 1 conventional simulation cell was displaced in the x direction by ∆R α = 0.005 Bohr. We obtained 0.9999 and 0.9726 for κ(S) −1 = 10 −7 and κ(S) −1 = 10 −15 , respectively. This clearly indicates that the linear dependency of the basis set deteriorates the overlap Ψ R+∆R α |Ψ R , thus increasing the statistical error of forces. The deterioration is explained as follows: Here, a simple Slater wavefunction without Jastrow factor is considered for the sake of clarity. In this case, the overlap Ψ R |Ψ R reads . When c i,{a,l} ≪ 1, the perturbation effect is rather small, and the orthonormalization condition almost holds, while c i,{a,l} ≫ 1 makes the perturbation effect significant, and by consequence the orthonormalization condition is certainly deteriorated. This is why the linear dependency of an atomic basis set deteriorates the overlap and, thus, induces a large error bar in forces. Thus, a complete all-electron and pseudopotential basis set database suitable for QMC solid state calculations will be quite useful for the application of QMC-forces to the calculations of phonons in realistic materials.
In summary, we report the first VMC determination of the momentum-resolved phonon dispersion of diamond. Our approach combines the ab initio quantum Monte Carlo framework with the so-called frozen phonon technique. It gives results in very good agreement with experiments and provides renormalized harmonic optical frequencies consistent with the experimental findings. We estimated the purely harmonic contribution to the phonon spectrum, by evaluating q-dependent anharmonic corrections by means of a path integral molecular dynamics driven by PBE forces. After including these corrections, the VMC phonon spectrum is significantly more accurate than the one computed by PBE. We found that alleviating the atomic basis-set redundancy of the trial wavefunction is key to reduce the statistical error of atomic forces and, thus, to make the VMC-phonons calculation feasible over the full Brillouin zone. This achievement paves the way to new relevant applications, particularly in correlated materials, where sometimes phonons are poorly reproduced within the DFT framework.  28,29 calculations in this study were performed by a SISSA quantum Monte Carlo package, called TurboRVB 30 . The package employs a manybody WF ansatz Ψ which can be written as the product of two terms, i.e., Ψ = Φ AS × exp J , where the term exp J and Φ AS are conventionally called Jastrow and antisymmetric parts, respectively. The antisymmetric part is denoted as the Antisymmetrized Geminal Power (AGP) that reads: , whereÂ is the antisymmetrization operator, and Φ r ↑ , r ↓ is called the paring function 56 . The spatial part of the geminal function is expanded over the Gaussian-type atomic orbitals: . When the paring function is expanded over M molecular orbitals where M is equal to the half of the total number of electrons (N/2), the AGP coincides with the Slater-Determinant ansatz 36,57 . In this study, we restricted ourselves to a Jastrow-Slater determinant (JSD) by setting M = 1 2 · N, wherein the coefficients of atomic orbitals, i.e., c i,k , were obtained by a Density Functional theory (DFT) calculation, and were fixed during a VMC optimization.
The Jastrow term is composed of one-body, two-body   . S-1 and three/four-body factors (J = J 1 J 2 J 3/4 ). The onebody and two-body factors are essentially used to fulfill the electron-ion and electron-electron cusp conditions, respectively, and the three/four-body factor is employed to consider further electron-electron correlations (e.g., electron-nucleuselectron). Since we employed an energy-consistent pseudo potential developed by Burkatzki et al. 31 , we used only the inhomogeneous part of J inh 1 , J 2 and J 3/4 in this study (i.e., no electron-ion cusp corrections are needed). The inhomogeneous one-body Jastrow factor J inh 1 reads J inh 1 (r 1 , . . . , N atom a=1 l M a,l χ a,l (r i ) , where r i are the electron positions, R a are the atomic positions with corresponding atomic number Z a , l runs over atomic orbitals χ a,l (e.g., GTO) centered on the atom a, N atom is the total number of atoms in a system, and {M a,l } are variational parameters. The two-body Jastrow factor is defined as: and F is a variational parameter. The three-body Jastrow factor is: J 3/4 (r 1 , . . . r N ) = exp i< j Φ Jas r i , r j , and Φ Jas r i , r j = l,m,a,b g a,l,m,b χ Jas a,l (r i ) χ Jas b,m r j , where the indices l and m again indicate different orbitals centered on corresponding atoms a and b. In this study, the coefficients of the three/four-body Jastrow factor were set to zero for a b because it significantly decreases the number of variational parameters while rarely affects variational energies . A basis set (3s2p) was employed for the atomic orbitals of the Jastrow part. The variational parameters in the Jastrow factor were optimized by the so-called stochastic reconfiguration 32 and/or the modified linear method 30,33 implemented in TurboRVB. Using the optimized wavefunction, energies and forces are calculated at the VMC and the LRDMC levels. All LRDMC calculations were performed by the original single-grid scheme 28 with the discretization grid size a = 0.20 Bohr. For the Raman frequency calculation, several lattice spaces, i.e., a = 0.20, 0.30, 0.35, 0.40, and 0.50 Bohr were also used to extrapolate the energies (i.e., a → 0) with a quadratic function (i.e., E(a) = E 0 + k 1 · a 2 + k 2 · a 4 ). Figure S-1 shows the phonon dispersions and Raman frequencies of diamond calculated using 2 × 2 × 2 conventional supercell at the VMC level. "Raw data" denotes the phonon dispersion obtained by the finite-displacement method, while "Corrected" is the phonon dispersion where the one-body finite-size (see sec. B) and anharmonic (see sec. C) corrections are included. There are two major implementations to compute a harmonic phonon dispersion, i.e., the finite displacement method FIG. S-3. Phonon frequencies calculated by DFT at q = Γ and L points v.s. employed displacements in the finite displacement method, where the GGA-PBE functional was employed. The phonon dispersions were calculated using 2 × 2 × 2 conventional supercell with k twist = L (π,π,π).
(FDM) and the linear response method. The former is also called frozen phonon method. Their results should be consistent as far as employed hyperparameters are correct.  58 were employed for k and q integrations. The DFT calculations were performed with GGA-PBE functionals with Ultrasoft (US) pseudo potential taken from PS-Library 59 . Comparison of (a)-red with (c)-green proves that the one-body finite size error is very small in a phonon calculation thanks to the error cancellation when L twist for 2 × 2 × 2 conventional supercell is employed. In detail, the one-body finite size errors are − 0.16 THz, − 0.18 THz, and − 0.23 THz for q = Γ, X, and L, respectively. Comparison of (b)-blue with (c)green, which are almost overlapped, suggests that the finite difference method does not bias the result. As mentioned in the main text, 0.15 Bohr was employed for the displacement of a carbon atom in the VMC frozen-phonon calculation to decrease the statistical error. Figure " FIG. S-4. Plots of Raman Frequencies vs. lattice parameters of diamond. They were calculated using the linear response theory implemented in Quantum Espresso at a single q point, i.e., q = (0,0,0) with the LDA-PZ exchange-correlation functional. The primitive lattice was employed with k-grid = 8 × 8 × 8.
this underestimation as mentioned in the main text. Fig. S-4 shows a plot of Raman Frequencies vs. lattice parameters of diamond. They were calculated using the linear response theory implemented in Quantum Espresso at a single q point, i.e., q = (0,0,0) with the LDA-PZ exchange-correlation functional. The primitive lattice was employed with 8 × 8 × 8.
Appendix C: Anharmonic renormalizations from path integral molecular dynamics The anharmonic renormalizations of phonon frequencies reported in Table I with the superscript ' f ' were computed through the displacement-displacement zero-time Kubo-correlator built upon path integral molecular dynamics (PIMD) trajectories. The evaluation of the displacementdisplacement correlator over classical molecular dynamics trajectories is a standard method to predict force constant matrices and phonon spectra, but it misses nuclear quantum effects. The details of the extension of such a method to extract phonon dispersions from path integral molecular dynamics will be published soon in another paper 48 .
In Fig S-5, we report the whole spectrum that we have obtained using this new approach compared to the linear response DFT curve. In particular, the ion equations of motion were integrated using the PIOUD algorithm 60 , while at each time step the electronic potential energy surface was evaluated at DFT level using the Quantum Espresso package 37 and GGA-PBE 39 exchange-correlation functionals. The simulation was carried out at 300 Kelvin for 34 picoseconds employing 12 beads. We used the same 64-atoms supercell as for the QMC calculations, with a 2 × 2 × 2 k-grid for electronic integration and a cutoff for wavefunctions equal to 60 Ry.

Appendix D: Finite-size errors in QMC calculations
The effect of the finite-size errors on an EOS calculation has been investigated. Figure S-6 shows that finite-size extrapolations for two different lattice parameters (a = 6.6107 and 6.7407 Bohr). The differences of the energies are constant regardless of the simulation cell size, indicating that the finite-size error is negligible for an EOS calculation thanks to the error cancellation. Since diamond is an insulator, raw QMC energies obtained at L point can be smoothly extrapolated by a linear function. See. Ref. 34.

Appendix E: Vinet exponential function
In this study, the obtained energies were fitted by the Vinet exponential function 44 that reads: where E(V) is the total energy per atom. V is volume per atom, and E 0 , V 0 , B 0 , and B ′ 0 are parameters. The obtained stresses were fitted by its derivate form, i.e., where P(V) is the pressure. V is volume per atom, and V 0 , B 0 , and B ′ 0 are parameters. The non-linear fittings were performed using SciPy module 61 implemented in Python.

Appendix F: Atomic forces
TurboRVB evaluates atomic forces in a differential expression 45 : where J is the Jacobian of the space warp coordinate transformation (SWCT) 45 , e L is the local energy, and the brackets indicates a Monte Carlo average over the trial WF. The first term is called "Hellmann-Feynman" force (HF) and the second and third terms are called "Pulay" force (PF). All the terms above can be written by the partial derivatives of the local energy and those of the logarithm of the WF 45 . These differential expressions are efficiently computed in TurboRVB, by using the adjoint algorithmic differentiation technique 45 . TurboRVB also employs the so-called reweighting technique developed by Attaccalite and Sorella 62 to avoid divergence of the derivatives in the vicinity of the nodal surfaces. Fig. S-7 shows VMC forces, its Hellmann-Feynman, and its Pulay contributions v.s. inverse of the condition number of the overlap matrix, corresponding to Fig. 3 in the main text. Fig. S-7 clarifies that the erratic behavior of the forces at κ(S) −1 = 10 −12 comes from its Pulay contribution, indicating that the linear-dependency of a basis-set also deteriorates absolute values of forces. FIG. S-7. VMC forces, its Hellmann-Feynman, and its Pulay contributions v.s. inverse of the condition number of the overlap matrix. The forces were calculated with the space warp coordinate transformations (SWCT) 45 . The same VTZ atomic basis set and the pseudo potential 31 as in the phonon and EOS calculations were adopted, while 1 × 1 × 1 conventional cell (8 carbon atoms in a simulation cell) with k = Γ twist was employed. Only one carbon atom was displaced in x direction by 0.15 Bohr.