*Ab initio* molecular dynamics with quantum Monte Carlo

- International School for Advanced Studies (SISSA) and CRS Democritos, CNR – Istituto Nazionale per la Fisica della Materia (INFM), Trieste, Italy

Quantum Monte Carlo (QMC) is one of the most accurate electronic structure methods for *ab initio* many-body calculations in a broad range of electronic systems. Thanks to fast evolving massively parallel computers, much larger simulations of realistic systems become affordable. We have developed an *ab initio* simulation scheme at finite temperature based on molecular dynamics (MD) and QMC during the past few years. In this approach, a second order Langevin MD is employed by using a statistical evaluation of the forces acting on each atom by means of QMC. Moreover, the corresponding statistical noise is also used to drive and accelerate the dynamics of ions. The accuracy and the reliability of our *ab initio* MD have been studied systematically and it has been already successfully applied to study the molecular and atomic phases of liquid hydrogen under pressure and liquid water at ambient conditions. Since this *ab initio* method provides a better resolution of the thermodynamic properties of materials, we believe that it represents a very promising tool for the most challenging applications in physics, chemistry, and biology.

Since the Car-Parrinello molecular dynamics (Car and Parrinello, 1985) was proposed three decades ago, *ab initio* molecular dynamics (MD) based on the density functional theory (DFT) has been widely accepted as an accurate and powerful tool to study the thermodynamic properties of systems at ambient and extreme conditions in the fields of physics, chemistry, and biology. Even though DFT is in principle exact, its accuracy, in practice, is affected by the quality of the approximation used in the exchange correlation functional term whose universal exact expression is not accessible, as the simplest approximate functionals, such as LDA, PBE, BLYP, and several other ones, so far cannot be systematically improved in a computationally efficient way. Quantum Monte Carlo (QMC) is the most promising successor of DFT in *ab initio* molecular dynamics. Indeed, the scaling behavior of QMC with the number of electrons is very good, from the second power (Ceperley et al., 1977) for thermodynamic properties – e.g., energy per electron, because in this case the number of samples to obtain a fixed accuracy in this quantity can be decreased when increasing the number *N*_{e} of electrons, implying ${N}_{\mathrm{e}}^{3}\u2215{N}_{\mathrm{e}}={N}_{\mathrm{e}}^{2}$ scaling with the standard algorithm – to the fourth power of the number of electrons when chemical accuracy in the total energy is required – in this case, instead the number of samples required for fixed accuracy has to increase with *N*_{e}, yielding ${N}_{\mathrm{e}}^{3}\times {N}_{\mathrm{e}}={N}_{\mathrm{e}}^{4}$ scaling. For this reason (a cost similar to DFT, albeit with a considerably larger pre-factor), QMC is also by far more efficient compared with other quantum chemistry methods based on multi-determinant expansions or many-body perturbation theory around the Hartree–Fock, such as, for instance, configuration interaction and coupled cluster theory, respectively. These methods provide a similar or even better accuracy than QMC for systems consisting of a few atoms, but they become computationally prohibitive when the system sizes increase. Another very important reason favoring the QMC algorithm is its well established efficient parallelization in current and upcoming supercomputers. QMC codes can be easily adapted to millions of CPU cores while DFT codes are still struggling with its parallel performance over several thousand cores. For the above reasons, it is probably not too far from reality to claim that QMC represents the rising star in the *ab initio* MD simulations.

This opinion paper focuses on the construction of the QMC-based MD by showing the key ingredients and the main recent developments of the building blocks necessary for the efficient implementation of this method. There are four main components: the basics of QMC, the wavefunction optimization, the evaluation of forces, and the molecular dynamics scheme. The last three ingredients were efficiently developed within only 10 years, with substantial contribution given also by us (Casula et al., 2004; Sorella et al., 2007; Umrigar et al., 2007; Attaccalite and Sorella, 2008; Sorella and Capriotti, 2010; Luo et al., 2014; Mazzola et al., 2014; Zen et al., 2014). Thus, the QMC-based MD is no longer a dream.

