# Stable liquid Hydrogen at high pressure by a novel ab-initio molecular dynamics

## Abstract

We introduce an efficient scheme for the molecular dynamics of electronic systems by means of quantum Monte Carlo. The evaluation of the (Born-Oppenheimer) forces acting on the ionic positions is achieved by two main ingredients: i) the forces are computed with finite and small variance, which allows the simulation of a large number of atoms, ii) the statistical noise corresponding to the forces is used to drive the dynamics at finite temperature by means of an appropriate Langevin dynamics. A first application to the high-density phase of Hydrogen is given, supporting the stability of the liquid phase at and .

###### pacs:

47.11.Mn, 02.70.Ss, 61.20.Ja, 62.50.+pThe phase diagram of Hydrogen at high pressure is still under intense study from the experimental and theoretical point of view. In particular in the low temperature high-pressure regime there is yet no clear evidence of a metallic atomic solid, and either a molecular solid phase can be favoured even in this regimeneeds (), or the liquid phase can be stabilized at low temperature by increasing the pressure (see Ref.(galli ())).

Indeed, for high pressures around , a two fluid (proton and electron) superconducting phase, induced by the strong electron-phonon coupling, has been conjecturedash1 () and unusual quantum properties have been later predictedashcroft (). In this work we use an improved ab-initio molecular dynamics (AMD) by using accurate forces computed by Quantum Monte Carlo (QMC). We present preliminary results, showing that the liquid phase is energetically stable, due to the strong electron correlation, at least within the Resonating Valence Bond (RVB) variational approachvanilla (), which is very accurate also in the solid phase.

AMD is well established as a powerful tool to investigate many-body condensed matter systems. Indeed, previous attempts to apply Quantum Monte Carlo (QMC) for the dynamics of ionsmitas () or for their thermodynamic propertiespenalty () are known, but they were limited to small number of electrons or to total energy corrections of the AMD trajectories, namely without the explicit calculation of the forces. Indeed, the technical achievements that we are going to present in this letter are particularly important for the simulation of liquid or disordered phases by QMC.

Calculation of forces with finite variance. The simplest method for accurate calculations within QMC, is given by the so called variational Monte Carlo (VMC), which allows to compute the variational energy expectation value of a highly accurate correlated wave function (WF) by means of a statistical approach: electronic configurations , with given electron positions and spins for , are usually generated by the Metropolis algorithm according to the probability density . Then is computed by averaging statistically over the so called local energy , namely , where indicates conventionally the multidimensional integral over the electronic coordinates weighted by . In the present work we assume that the WF is given by a correlated Jastrow factor times a determinant of a matrix , such as for instance a Slater determinant. The main ideas of this approach can be straightforwardly generalized to more complicated and more accurate WF’s, as well as to projection QMC methodsreptation () more accurate than VMC.

The efficient calculation of the energy derivatives, namely the forces , for , where is the number of atoms, is the most important ingredient for the AMD. Within VMC they can be computed by simple differentiation of , using that not only the Hamiltonian but also depend explicitly on the atomic positions . This leads to two different contributions to the force , the Hellmann-Feynman and the Pulay one , where:

(1) | |||||

(2) |

However in order to obtain a statistically meaningful average, namely with finite variance, some manipulations are necessary because the first integrand may diverge when the minimum electron-atom distance vanishes, whereas the second integrand is analogously unbounded when an electronic configuration approaches the nodal surface determined by . By defining with () the distance of from the nodal region (the minimum electron-atom distance), ( ), whereas (), leading to an unbounded integral of the square integrand in Eq.(2) (Eq.1), namely to infinite variance. The infinite variance problem in Eq.(1) was solved in several ways. Here we adopt a very elegant and efficient scheme proposed by Caffarel and Assarafcaffarel (). Instead the infinite variance problem in Eq.(2) was not considered so far, and this is clearly a problem for a meaningful definition of ionic AMD consistent with QMC forces.

In this letter we solve this problem in the following simple way, by using the so called re-weighting method. We use a different probability distribution , determined by a guiding function :

(3) |

where for is a ”measure” of the distance from the nodal surface . By assumption may vanish only when () and therefore is chosen to depend only on . For reasons that will become clear later on we have adopted the following expression:

(4) |

Then the guiding function is defined by properly regularizing , namely:

(5) |

The non obvious regularization for instead of e.g. was considered in order to satisfy the continuity of the first derivative of when , thus ensuring that remains as close as possible to the trial function . In this way the Metropolis algorithm can be applied for generating configurations according to a slightly different probability and the exact expression of can be obtained by the so called umbrella average:

(6) |

