Multithermal-Multibaric Molecular Simulations from a Variational Principle

We present a method for performing multithermal-multibaric molecular dynamics simulations that sample entire regions of the temperature-pressure (TP) phase diagram. The method uses a variational principle [Valsson and Parrinello, Phys. Rev. Lett. 113, 090601 (2014)] in order to construct a bias that leads to a uniform sampling in energy and volume. The intervals of temperature and pressure are taken as inputs and the relevant energy and volume regions are determined on the fly. In this way the method guarantees adequate statistics for the chosen TP region. We show that our multithermal-multibaric simulations can be used to calculate all static physical quantities for all temperatures and pressures in the targeted region of the TP plane. We illustrate our approach by studying the density anomaly of TIP4P/ice water.

Figure 1

Illustration of the algorithm to find the energy interval $E_1–E_2$ from the temperature range ${\beta }_1-{\beta }_2$. Free energies at temperatures $\beta$, ${\beta }_1$, and ${\beta }_2$ are shown in blue, green, and orange solid lines. $p(E)$ is shown with a gray solid line. The limits of the energy interval and the energy threshold $\epsilon$ are shown in black dashed lines.

Figure 2

Density of TIP4P/ice liquid water at 1 bar as a function of temperature. The blue circles with error bars are the results of individual isothermal-isobaric simulations. Errors have been calculated using block averages [25, 26]. The orange dashed line is the mean density calculated using Eq. (8) from a single multithermal simulation. The orange filled region is the error of the mean. The melting temperature $T_m$ and the temperature of maximum density of the model $T_{\max }$ are also indicated.

Figure 3

Idealized illustration of the mapping from a set of (inverse) temperatures and pressures to the set of relevant energies and volumes. The target distribution $p(E,\mathcal{V})$ corresponds to the latter set. A realistic version of this plot is shown in the Supplemental Material [12].

Figure 4

Density of TIP4P/ice liquid water as a function of temperature for different pressures. The solid lines are the result from a single multithermal-multibaric simulation. The circles are references from individual isothermal-isobaric simulations at 1 bar and 300 MPa.

Figure 5

Density of TIP4P/ice liquid water as a function of temperature and pressures. The temperature of maximum density $T_{\max }$ for different pressures is shown in black circles. $T_{\max }$ is also represented using a $B$-spline approximation to the points. The region were liquid water is undercooled with respect to ice Ih has been shaded in gray. We stress that only the liquid state was sampled in the simulation.

Figure 6

Oxygen-oxygen radial distribution functions $g_{OO}(r)$ of TIP4P/ice liquid water for different temperatures and pressures. Results from the multithermal-multibaric simulations are shown in solid blue, orange, and green lines for the temperatures of 260, 300, and 350 K, respectively. The columns correspond to pressures of 0.1 (1bar), 100, and 300 MPa. References calculated during individual isothermal-isobaric simulations are shown in dashed gray lines.