Free energy methods in Coupled Electron Ion Monte Carlo
Abstract
Recent progress in simulation methodologies and in computer power allow first principle simulations of condensed systems with BornOppenheimer electronic energies obtained by Quantum Monte Carlo methods. Computing free energies and therefore getting a quantitative determination of phase diagrams is one step more demanding in terms of computer resources. In this paper we derive a general relation to compute the free energy of an abinitio model with Reptation Quantum Monte Carlo (RQMC) energies from the knowledge of the free energy of the same abinitio model in which the electronic energies are computed by the less demanding but less accurate Variational Monte Carlo (VMC) method. Moreover we devise a procedure to correct transition lines based on the use of the new relation. In order to illustrate the procedure, we consider the liquidliquid phase transition in hydrogen, a first order transition between a lower pressure, molecular and insulating phase and a higher pressure, partially dissociated and conducting phase. We provide new results along the isotherm across the phase transition and find good agreement between the transition pressure and specific volumes at coexistence for the model with RQMC accuracy between the prediction of our procedure and the values that can be directly inferred from the observed plateau in the pressurevolume curve along the isotherm. This work paves the way for future use of VMC in first principle simulations of high pressure hydrogen, an essential simplification when considering larger system sizes or quantum proton effects by Path Integral Monte Carlo methods.
onte Carlo Methods; Quantum Monte Carlo; High pressure hydrogen, Free energies and phase diagrams
1 Introduction
Firstprinciples simulation has become an essential method to investigate the physical behavior of condensed matter systems, in particular when the chemical nature of the interaction among the nuclear degrees of freedom depends on the external conditions, like for instance in systems under extreme conditions of pressure and temperature. Those methods generally treat systems of classical nuclei at finite temperature and electrons in their instantaneous ground state in the framework of the BornOppenheimer (BO) approximation, although extensions to treat quantum nuclei at finite temperature within the Path Integral formalism, and/or electron at finite temperature with fractional (fermi) occupation have been developed. The electronic structure at fixed nuclear positions is generally computed within the Density Functional Theory (DFT) framework [1] or, when higher accuracy is required, by Quantum Monte Carlo methods (QMC) [2]. DFT is, in principle, an exact theory but, in practice, it is based on uncontrolled approximations for the exchangecorrelation functional. Despite some well documented limitations, DFT is often accurate and fast enough to be used in conjunction with Molecular Dynamics to perform dynamical studies for systems of several hundred nuclei and extract physical information (firstprinciples Molecular Dynamics, FPMD) [3, 4, 5]. Ground state QMC comes in two different flavors: Variational Monte Carlo (VMC) and Projection Monte Carlo (either Diffusion Monte Carlo (DMC) or Reptation Quantum Monte Carlo (RQMC)). VMC is faster but the accuracy of its predictions is limited by the skill in designing and optimizing suitable manybody trial wave functions. RQMC is a method to automatically optimize trial wave functions; it provides the most accurate description of the electronic properties but has larger computational requirements (about an order of magnitude larger than VMC, not counting the VMC optimization process). It is an exact method for bosons since it is able to provide the properties of the bosonic ground state in a time which is polynomial in the number of particles. In the case of fermions (electrons) the scaling is exponential and for practical calculations we have to resort to the ”fixed node approximation”. This approximation makes the method variational with respect to the trial nodes; usually the results remain highly accurate if good trial nodal surfaces are employed [2]. With recent progress in optimization procedures [6] the location of the trial nodal surfaces, suitably parametrized, could also be optimized to achieve higher accuracy.
In recent years we have developed an abinitio method in which electronic energies are computed by ground state Quantum Monte Carlo (QMC) methods while nuclear degrees of freedoms are sampled by Metropolis Monte Carlo. This method, called Coupled ElectronIon Monte Carlo (CEIMC) [7], has been successfully applied to high pressure hydrogen. The results are a benchmark of DFT calculations in the region of phase diagram of interest for planetary physics [8, 9] and across the metalinsulator transition region in the fluid phase [10, 11].
The general picture emerging is that FPMD for high pressure hydrogen is generally accurate away from the metalinsulator transition while QMC accuracy is required when approaching the metallization of the system.
In CEIMC we sample the equilibrium distribution for the protonic degrees of freedom by a Metropolis based Monte Carlo scheme employing either VMC or RQMC to compute the difference between the BO energies of two nearby protonic configurations in the Metropolis method. The electronic calculation needs to be repeated for each attempted move of the protons and it represents by far the largest part of the computational load of the method. It is therefore very tempting to exploit VMC rather than RQMC in CEIMC. On the other hand if the trial wave function is not accurate enough, the use of VMC could provide biased results which can lead to inaccurate predictions of the physical behavior of the system [10].
The key quantity in tracing transition lines is the free energy difference (either Helmholtz or Gibbs) between coexisting phases which is usually obtained by computing the absolute free energies of the coexisting phases. This can be achieved by the Coupling Constant Integration (CCI) method, a general strategy introduced by Kirkwood [12] to adiabatically transform the system into another system of known free energy. This strategy is often used in abinitio simulations to compute the free energy of the abinitio model starting from the knowledge of the free energy of an effective classical representation of the same system, for which standard free energy methods can be easily applied [13, 14, 15, 16, 11].
In this paper, within the framework of the CCI and of the CEIMC method, we derive a relation which allows us to obtain the free energy of the system with RQMC electronic energies from the knowledge of the free energy of the system with VMC energies. The simple relation, derived in the next section, is of practical relevance when tracing transition lines since it allows us to obtain the transition lines with RQMC accuracy by performing free energy calculations with VMC electronic energies. This strategy improves the efficiency of the method by roughly times, depending on the accuracy of the trial wave function.
To illustrate the use of the new relation we consider the liquidliquid phase transition in high pressure hydrogen, predicted earlier [17, 18, 19] using the socalled chemical model, but only recently put on firm basis by firstprinciple simulations [11, 15]. This is a first order transition between an insulating, molecular liquid phase at lower pressure and a metallic, partially dissociated liquid phase at higher pressure. It is predicted to occur below 2000K; on increasing pressure the specific volume is discontinuous, corresponding to a partial dissociation of molecules and to an abrupt jump of roughly four orders of magnitude of the DC conductivity [11]. Above T=2000K the process is observed to become continuous and a critical point is predicted to be located at and , . The transition line T(P) has a negative slope, meaning that lowering the temperature requires higher pressure to cross the transition, and it is predicted to meet the melting line of the molecular insulating crystalline structure (Phase I) at a triple point located at and within the DFT framework and and using CEIMC [11]. The latter prediction was based on an extrapolation to high pressure of lower pressure results.
In this paper we performed new CEIMC calculations of high pressure hydrogen for T=600K and across the liquidliquid phase transition. We performed calculations using both VMC and RQMC with the aim of testing the accuracy of VMC across this particular phase transition. Having at our disposal the equation of state obtained by the two methods we can test the reliability of our new relation for correcting free energies.
The paper is organized as follows. In the next section we will set the formal framework of RQMC and we will derive the formula to correct the free energy. In section 3, as an application of our strategy, we will present recent results obtained for the liquidliquid transition line of high pressure hydrogen. Finally in section 4, we will draw some conclusions.
2 Theoretical framework
Let us consider a system of monovalent ions of mass and electrons of mass in a volume and at thermal equilibrium with a bath at temperature . If we assume the validity of the BornOppenheimer approximation and consider electrons in their ground state, the electrons provide a manybody potential for the ionic motion. The ionic “Hamiltonian” is
(1) 
where and are the set of momentum and coordinate operators of the ions, respectively, and is the electronic energy which depends parametrically on the nuclear positions. The thermodynamics of the system is derived from the knowledge of the Helmholtz free energy
(2) 
where and the trace is over nuclear degrees of freedom. The term can be computed at various levels of approximation, for instance using empirical effective potentials or Density Function Theory (DFT) or Quantum Monte Carlo (QMC) methods. Within the latter class of methods we can distinguish between Variational Monte Carlo (VMC) and Reptation Quantum Monte Carlo (RQMC) which provide different level of accuracy.
In VMC we make an ansatz for the manybody wave function and compute the electronic energy as a statistical average of the local energy over the distribution
(3) 
where is the vector of the 3Nelectronic positions and indicates the electronic hamiltonian
(4) 
In an electronion system, although not indicated explicitly, the trial wave function depends parametrically on the ionic coordinates and generally on parameters that need to be optimized by the use of the variational theorem.
In RQMC the trial wave function is automatically optimized by the application of the projection operator which filters out its excited states components. All excited states will be suppressed exponentially fast with increasing , the rate of convergence to the ground state depending on the energy gap between the ground and the first excited state nonorthogonal to the trial function. Here t is a parameter with the units of inverse energy; we will refer to it as the projection time because of the analogy between the projection operator and the propagator of the real time dynamics of the system, . The total energy function of time is defined as
(5) 
where is the exact ground state ^{1}^{1}1In case of fermions one must employ the fixed node approximation to avoid the sign problem and is the fixed node ground state energy which is an upper bound of the true ground state energy [2].
Similar to a thermal partition function, let us define the generating function of the moments of at time as
(6) 
The total energy at time is simply the derivative of the logarithm of
(7) 
and the variance of the energy is the second derivative
(8) 
which is nonnegative by definition and implies that the energy decreases monotonically with time. The (fixed node) ground state is reached at large time (much larger than the inverse gap) and
(9)  
(10) 
The last relation is the generalization of the zero variance principle to RQMC. Note that the average in eq.(8) is over the electronic ground state, not to be confused with the trace over ionic degrees of freedom in eq. (2).
The VMC level of description is obtained at . This observation allows us to derive an expression to connect the free energy of the system with VMC electronic energy to the free energy of the system with RQMC electronic energy. Let us call the nuclear hamiltonian of eq. (1) with electronic energy obtained after a time projection, and the corresponding free energy. The time derivative of the free energy is
(11) 
where indicates the trace over ionic degrees of freedom with statistical weight . The free energy at any positive time t is obtained from the knowledge of the VMC free energy as
(12) 
where the second term on the r.h.s. is the integral of the average, over ionic sampling, of the variance of the electronic energy. Since the variance is positive by definition, the time derivative of the free energy is negative, implying that the free energy is a monotonically decreasing function of the projection time : , that is the free energy of the VMC system is an upper bound to the fixed node free energy which is obtained at infinite projection time.
Eq. (12) is the main result of our work. It is particularly useful to increase the level of accuracy of transition lines from VMC to RQMC. Indeed once a transition point has been detected with the VMC energies, RQMC simulations at fixed thermodynamic points inside the two competing phases can be run for increasing projection time and the corrections to the free energy of both phases can be computed as we will show in the next section. This correction will determine how the transition point moves in going from VMC to RQMC accuracy. This method is much faster than performing thermodynamic integration directly within RQMC.
Note that, being a generalization of the CCI method, eq. (12) holds only if the system does not cross a phase boundary during the time projection.
3 Results for high pressure hydrogen
To test the use of eq. (12), we consider the first order liquidliquid phase transition in high pressure hydrogen. We have extended previous CEIMC calculations [11] to the isotherm, in the range of pressures corresponding to a range of specific volumes and of where we expect the transition to be located. We performed CEIMC simulations in the NVT ensemble, on systems of classical protons and electrons. We adopted backflowSlaterJastrow electronic trial wave functions [20, 21], with analytical RPA and numerically optimized twobody correlations and backflow terms. The single particle orbitals of the Slater determinants were obtained by selfconsistent DFTLDA calculations. To reduce the finite size effects we adopted TwistAveraged Boundary Conditions [22] with a fixed grid of twist angles. Moreover, finite size corrections were added to energy and pressure [23, 24]. More details can be found in ref. [11].
On the left panel of Fig. 1 we show the pressure versus the specific volume along the isotherm, as obtained by CEIMC with VMC energies (blue squares) and RQMC energies (pink circles). In the RQMC calculations the total projection time () and the time step () were chosen to provide converged energies, while pressures were corrected for finite and nonzero . Numerical values are summarized in Table 1. The two curves have the same qualitative behavior and are very close especially at large specific volume, where pressures from the two methods agree within statistical errors. Both curves have a plateau at , signifying a weakly first order phase transition (divergence of the isothermal compressibility). A peculiarity of this transition, probably related to its weakly first order character, is the absence of metastable states which typically appear in other first order transitions, such as, for instance, the gasliquid transition, and makes difficult to obtain a quantitative location of the transition point without computing the free energies of the competing phases. It follows that the coexistence pressure and the specific volume of the coexisting phases can be qualitatively located directly from the behavior without resorting to the free energy methods.
The values of the coexistence pressure, inferred from the position of the plateau in the curves in Figure 1, are of about and for the RQMC and the VMC data respectively. The EOS for the two different phases can be fitted with polynomial functions. At pressures below the transition, the system is in its normal phase (phase I), molecular and insulating. The curves can be well represented by a second order polynomial , with , and for RQMC data and , and for VMC data. For the phase II, the conducting fluid, a first order polynomial is sufficient to describe the pressure behavior; the values of the parameters are and for RQMC data and and for VMC data. The specific volumes of the two phases at coexistence can be estimated by extrapolating those fits up to the coexistence pressure, as shown in the right panel of Figure 1. The volume change at transition is in both methods . The specific volumes of the phases at coexistence are and for RQMC data and and for VMC data.
To test the ability of the correction formula (12) to predict the RQMC coexistence pressure of the metallization transition at we proceed as follows.