Now, the re-weighting factor , cancels out the divergence of the integrand, that was instead present in Eq.(2). Hence the mentioned integrands in the numerator and ( ) in the denominator of Eq.(6) represent bounded random variables and have obviously finite variance. In this way the problem of infinite variance is definitely solved within this simple re-weighting scheme.

We show in Fig.(1) the efficiency of the method for computing the Pulay force component acting on a Hydrogen proton at in a bcc lattice. As it is clear in the plot for the case, the difference between a method with finite variance and the standard one with infinite variance is evident. In this way the correlation matrix , defining the statistical correlation between the force components, can be efficiently evaluated:

(7) |

where the brackets indicate the statistical average over the QMC samples. The correlation matrix , that within the conventional method is not even defined, will be a fundamental ingredient for a consistent AMD with QMC forces and therefore the solution of the infinite variance problem is crucial for this purpose.

Langevin dynamics. In the following derivation we assume that ions have unit mass, that can be generally obtained by e.g. a simple rescaling of lengths for each ion independently. For clarity and compactness of notations, we also omit the ionic subindices when not explicitly necessary. Moreover matrices (vectors) are indicated by a bar (arrow) over the corresponding symbols, and the matrix-vector product is also implicitly understood. We start therefore by the following AMD equations for the ion coordinates and velocities :

(8) | |||||

(9) |

By using the fluctuation-dissipation theorem the friction matrix is related to the temperature (henceforth the Boltzmann constant ) by:

(10) |

where is generally a symmetric correlation matrix:

(11) |

It is important to emphasize that, as a remarkable generalization of the standard AMD used in parrforce (), in the present approach the friction matrix , may depend explicitly on the ion positions , so that Eq.(10) can be satisfied even for a generic correlation matrix . In fact we have the freedom to consider a QMC contribution in :

(12) |

where and is the identity matrix up to another positive constant , , and can be estimated by Eq.(7). In the following we will show that, for appropriate , it is possible to follow the Langevin dynamics by means of noisy QMC forces.

Integration of the Langevin dynamics. Henceforth the velocities are computed at half-integer times , whereas coordinates are assumed to be defined at integer times . Then, in the interval and for small, the positions are changing a little and, within a good approximation, the dependence in the Eq.(8) can be neglected, so that this differential equation becomes linear and can be solved explicitly. The closed solution can be recasted in the following useful form, where the force components appear corrected by appropriate noisy vectors :

(13) | |||||

(14) | |||||

(15) | |||||

(16) |

By using that from Eq.(10), and that its dependence on can be consistently neglected in this small time interval, the correlator defining the discrete (time integrated) noise can be computed explicitly:

(17) |

This means that the QMC noise has to be corrected in a non trivial way as explained in the following.

Noise correction. The QMC noise is given during the simulation, and therefore in order to follow the correct dynamics another noise has to be added to the noisy force components in a way that the total integrated noise is the correct expression (17), i.e. By using that the QMC noise in Eq.(7) is obviously independent of the external noise, we easily obtain the corresponding correlation matrix:

(18) |

On the other hand, after substituting the expression (12) in Eq.(10) and using the expression (17) for , we obtain a positive definite matrix in Eq.(18) for proof (). Hence is a generic Gaussian correlated noise that can be easily sampled by standard algorithms. After that the random vector is added to the force obtained by QMC, and replaces in Eq.(13). This finally allows to obtain an accurate AMD with a corresponding small time step error. The main advantage of this technique is that, at each iteration, by means of Eq.(10), the statistical noise on the total energy and forces (see Fig.2) can be much larger than the target temperature , and this allows to improve dramatically the QMC efficiency. Optimization of the WF. In the following examples we consider a cubic box with periodic boundary conditions, and use a variational WF that is able to provide a very accurate description of the correlation energy, due to a particularly efficient choice of the determinant factor, that allows to describe the RVB correlations casulamol (); benzenedim (). The WF contains several variational parameters, indicated by a vector , that have to be consistently optimized during the AMD. The Jastrow factor used here depends both on the charge and spin densities, and it is expanded in a localized atomic basis. As it is shown in the table, the accuracy of our WF is remarkable. Indeed the small difference between the so called DMC -providing the lowest possible variational energy within the same nodal surface of - and the VMC energies clearly supports the accuracy of our calculation. The size effects are very large for the metal and have been estimated () following Ref.oldh, .

N | pier () | pier () | ||
---|---|---|---|---|

16 | -0.48875(5) | -0.4878(1) | -0.49164(4) | -0.4905(1) |

54 | -0.53573(2) | -0.5353(2) | -0.53805(4) | -0.5390(5) |

