Acceleratingab initioMolecular Dynamics and Probing the Weak Dispersive Forces in Dense Liquid Hydrogen

Abstract

We propose an ab-initio molecular dynamics method, capable to reduce dramatically the autocorrelation time required for the simulation of classical and quantum particles at finite temperature. The method is based on an efficient implementation of a first order Langevin dynamics modified by means of a suitable, position dependent acceleration matrix S. Here we apply this technique, within a Quantum Monte Carlo (QMC) based wavefunction approach and within the Born-Oppheneimer approximation, for determining the phase diagram of high-pressure Hydrogen with simulations much longer than the autocorrelation time. With the proposed method, we are able to equilibrate in few hundreds steps even close to the liquid-liquid phase transition (LLT). Within our approach we find that the LLT transition is consistent with recent density functionals predicting a much larger transition pressures when the long range dispersive forces are taken into account.One of the most important problems in ab-initio molecular dynamics (MD) is the coexistence of much different time scales underlying physical processes, even when computing equilibrium finite temperature properties. This has been a challenge for most ab-initio simulations in systems of biophysical interest, the most striking one being protein folding [1,2], where the microscopic time scale of molecular vibration is of the order of f s, while the macroscopic time scale of the folding exceeds the µs. Also in water system simulations, the weak Van der Waals (vdW) interactions and the related hydrogen bond imply very difficult equilibration properties at ambient conditions, so that the most accurate simulations are based on force field fitting forces. We will show here that, a computationally expensive ab-initio method for dense liquid hydrogen can be combined with a sampling technique capable to reduce drastically the autocorrelation times. This method has the potential, due to its simplicity, to be applied with success also to much more complex systems and other ab-initio techniques.At variance of all previous attempts[3-10], we propose that an optimal way to get rid of different time scales is based on the use of first order Langevin dynamics:where T is the temperature, f Ri = −∂ Ri V ( R) and V R is an energy potential of the classical coordinates R, e.g. atomic positions. For S ij = δ ij this is the conventional first order Langevin dynamics (CFOLD). Preconditioned Langevin equations, formally similar to Eq. 40, have been already introduced in several fields different from molecular dynamics, such as gauge field theories[11], statistics [12] and machine learning [13], with different definitions of the preconditioning matrix S. Instead, the main idea of this work is based on the simple solution of CFOLD for the harmonic potential V (R) = −K R, where K is the elastic matrix term. In this case one is able to show that the i) the autocorrelation time is independent of temperature τ −1 corr = K min , where K min (K max ) is the lowest (largest) non-zero eigenvalue of the elastic matrix...

Phase Diagram of Hydrogen and a Hydrogen-Helium Mixture at Planetary Conditions by Quantum Monte Carlo Simulations

“…However, the few QMC results for the H phase diagram, until now, have not agreed well. In particular, while all QMC simulations agreed qualitatively on a larger dissociation and metallization pressure for pure dense liquid H, compared to PBE, the precise location was not welldetermined due to different QMC implementations, variational wave-function, and finite size effects errors 9,[20][21][22][23] . We perform simulations with 64 and 128 H atoms for the H compound and with 118 H and 10 He atoms for the H-He mixture (see Methods).…”

“…We plot g(r) at different particle densities (defined using the Wigner-Seitz radius defined above) at each temperature in order to constrain the metallization/dissociation crossover. The single Z basis refers to the one adopted in the previous work 23 , where the accelerated molecular dynamics was introduced. 2Z and 3Z basis sets refer instead to the more conventional "correlation consistent polarized valence x−zeta" (ccpVxZ , with x = 2, 3, respectively) sequence 40 .…”

“…The plane wave basis set depends only on the kinetic energy cutoff, taken large enough for converged results, and is obtained by a standard electronic structure package for 64 atoms. ∆0 here is a parameter that is used to optimize the time step error ∆t of the molecular dynamics 23 . For small enough ∆t all results should be unbiased.…”

“…We have used the recently introduced accelerated molecular dynamics 23 to sample the classical canonical distribution in the NVT ensemble. In this dynamics a crucial role is played by the parameter ∆ 0 that is used to scale the covariance matrix S used for the acceleration in order to represent as close as possible the Hessian matrix.…”

“…allization pressure: the electronic size effects errors, not adequately removed in Mazzola et al 20,22 and a localized basis set 29 that was too small (1Z) to describe the metal and the insulator with the same accuracy 23 . In addition, our previous studies used a less efficient optimization method, indeed, the so called "linear method" 30 requires a careful generalization to the case of complex wavefunctions (see Methods for details).…”