We assume that the coexistence value for the pressure and the specific volumes from the qualitative visual analysis above are correct and provide our reference results for the VMC system. With this assumption and using the fits to the EOS’s above, we reconstruct the Helmholtz free energy curves for the two phases up to a constant value, by integrating the VMCEOS’s. For a given phase the free energy per particle at the specific volumes is related to the free energy per particle at coexistence volume by
(13) where the final expression has been obtained by integrating a second order polynomial fit of the equation of state. The free energies of two coexisting phases can be referred to a single reference value noticing that at coexistence
(14) The resulting curves are plotted in Figure 2. The magenta open squares represent the free energy per particle of phase II and the blue open circles the free energy per particle of phase I. In order to enhance their curvature and make the visualization of the common tangent easier,we subtracted the same linear common tangent from each curves. This common tangent (represented by the horizontal black line in Figure 2) is, by construction, the coexistence pressure corresponding to the plateau in the VMC pressure curve as discussed above.

Once the free energy curves for the VMC system are known we can apply the correction term in eq. (12) to obtain the RQMC free energy curves. The most convenient strategy is to compute the correction at a single fixed specific volume for each phase and to reconstruct the RQMCfree energy curves by integrating the RQMCEOS of the two phases. To this aim, at a single fixed density in each phase we perform several RQMC calculations for increasing projection time to numerically compute the correction term
(15) The values of the reference volumes must be taken far enough from the coexistence region to avoid crossing the phase boundary during the time integration process. To satisfy this requirement we have chosen for the metallic liquid (phase II) and for the insulating liquid (phase I). In Figure 3 we compare the protonproton radial distribution functions at the two densities as obtained by CEIMC with VMC and RQMC electronic energies. We can see that at both densities the system is in the same phase regardless of the electronic method used. At both present a sharp peak centered at , which indicates that the system is in the molecular phase; at the molecular peak has completely disappeared and the form of is typical of a monoatomic fluid.
In order to evaluate the integral in eq. (12), at each reference specific volume we performed additional CEIMCRQMC simulations with projection times and a time step . The behavior of with , corrected by finite time step error, is shown in Figure 4 and is well fitted by a stretched exponential function: and its time integral is obtained by using the relation
(16) where is the Gamma function. The values of the fitting parameters in the two phases are , , at and , , at and the RQMCVMC Helmholtz free energy difference per particle results to be for the molecular insulating phase (phase I) and for the metallic phase (phase II).