128 | -0.49495(1) | -0.4947(2) | -0.49661(3) | -0.4978(4) |

250 | -0.49740(2) | - | -0.49923(2) | - |

432 | -0.49943(3) | - | - | - |

-0.501(1) | - | -0.503(1) | - |

In order to optimize the WF we use the recent method introduced in Ref.hesscyrus, , devised here in an appropriate way to optimize a large number of parameters during the AMD simulation, as described in Ref.benzenedim, . This allows to remain efficiently within the Born-Oppenheimer energy surface each time the ionic positions are changed according to Eq.(13).

Application to high-pressure Hydrogen. We show in Fig.(2) the evolution of the internal energy and corresponding temperature as a function of time with the proposed AMD with QMC forces, starting from the bcc Hydrogen solid at , considered henceforth. Although we have not studied the possible stability of all other solid phases yet, for also the simple hexagonal structure (SH) melts. This already provides a clear support to the liquid phase because these two atomic solids are the most stable ones, so far proposed at zero temperature. The proton-proton correlation function obtained starting from the two atomic solids is shown in the inset.

We performed a finite size scaling analysis based on the comparison of our QMC results with the LDA onesoldh (). We found that LDA favours atomic solid phases respect to QMC. For instance the internal energy difference between the liquid phase and the BCC solid one is only -0.0004 Ha in LDA while in QMC is -0.011 Ha. On the other hand the molecular solid (mhpc-coldh ()) is also similarly preferred to the atomic solid by QMC, but appears to have much larger zero point energy corrections compared to the liquid.russo () Indeed we have studied quantum effects on protons by using the Wigner-Kirkwood expansion on the free energy and obtained for the liquid, namely a much smaller correction than the SH [] and the molecular solid [] ones. Finally even the more accurate DMC does not affect the liquid stability, because it lowers the internal energy of all the phases studied by about the same amount ().

In conclusion we have shown that it is possible to make a realistic and accurate AMD simulation with QMC forces. We found that the bcc and SH solid structures appear clearly unstable even at low temperatures where a molecular liquid has much lower internal energy, and is further stabilized by considering quantum effects on protons. This important finding highlights the present QMC technique as a possible and accurate alternative to study phase diagrams of materials and as a benchmark for other approximate methods.

We acknowledge partial support by PRIN MIUR and CNR. We thank D.M. Ceperley, C. Pierleoni, R. Car, F. Becca, M. Casula, M. Fabrizio for useful discussions, and the excellent stability of SP5 in CINECA.

### References

- S. A. Bonev, E. Schwegler, T. Ogitsu and G. Galli, Nature, 431, 669 (2004).
- N.W. Ashcroft , J.Phys. A 12, A129-137 (2000).
- E. Babaev, A. Sudbe, and N. W. Ashcroft, Nature, 431, 666 (2004),
- C. J. Pickard1 and R. J. Needs, Nature Physics, 3, 473 (2007).
- see e.g. P.W. Anderson et al. J. Phys. Cond. Mat. 16 R755-R769 (2004) and references therein.
- J. C. Grossman and L. Mitas Phys. Rev. Lett. 94, 056403 (2005)
- C. Pierleoni, D.M. Ceperley, and M. Holzmann Phys. Rev. Lett. bf 93, 146402 (2004), D.M. Ceperley and M. Dewing J. Chem. Phys. 110, 9812 (1999).
- S. Moroni and S. Baroni Phys. Rev. Lett. 82, 4745 (1999).
- R. Assaraf and M. Caffarel J. Chem. Phys. 113, 4028 (2000).
- F. R. Krajewski and Michele Parrinello, Phys. Rev. B73, 041105(R) (2006).
- By definition commutes with , and therefore they have common eigenvectors. On the other hand each eigenvalue of is positive, and the corresponding one for is , namely a monotonically increasing function of . Thus this eigenvalue is greater than the one evaluated for the minimum value, namely for and . This bound remains positive, vanishing only for .
- M. Casula, C. Attaccalite and S. Sorella J. Chem. Phys. 121 7110 (2004).
- S. Sorella, M. Casula and D. Rocca, J. Chem. Phys. in press.
- M. Holzmann, D. M. Ceperley, C. Pierleoni, and K. Esler, Phys. Rev. E68, 046707 (2003).
- C. Umrigar et al., Phys. Rev. Lett. 98, 110201 (2007).
- V. Natoli, R. M. Martin, and D. M. Ceperley Phys. Rev. Lett. 70 1952 (1993), ibidem 74 1601 (1995).
- V. V. Kechin, JETP Letters 79, 40 (2004).