## Basics of QMC

In QMC, the wavefunction ansatz is of fundamental importance for the efficiency, reliability, and accuracy of both variational and diffusion Monte Carlo calculations (VMC/DMC), namely the two most popular methods in QMC, the latter being more accurate and based on an approximate evaluation of the ground state properties with the restriction of having the same phases of the best reference VMC wavefunction – i.e., the so-called fixed-node approximation (Anderson, 1976; Ceperley and Alder, 1984; Umrigar et al., 1993), which is unavoidable for fermionic systems. The repeated evaluation of the wavefunction with also its gradients and Laplacian is one of the most consuming parts of the computation; the optimal wavefunction ansatz for MD should be also compact. The most popular choice is the Jastrow factor (Jastrow, 1955) correlated Slater determinant, which is simple and generally accurate. Pseudopotentials (Trail and Needs, 2005; Burkatzki et al., 2007) are often employed like in DFT to replace chemically in-active core electrons and keep only the valence electrons in QMC calculations.

## Wavefunction Optimization

In the MD simulation, after each ionic move, the parameters in the VMC wavefunction need to be optimized in order to satisfy the Born–Oppenheimer approximation. Efficient optimization methods ensure that the energy converges to the right minimum within a sufficient accuracy in just a few steps. At present, the variance minimization method (Umrigar et al., 1988) has been recently abandoned because sometimes it failed to yield the best possible wavefunction (Snajdr and Rothstein, 2000; Bressanini et al., 2002). This method has been replaced by much more efficient and robust energy minimization methods. In our previous works, we introduced an energy minimization method, called stochastic reconfiguration (SR), which is based only on energy first derivatives, and allows a much more efficient optimization compared with other more conventional methods of this type, such as the steepest decent. In addition, the trial step length can be automatically determined by the SR with the Hessian accelerator (SRH) (Sorella, 2005) method, which makes use of the second derivatives of the energy. This method can be further improved by the recently introduced “linear method” (Umrigar and Filippi, 2005; Umrigar et al., 2007), which provides a better sampling of the Hessian matrix and represents at present the state of the art for the QMC optimization of wavefunction in realistic systems. For the application of the above techniques, to large systems, it is important to perform in a very efficient way the expensive inversion operation of a big matrix with its leading dimension (10^{4} is possible on current supercomputers) equal to the number of variational parameters (Sorella et al., 2007; Neuscamman et al., 2012).

## Evaluation of Forces

In order to study thermodynamic properties, classical Monte Carlo and molecular dynamics are both widely used simulation methods. The former has a significant advantage that it does not necessarily require the computation of forces. Unfortunately, it is not straightforward to combine classical Monte Carlo and quantum Monte Carlo because all the quantities evaluated by QMC, including the energy, carry statistical errors that are very expensive to eliminate.^{1} In order to move the ions according to the QMC energy and sample the canonical ensemble correctly, the penalty method (Ceperley and Dewing, 1999) was introduced to remove the bias due to the statistical uncertainty in the knowledge of the energy by rejecting the proposed moves more frequently than the standard Metropolis algorithm. For this reason, this method is very expensive, particularly in the low temperature regime, because of too many rejected moves, and so far applications have been limited to hydrogen with up to 54 protons in this regime (Pierleoni and Ceperley, 2006; Morales et al., 2010a,b).

Generally speaking, it is clear that, when the computational cost for the calculation of the nuclear forces is comparable to that of the energy, MD should be more efficient, because with the same cost all the atoms are moved without any rejection. For instance, in DFT, where the forces are obtained almost for free by applying the Hellmann–Feynman theorem, MD is a common practice to sample the canonical distribution, and, to our knowledge, only hybrid methods based on Monte Carlo and molecular dynamics (Rossky et al., 1978; Duane et al., 1987) can be competitive.