finally we use the RQMCEOS to construct the RQMC free energy curves and to determine the coexistence pressure and volumes between the phases I and II by the common tangent construction. In figure 2 the RQMC free energy curves of phase I and II are represented by the closed turquoise circles and the closed purple squares, respectively. Despite the fact that the correction from the VMC free energy is quite small^{2}^{2}2We can have an idea of the relative amount of the free energy corrections by comparing the difference with the VMC free energy difference between phase I and II at the same thermodynamic points, : the contribution of the correction is then quite small as expected, about of the VMC free energy differences. , less than per particle over the range of volumes under analysis, it has a significant effect on the coexistence: the common tangent construction applied to those curves provides a coexistence pressure of and specific volumes of the two phases at coexistence of and . Those values are in excellent agreement with the values directly inferred from the RQMC pressure curve (, and ), being the differences between the two predictions within error bars.
4 Discussion and Conclusion
Within the Coupling Constant Integration strategy and for the Coupled ElectronIon Monte Carlo method, we have derived an original relation to obtain the free energy of an abinitio system with Reptation Quantum Monte Carlo electronic energies from the knowledge of the free energy of the same system but with Variation Monte Carlo electronic energies, a considerably less demanding calculation. This relation will be useful in improving the accuracy of the CEIMC method in tracing transition lines, therefore predicting phase diagrams, without increasing too much the already large computational requirements of the method. In order to test the use of the new relation and to illustrate the procedure to correct transition points we have considered the liquidliquid transition in high pressure hydrogen, a first order transition between a lower pressure, molecular and insulating phase (I) and an higher pressure, partially dissociated and conducting phase (II). We have reported new CEIMC results along the isotherm at in the density range corresponding to across the transition, an isotherm not considered before. Not presenting metastable states, this transition can be detected by the appearance of a plateau in the pressurevolume curve and the values for coexistence pressure and for the specific volumes of the phases at coexistence can be inferred without performing free energy calculations. Nonetheless, as detailed in the paper, it can still be used to illustrate our procedure for correcting transition lines. In particular, we get excellent agreement between the RQMC critical pressure inferred from the data () and the value obtained by our procedure (). A similar agreement is obtained for the specific volumes of the two phases at coexistence.
The close agreement between VMC and RQMC based Equation of States (EOSs) and transition points results from the high quality of our trial wave function for high pressure hydrogen even across the metalinsulator transition, a notoriously difficult region for DFT based first principle methods. This agreement strongly supports the use of VMC in CEIMC investigations of the phase diagram of hydrogen at high pressure, allowing calculations on larger systems or systems with quantum protons within the Path Integral formalism [7]. Moreover the formula derived in the present work can be used to improve the CEIMCVMC transition lines to RQMC accuracy.
Acknowledgments
CP is supported by the Italian Institute of Technology (IIT) under the SEED project grant n 259 SIMBEDD  Advanced Computational Methods for Biophysics, Drug Design and Energy Research. DMC is supported by DOE grant DEFG5209NA29456. This work was performed in part under the auspices of the US DOE by LLNL under Contract DEAC5207NA27344. Financial support from the Erasmus Mundus ProgramAtosim is acknowledged. Computer resources were provided by CASPUR (Italy) within the Competitive HPC Initiative, grant number: cmp09837, and by the DEISA Consortium (www.deisa.eu), cofunded through the EU FP6 project RI031513 and the FP7 project RI222919, through the DEISA Extreme Computing Initiative (DECI 2009).
References
References
 [1] R. M. Martin, Electronic Structure: Basic Theory and Practical Methods, Cambridge University Press ( Cambridge, 2004)
 [2] M. W. C. Foulkes, L. Mitas, R. J. Needs and G. Rajagopal, Rev. Mod. Phys. , 33 (2001).
 [3] R. Vuilleumier, Lect. Notes Phys. , 223 (2006).
 [4] J. Kohanoff, Electronic Structure Calculation for Solids and Molecules, Cambridge University Press (Cambridge, 2006).
 [5] D. Marx and J. Hutter, Ab Initio Molecular Dynamics: basic theory and advanced methods, Cambridge University Press (Cambridge, 2009).
 [6] C. J. Umrigar, J. Toulouse,C. Filippi, S. Sorella, and R. G. Hennig, Phys. Rev. Lett. , 110201 (2007).
 [7] C. Pierleoni and D. M. Ceperley, Lect. Notes Phys. , 641 (2006).
 [8] C. Pierleoni, D. M. Ceperley and M. Holzmann, Phys. Rev. Lett. , 146402 (2004).
 [9] M. A. Morales, C. Pierleoni and D. M. Ceperley, Phys. Rev. E , 021202 (2010).
 [10] K. T. Delaney, C. Pierleoni, and D. M. Ceperley Phys. Rev. Lett. , 235702 (2006).
 [11] M. A. Morales, C. Pierleoni, E. Schwegler and D.M. Ceperley, PNAS , 12799 (2010).
 [12] J. G. Kirkwood, J. Chem. Phys. , 315 (1935).
 [13] D. Frenkel and B. Smit, Understanding Molecular Simulation, ed., Academic Press, London (2002).
 [14] D. Alfé, M. J. Gillan, M. Towler and R. Needs, Phys. Rev. B , 161101 (2004).
 [15] W. Lorenzen, B. Holst, and R. Redmer, Phys. Rev. B , 195107 (2010).
 [16] M. A. Morales, E. Schwegler, D.M. Ceperley, C. Pierleoni, S. Hamel and K. Caspersen PNAS , 1324 (2090).
 [17] L. Landau and J. Zeldovich, Acta Physicochim. URSS , 194 (1943).
 [18] W. Ebeling and W. Richert, Phys. Lett. A , 85 (1985).
 [19] D. Saumon and G. Chabrier, Phys. Rev. A , 2084 (1992).
 [20] M. Holzmann, D.M. Ceperley, C. Pierleoni and K. Esler, Phys. Rev. E , 046707 (2003).
 [21] C. Pierleoni, K. T. Delaney, M. A. Morales, D. M. Ceperley, and M. Holzmann, Comput. Phys. Commun. , 89 (2008).
 [22] C. Lin, F. H. Zong, and D. M. Ceperley, Phys. Rev. E , 016702 (2001)
 [23] S. Chiesa, D. M. Ceperley, R. M. Martin, and M. Holzmann, Phys. Rev. Lett. , 076404 (2006)
 [24] N. D. Drummond, R. J. Needs, A. Sorouri and W. M. C. Foulkes, Phys. Rev. B , 125106 (2008)
1.310  9.417  331.6(6)  327.2(6) 
1.320  9.634  311.4(4)  305.2(6) 
1.330  9.855  293.9(6)  284.0(6) 
1.340  10.08  283.8(7)  275.9(6) 
1.350  10.31  277.2(7)  275.0(6) 
1.355  10.42  277.8(8)   
1.360  10.54  271.5(6)  262.8(6) 
1.370  10.77  253.1(6)  252.6(6) 
1.380  11.01  240.7(6)  239.2(5) 
1.390  11.25  231.0(6)  229.6(6) 
1.400  11.49  220.4(6)  217.8(6) 