Path-integral Monte Carlo simulations of solid parahydrogen using two-body, three-body, and four-body ab initio interaction potential energy surfaces

Abstract

We present path integral Monte Carlo simulation results for the equation of state of solid parahydrogen between $ 0.024 \, {\r{A}}^{-3} $ and $ 0.1 \, {\r{A}}^{-3} $ at $ T = 4.2 \, $ K. The simulations are performed using non-additive isotropic ab initio two-body, three-body, and four-body potential energy surfaces (PES). We apply corrections to account for both the finite size simulation errors and the Trotter factorization errors. Simulations that use only the two-body PES during sampling yield an equation of state similar to that of simulations that use both the two-body and three-body PESs during sampling. With the four-body interaction energy, we predict an equilibrium density of $ 0.02608 \, {\r{A}}^{-3} $, very close to the experimental result of $ 0.0261 \, {\r{A}}^{-3} $. The inclusion of the four-body interaction energy also brings the simulation results in excellent agreement with the experimental pressure-density data until around $ 0.065 \, {\r{A}}^{-3} $, beyond which the simulation results overestimate the pressure. These PESs overestimate the average kinetic energy per molecule at the equilibrium density by about $ 7 \% $ compared to the experimental result. Our findings suggest that, at higher densities, we require five-body and higher-order many-body interactions to quantitatively improve the agreement between the pressure-density curve produced by simulations, and that of experiment. Using the four-body PES during sampling at excessively high densities, where such higher-order many-body interactions are likely to be significant, causes an artificial symmetry breaking in the hcp lattice structure of the solid.