Unfortunately, computing forces with QMC is substantially more demanding than DFT. To obtain the forces in an accurate and efficient way within the VMC framework, several very important improvements have been done, concerning both theories and algorithms. The correlated sampling (Filippi and Umrigar, 2000) method allows the computation of the energy derivatives with errors much smaller than those obtained with a straightforward finite difference method. Together with the space warp coordinate transformation (Umrigar, 1989), forces computed by QMC were used for structural optimization (Casula et al., 2004). Another recent and important development in QMC was the solutions (Assaraf and Caffarel, 2000, 2003; Attaccalite and Sorella, 2008; Sorella and Capriotti, 2010; Zen et al., 2013) of the intrinsic infinite variance problem occurring in the calculation of nuclear forces in the simplest VMC scheme. Moreover, thanks to the algorithmic differentiation algorithm (Sorella and Capriotti, 2010), the cost of computing all the force components in a system can be afforded with a computational time at most a factor four larger than the one corresponding to the energy. This progress has led to several works, where structural optimization and highly accurate evaluations of the equilibrium configurations as well as related properties were possible even for quite large systems (Barborini and Guidoni, 2012; Coccia et al., 2013, 2014; Guareschi and Filippi, 2013; Zen et al., 2013).

## Molecular Dynamics with QMC

Even though the efficient wavefunction optimization and fast evaluation of forces are achieved, the application of QMC for MD simulation of materials remains challenging due to the theoretical difficulties in applying the Newton’s equations of motion when the forces are given with a statistical uncertainty. For instance, the basic law of energy conservation cannot be met at all, when the forces are not exactly given at each step. Since it is prohibitive to eliminate the statistical error (for the same reason, see footnote 1), MD methods that allow the presence of random noise in the force components are certainly more suitable for QMC.

In our works, we have adopted a second order Langevin dynamics (SLD) for the sampling of the ionic configurations within the ground state Born–Oppenheimer approach. In this scheme, the statistical noise corresponding to the forces, is used to drive the dynamics at finite temperature by means of an appropriate generalized Langevin dynamics (Attaccalite and Sorella, 2008), while, during the dynamics, a noise correction is added to compensate the extra energy pumped into the system by the noise of the QMC forces. Later works also introduced two integration methods (Luo et al., 2014; Mazzola et al., 2014), which allowed us to work with noisy forces and large time steps, thus reducing the computational cost by a substantial fraction. The price to pay is the limitation of this approach to the evaluation of static quantities at finite temperature, because dynamical quantities, such as the diffusion constant, cannot be evaluated with the above stochastic implementation of SLD. Similar approaches (Krajewski and Parrinello, 2006; Kühne et al., 2007) have been proposed where an SLD algorithm has been devised also at the DFT level.

Last but not least, our proposed MD scheme with QMC also includes an automatic optimal tuning for the damping of the slow and the fast modes, which greatly improves the efficiency of the sampling. In the SLD equations, the damping coefficients matrix can be arbitrarily chosen to speed up the sampling, provided the fluctuation–dissipation theorem is satisfied by adding appropriate correlated noise to the force components. In our scheme, we have defined this correlation matrix in terms of the one corresponding to the noisy forces obtained with QMC. This is because we have empirically verified that the latter matrix is almost proportional to the dynamical matrix (see Figure 1). Therefore, the energy modes with higher frequency vibrations can be systematically damped much more than the slow modes, and this clearly allows a faster propagation with larger time steps.

**Figure 1. The eigenvectors of the covariance matrix obtained with the noisy force components evaluated by VMC**. These eigenvectors correspond to the three vibrational modes of the water monomer: bending (red), symmetrical (green), and asymmetrical (blue) stretching. The smallest eigenvalue corresponds to the lowest frequency vibrational mode (red).

The VMC-based MD has been successfully applied to various systems from small molecules (Luo et al., 2014) to liquids, for instance, the liquid hydrogen at high pressure (Attaccalite and Sorella, 2008; Mazzola et al., 2014) and the liquid water at ambient conditions (Zen et al., 2014) for the properties of vibrational frequencies (Luo et al., 2014), atomic structures (Zen et al., 2014), and phase diagrams (Mazzola et al., 2014). It is implemented in the `TurboRVB`

