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

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.