Very fast averaging of thermal properties of crystals by molecular simulation

Knowledge of approximate harmonic behavior of crystals is introduced into a new “mapped averaging” framework to yield alternative expressions for the thermodynamic properties of crystalline systems. The expressions separate the known harmonic behavior from residual averages, which thus encapsulate anharmonic contributions to the properties. With harmonic contributions removed, direct measurement of these anharmonic contributions by molecular simulation can be accomplished without contamination by noise produced by the already-known harmonic behavior. We show with application to the Lennard-Jones model that first-derivative properties (pressure, energy) can be obtained to a given precision via this harmonically mapped averaging at least 10 times faster than by using conventional averaging, and second-derivative properties (e.g., heat capacity) are obtained at least 100 times faster; in more favorable cases, the speedup exceeds a millionfold. Free-energy calculations are accelerated by 50 to 1000 times. Data obtained using these formulations are rigorous and not subject to any added approximation, and in fact are less sensitive to inaccuracies relating to finite-size effects, potential truncation, equilibration, and similar considerations. Moreover, the approach does not require any alteration in how sampling is performed during the simulation, so it may be used with standard Monte Carlo or molecular dynamics methods. However, the mapped averages do require evaluation of first and second derivatives of the intermolecular potential, for evaluation of first and second thermodynamic-derivative properties, respectively. Apart from its usefulness to simulation, the formalism developed here may constitute a basis for new theoretical treatments of crystals.

Figure 1

Views of the precision of harmonically mapped averaging in comparison to conventional averaging. Values plotted in (a) are for each configuration; (b) and (c) are averages over $1\times {10}^7$ MC steps. Lines in (b) and (c) simply join the points, and bars are 68% confidence intervals. Difficulty ratio in (d), (e), and (f) is ratio of $D\equiv t^{1/2}\sigma$ for conventional versus mapped averaging, where $t$ is CPU time required to obtain a result with uncertainty $\sigma$ [15]. To give a physical scale, using argonlike parameters, the pressure for states in (d) ranges from about $-150$ to $+150$ MPa, and for (e) about 1 TPa to 200 MPa. Data in (f) show performance along the melting line.

Figure 2

Comparison of the uncertainties from the conventional and mapped measurements of heat capacity using “direct” method (i.e., formulas in Table 1) vs differentiation of a polynomial fit of $U_{\mathrm{anh}}/T^2$.

Figure 3

Difficulty in measuring the free energy. The small jumps in the TI results are due to changing the order of the polynomial fit (see text).

Figure 4

Finite-size effects of measuring the pressure directly (no correction) and using the HarmFSC correction $[P_{\mathrm{harm}}\left(\infty \right)-P_{\mathrm{harm}}\left(N\right)]$. All data were computed via harmonically mapped averaging with $r_c=3.0$. Inset: Effect of temperature on the finite-size effects of the corrected pressure. Each line is the HarmFSC-corrected pressure for the indicated temperature, minus the value for $N\to \infty$ (i.e., the $y$ intercept), which is subtracted to allow a common scale for all temperatures. Lines are linear fits weighted by uncertainties.

Figure 5

Infinite-range pressure as estimated by three choices for the correction to the finite-range simulation average: the standard LRC (assuming homogeneous medium) [11]; LatLRC $\equiv P_{\mathrm{lat}}\left(\infty \right)-P_{\mathrm{lat}}\left(r_{\mathrm{c}}\right)$; and HarmLRC $\equiv P_{\mathrm{harm}}\left(\infty \right)-P_{\mathrm{harm}}\left(r_{\mathrm{c}}\right)$. Pressure estimates are plotted as a function of the potential-cutoff radius used in the simulation. Some points from LRC and LatLRC with short $r_{\mathrm{c}}$ are out of the scale (as indicated by their absence where other points are plotted at a given $r_c$). All data were obtained via harmonically mapped averaging of a system of $N=4000$ LJ atoms.

Figure 6

Normalized autocorrelation function of the conventional and mapped measurements of the energy. A trajectory of $2\times {10}^5$ cycles is used.

Figure 7

Number of independent energy samples out of $2\times {10}^5$ measured samples. The number of independent samples is estimated as $\mathrm{Var}/{\sigma }^2$ with $\mathrm{Var}$ and $\sigma$ the sample variance and uncertainty, respectively.

Figure 8

Equilibration rate of conventional and harmonically mapped averages. All data are normalized to the initial energy $U_0$ (off the average $U_{\mathrm{avg}}$). The dotted lines are the normalized uncertainties from a production run of $1\times {10}^4$ MC steps.

Figure 9

Convergence of the average energy with respect to MD time-step size. Red line simply joins the points. All points are generated using the same integration time, $t=5000$ (LJ units).