QMC package (Sorella, 2015) with a hybrid MPI and OpenMP paradigm and its performance in the production run reaches peta-scale, as it has been successfully tested up to 33K cores on BG/Q and 264K cores on the K-computer.

## Conclusion and Outlooks

From the high performance computing perspective, the four main components that we have briefly reviewed in this work show different features in computational intensity and communications. The statistical calculation for the energy and its derivatives required by the optimization and forces puts very heavy loads on the CPU but only needs very little communication at the end of a big piece of work, which is an ideal pattern for the parallelization. In addition, each optimization step involves inversion operations of big matrices. It can be a potential bottleneck if it is solved in parallel on a network with high latency. The computational and communication cost of updating ion positions is negligible compared to the other three parts for tens and hundreds of atoms.

In conclusion, we have briefly described the main computational achievements that have allowed recently the first *ab initio* simulations of realistic materials with an efficient molecular dynamics based on accurate QMC forces. We believe that our work represents a first step toward an efficient tool for *ab initio* simulations of material properties much more reliable than DFT. Indeed, our evaluation of forces is limited to the VMC approach, and we expect some progress can be done by adopting the more reliable DMC method. Similarly, it should be possible to extend the method to more accurate wavefunctions, containing, for instance, many-body and backflow correlations, as well as implicit multi-determinant wavefunctions, such as the antisymmetrized geminal power (AGP) and the Pfaffian wavefunction. For all the above reasons, we remark once again our belief that QMC should be considered as the raising star for *ab initio* simulations of materials.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Footnote

**^**The statistical error decays as $1\u2215\sqrt{N}$, where*N*is the number of uncorrelated samples.

## References

Anderson, J. B. (1976). Quantum chemistry by random walk. *J. Chem. Phys.* 65, 4121–4127. doi: 10.1063/1.432868

Assaraf, R., and Caffarel, M. (2000). Computing forces with quantum Monte Carlo. *J. Chem. Phys.* 113, 4028–4034. doi:10.1063/1.1286598

Assaraf, R., and Caffarel, M. (2003). Zero-variance zero-bias principle for observables in quantum Monte Carlo: application to forces. *J. Chem. Phys.* 119, 10536–10552. doi:10.1063/1.1621615

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Attaccalite, C., and Sorella, S. (2008). Stable liquid hydrogen at high pressure by a novel ab initio molecular-dynamics calculation. *Phys. Rev. Lett.* 100, 114501. doi:10.1103/PhysRevLett.100.114501

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Barborini, M., and Guidoni, L. (2012). Reaction pathways by quantum Monte Carlo: insight on the torsion barrier of 1,3-butadiene, and the conrotatory ring opening of cyclobutene. *J. Chem. Phys.* 137, 224309. doi:10.1063/1.4769791

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Bressanini, D., Morosi, G., and Mella, M. (2002). Robust wave function optimization procedures in quantum Monte Carlo methods. *J. Chem. Phys.* 116, 5345–5350. doi:10.1063/1.1455618

Burkatzki, M., Filippi, C., and Dolg, M. (2007). Energy-consistent pseudopotentials for quantum Monte Carlo calculations. *J. Chem. Phys.* 126, 234105. doi:10.1063/1.2741534

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Car, R., and Parrinello, M. (1985). Unified approach for molecular-dynamics and density-functional theory. *Phys. Rev. Lett.* 55, 2471–2474. doi:10.1103/PhysRevLett.55.2471

Casula, M., Attaccalite, C., and Sorella, S. (2004). Correlated geminal wave function for molecules: an efficient resonating valence bond approach. *J. Chem. Phys.* 121, 7110–7126. doi:10.1063/1.1794632

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Ceperley, D., Chester, G. V., and Kalos, M. H. (1977). Monte Carlo simulation of a many-fermion study. *Phys. Rev. B* 16, 3081–3099. doi:10.1103/PhysRevB.16.3081

