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

We present path integral Monte Carlo simulation results for the equation of state of solid parahydrogen between 0.024 and 0.1Å-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 (PESs). 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Å-3, very close to the experimental result of 0.0261Å-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Å-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 the 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.