Accelerating ab initio Molecular Dynamics and Probing the Weak Dispersive Forces in Dense Liquid Hydrogen

We propose an ab initio molecular dynamics method, capable of dramatically reducing the autocorrelation time required for the simulation of classical and quantum particles at finite temperatures. 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 to both Lennard-Jones models, to demonstrate the accuracy and speeding-up of the sampling, and within a quantum Monte Carlo based wave function approach, 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 a few hundred 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 pressure when the long range dispersive forces are taken into account.

Figure 1

Structural optimization of linear structures of length $L_c$. We use both the standard steepest descent method (black) and the accelerated one [red, Eq. (3)] with different choices for the preconditioning matrix, depending on the system. The approach to the stable minimum energy configuration is shown as a function of the number of iterations at fixed optimal time step $\mathrm{\Delta }$ for both cases. (a) System of five harmonic dimers (having an equilibrium distance of 1.4 a.u. and spring stiffness $k=10$, see SM [14]), interacting via a Lennard-Jones potential $V=4{\epsilon }_0[({\sigma }_0/r{)}^{12}-({\sigma }_0/r{)}^6]$, with ${\epsilon }_0=0.03$ and ${\sigma }_0=2$. The starting configuration is given by equally spaced particles at distance $r=1.2$ a.u. We use the standard Newton method with a regularized Hessian matrix in the accelerated case (red). (b) Here, we have ten hydrogen atoms, whose interactions are given by an ab initio QMC calculation. The starting configuration has equally spaced atoms at distance 1 a.u. The preconditioning matrix is given by Eq. (7) and $K_{\mathrm{cond}}\simeq 200$ in this case.

Figure 2

Average potential energy [in atomic units (a.u.)] as a function of the integration time step $\mathrm{\Delta }$ in two LJ models. Black points refer to the CFOLD, while red points to the improved one [Eq. (4)], with $S$ given by the regularized Hessian matrix (see SM [14] for details). We numerically prove that the preconditioned MD does not introduce any bias in the sampled distribution [see, also, the converged radial pair distribution functions $g(r)$ in the insets (b), (d)] and drastically improves the time-step error. Continuous lines indicate linear fit of the data. (a) System of harmonic dimers interacting via a LJ potential. Here, $K_{\mathrm{cond}}\approx 25$. For a fair comparison, we tune the regularization parameters of the Hessian matrix in order to have the same particle diffusion—at fixed $\mathrm{\Delta }$—for the CFOLD and the preconditioned case. Notice that the drift-diffusion term [second line of Eq. (4)] is essential to obtain the correct distribution [blue line in (b)]. (c) Low temperature atomic solid of particles interacting via a pure LJ potential. Notice the logarithmic scale in the $x$ axis. Here, we normalize $\mathrm{\Delta }$ using the autocorrelation time $\tau$. $K_{\mathrm{cond}}\approx 100$ in this case.

Figure 3

Phase diagram of liquid hydrogen. We indicate experimental results of Refs. [41, 43] with solid symbols, while theoretical results (with previous VMC and Coupled Electron Ion Monte Carlo (CEIMC) simulations, with classical nuclei approximation) of Refs. [47, 52] with open symbols. The present results are plotted in red. We also plot DFT metallization lines using three popular exchange correlation functionals, PBE, vdW-DF1 and DF2. These results are taken from Ref. [43].