Ceperley, D., and Dewing, M. (1999). The penalty method for random walks with uncertain energies. *J. Chem. Phys.* 110, 9812–9820. doi:10.1063/1.478034

Ceperley, D. M., and Alder, B. J. (1984). Quantum Monte Carlo for molecules: Green’s function and nodal release. *J. Chem. Phys.* 81, 5833–5844. doi:10.1063/1.447637

Coccia, E., Varsano, D., and Guidoni, L. (2013). Protein field effect on the dark state of 11-cis retinal in rhodopsin by quantum Monte Carlo/molecular mechanics. *J. Chem. Theory Comput.* 9, 8–12. doi:10.1021/ct3007502

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Coccia, E., Varsano, D., and Guidoni, L. (2014). Ab initio geometry and bright excitation of carotenoids: quantum Monte Carlo and many body green’s function theory calculations on peridinin. *J. Chem. Theory Comput.* 10, 501–506. doi:10.1021/ct400943a

Duane, S., Kennedy, A., Pendleton, B., and Roweth, D. (1987). Hybrid Monte Carlo. *Phys. Lett. B* 195, 216–222. doi:10.1016/0370-2693(87)91197-X

Filippi, C., and Umrigar, C. J. (2000). Correlated sampling in quantum Monte Carlo: a route to forces. *Phys. Rev. B* 61, R16291–R16294. doi:10.1103/PhysRevB.61.R16291

Guareschi, R., and Filippi, C. (2013). Ground- and excited-state geometry optimization of small organic molecules with quantum Monte Carlo. *J. Chem. Theory Comput.* 9, 5513–5525. doi:10.1021/ct400876y

Jastrow, R. (1955). Many-body problem with strong forces. *Phys. Rev.* 98, 1479–1484. doi:10.1103/PhysRev.98.1479

Krajewski, F., and Parrinello, M. (2006). Linear scaling electronic structure calculations and accurate statistical mechanics sampling with noisy forces. *Phys. Rev. B* 73, 041105. doi:10.1103/PhysRevB.73.041105

Kühne, T., Krack, M., Mohamed, F., and Parrinello, M. (2007). Efficient and accurate Car-Parrinello-like approach to Born-Oppenheimer molecular dynamics. *Phys. Rev. Lett.* 98, 066401. doi:10.1103/PhysRevLett.98.066401

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Luo, Y., Zen, A., and Sorella, S. (2014). Ab-initio molecular dynamics with noisy forces: validating the quantum Monte Carlo approach with benchmark calculations of molecular vibrational properties. *J. Chem. Phys.* 141, 194112. doi:10.1063/1.4901430

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Mazzola, G., Yunoki, S., and Sorella, S. (2014). Unexpectedly high pressure for molecular dissociation in liquid hydrogen by electronic simulation. *Nat. Commun.* 5, 3487. doi:10.1038/ncomms4487

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Morales, M. A., Pierleoni, C., and Ceperley, D. M. (2010a). Equation of state of metallic hydrogen from coupled electron-ion Monte Carlo simulations. *Phys. Rev. E* 81, 021202. doi:10.1103/PhysRevE.81.021202

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Morales, M. A., Pierleoni, C., Schwegler, E., and Ceperley, D. M. (2010b). Evidence for a first-order liquid-liquid transition in high-pressure hydrogen from ab initio simulations. *Proc. Natl. Acad. Sci. U.S.A.* 107, 12799–12803. doi:10.1073/pnas.1007309107

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Neuscamman, E., Umrigar, C. J., and Chan, G. K.-L. (2012). Optimizing large parameter sets in variational quantum Monte Carlo. *Phys. Rev. B* 85, 045103. doi:10.1103/PhysRevB.85.045103

Pierleoni, C., and Ceperley, D. M. (2006). Computer simulations in condensed matter systems: from materials to chemical biology volume 1. *Lect. Notes Phys.* 703, 641–683. doi:10.1007/3-540-35273-2_18

Rossky, P. J., Doll, J. D., and Friedman, H. L. (1978). Brownian dynamics as smart Monte Carlo simulation. *J. Chem. Phys.* 69, 4628–4633. doi:10.1063/1.436415

Snajdr, M., and Rothstein, S. M. (2000). Are properties derived from variance-optimized wave functions generally more accurate? Monte Carlo study of non-energy-related properties of h2, he, and lih. *J. Chem. Phys.* 112, 4935–4941. doi:10.1063/1.481047

Sorella, S. (2005). Wave function optimization in the variational Monte Carlo method. *Phys. Rev. B* 71, 241103. doi:10.1103/PhysRevB.71.241103

Sorella, S. (2015). TurboRVB quantum Monte Carlo package. Available at: http://people.sissa.it/~sorella/web/index.html

Sorella, S., and Capriotti, L. (2010). Algorithmic differentiation and the calculation of forces by quantum Monte Carlo. *J. Chem. Phys.* 133, 234111. doi:10.1063/1.3516208

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Sorella, S., Casula, M., and Rocca, D. (2007). Weak binding between two aromatic rings: feeling the van der Waals attraction by quantum Monte Carlo methods. *J. Chem. Phys.* 127, 014105. doi:10.1063/1.2746035

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Trail, J. R., and Needs, R. J. (2005). Norm-conserving hartree-fock pseudopotentials and their asymptotic behavior. *J. Chem. Phys.* 122, 14112. doi:10.1063/1.1829049

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Umrigar, C. J. (1989). Two aspects of quantum Monte Carlo: determination of accurate wavefunctions and determination of potential energy surfaces of molecules. *Int. J. Quantum Chem.* 36, 217–230. doi:10.1002/qua.560360826

Umrigar, C. J., and Filippi, C. (2005). Energy and variance optimization of many-body wave functions. *Phys. Rev. Lett.* 94, 150201. doi:10.1103/PhysRevLett.94.150201

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Umrigar, C. J., Nightingale, M. P., and Runge, K. J. (1993). A diffusion Monte Carlo algorithm with very small time-step errors. *J. Chem. Phys.* 99, 2865–2890. doi:10.1063/1.465195

Umrigar, C. J., Toulouse, J., Filippi, C., Sorella, S., and Rhenning, H. (2007). Alleviation of the fermion-sign problem by optimization of many-body wave functions. *Phys. Rev. Lett.* 98, 110201. doi:10.1103/PhysRevLett.98.110201

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Umrigar, C. J., Wilson, K. G., and Wilkins, J. W. (1988). Optimized trial wave functions for quantum Monte Carlo calculations. *Phys. Rev. Lett.* 60, 1719–1722. doi:10.1103/PhysRevLett.60.1719

Zen, A., Luo, Y., Mazzola, G., Guidoni, L., and Sorella, S. (2014). Ab-initio molecular dynamics simulation of liquid water by quantum Monte Carlo. *J. Chem. Phys.* arXiv: 1412.2936.

Zen, A., Luo, Y., Sorella, S., and Guidoni, L. (2013). Molecular properties by quantum Monte Carlo: an investigation on the role of the wave function ansatz and the basis set in the water molecule. *J. Chem. Theory Comput.* 9, 4332–4350. doi:10.1021/ct400382m

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Keywords: quantum Monte Carlo, *ab initio* calculations, molecular dynamics, many-body physics, high performance computing

Citation: Luo Y and Sorella S (2015) *Ab initio* molecular dynamics with quantum Monte Carlo. *Front. Mater.* **2**:29. doi: 10.3389/fmats.2015.00029

Received: 11 February 2015; Paper pending published: 26 February 2015;

Accepted: 22 March 2015; Published online: 07 April 2015.

Edited by:

Simone Taioli, Bruno Kessler Foundation, Italy; Charles University, Czech RepublicCopyright: © 2015 Luo and Sorella. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: xw111luoye@gmail.com