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.

A framework for computing the anharmonic free energy (FE) of metallic crystals using BornOppenheimer ab initio molecular dynamics (AIMD) simulation, with unprecedented efficiency, is introduced and demonstrated for the hcp phase of iron at extreme conditions (up to ≈ 290 GPa and 5000 K). The advances underlying this work are: (1) A recently introduced harmonically-mapped averaging temperature integration (HMA-TI) method reduces the computational cost by order(s) of magnitude compared to the conventional TI approach. The TI path starts from zero Kelvin, where it assumes the behavior is given exactly by a harmonic treatment; this feature restricts application to systems that have no imaginary phonons in this limit. (2) A Langevin thermostat with the HMA-TI method allows use of a relatively large MD time step (4 fs, which is about eight times larger than the size needed for the Andersen thermostat) without loss of accuracy. (3) AIMD sampling is accelerated by using density functional theory (DFT) with a low-level parameter set, then the measured quantities of selected configurations are robustly reweighted to a higher level of DFT. This introduces a speedup of about 20-30× compared to directly simulating the accurate system. (4a) The temperature (T ) dependence of the hcp equilibrium shape (i.e. c/a axial ratio) is determined (including anharmonicity), with uncertainty less than ±0.001. (4b) Electronic excitation is included through Mermin's finite-temperature extension of the T = 0 K DFT. A simple FE perturbation method is introduced to handle the difficulty associated with applying the TI method with a Tdependent geometry and (due to electronic excitation) potential-energy surface. (5) The FE in the thermodynamic limit is attained through extrapolation of only the (computationally inexpensive) quasiharmonic FE, because the anharmonic FE contribution has negligible finite-size effects. All methods introduced here do not affect the AIMD sampling -results are obtained through postprocessing -so established AIMD codes can be employed without modification. Analytical formulas fitted to the results for the variation of the equilibrium c/a ratio and FE components with T are provided. Notably, effects of magnetic excitations are not included, and may yet prove important to the overall FE; if so, it is plausible that such contributions can be added perturbatively to the FE values reported here. Notwithstanding these considerations, FE values are obtained with an estimated accuracy and precision of 2 meV/atom, suggesting that the capability to compute the phase diagram of iron at Earth's inner core conditions is within reach.

Four methods for calculation of the classical free energy of crystalline systems are compared with respect to their efficiency and accuracy. Two of the methods involve thermodynamic integration along an unphysical path (λ integration, λI), and two involve integration in temperature from the low-temperature harmonic limit (T integration, TI). Specifically, the methods considered are (1) Frenkel-Ladd integration from a noninteracting Einstein crystal reference (ECR-λI); (2) conventional integration in temperature (Conv-TI); (3) integration from an interacting quasi-harmonic reference (QHR-λI); and (4) temperature integration using harmonically mapped averaging to evaluate the integrand (HMA-TI). The latter two methods are "harmonically assisted", meaning that they exploit the harmonic nature of the crystal to greatly reduce fluctuations in the relevant averages. This feature allows them to deliver a result of much higher precision for a given computational effort, compared to ECR-λI and Conv-TI, and with no less accuracy. Regarding the harmonically assisted methods, HMA-TI has several advantages over QHR-λI with respect to the simplicity of the integration path (which promotes a more accurate result), ease of implementation, and usefulness of the data recorded along the integration path.

A general framework is established for reformulation of the ensemble averages commonly encountered in statistical mechanics. This "mapped-averaging" scheme allows approximate theoretical results that have been derived from statistical mechanics to be reintroduced into the underlying formalism, yielding new ensemble averages that represent exactly the error in the theory. The result represents a distinct alternative to perturbation theory for methodically employing tractable systems as a starting point for describing complex systems. Molecular simulation is shown to provide one appealing route to exploit this advance. Calculation of the reformulated averages by molecular simulation can proceed without contamination by noise produced by behavior that has already been captured by the approximate theory. Consequently, accurate and precise values of properties can be obtained while using less computational effort, in favorable cases, many orders of magnitude less. The treatment is demonstrated using three examples: (1) calculation of the heat capacity of an embedded-atom model of iron, (2) calculation of the dielectric constant of the Stockmayer model of dipolar molecules, and (3) calculation of the pressure of a Lennard-Jones fluid. It is observed that improvement in computational efficiency is related to the appropriateness of the underlying theory for the condition being simulated; the accuracy of the result is however not impacted by this. The framework opens many avenues for further development, both as a means to improve simulation methodology and as a new basis to develop theories for thermophysical properties.

We present a comparative study of methods to compute the absolute free energy of a crystalline assembly of hard particles by molecular simulation. We consider all combinations of three choices defining the methodology: (1) the reference system: Einstein crystal (EC), interacting harmonic (IH), or r(-12) soft spheres (SS); (2) the integration path: Frenkel-Ladd (FL) or penetrable ramp (PR); and (3) the free-energy method: overlap-sampling free-energy perturbation (OS) or thermodynamic integration (TI). We apply the methods to FCC hard spheres at the melting state. The study shows that, in the best cases, OS and TI are roughly equivalent in efficiency, with a slight advantage to TI. We also examine the multistate Bennett acceptance ratio method, and find that it offers no advantage for this particular application. The PR path shows advantage in general over FL, providing results of the same precision with 2-9 times less computation, depending on the choice of a common reference. The best combination for the FL path is TI+EC, which is how the FL method is usually implemented. For the PR path, the SS system (with either TI or OS) proves to be most effective; it gives equivalent precision to TI+FL+EC with about 6 times less computation (or 12 times less, if discounting the computational effort required to establish the SS reference free energy). Both the SS and IH references show great advantage in capturing finite-size effects, providing a variation in free-energy difference with system size that is about 10 times less than EC. This result further confirms previous work for soft-particle crystals, and suggests that free-energy calculations for a structured assembly be performed using a hybrid method, in which the finite-system free-energy difference is added to the extrapolated (1/N→0) absolute free energy of the reference system, to obtain a result that is nearly independent of system size.

The precision and accuracy of the anharmonic energy calculated in the canonical (NVT) ensemble using three different thermostats (viz., Andersen, Langevin, and Nosé-Hoover) along with no thermostat (i.e., microcanonical, NVE) are compared via application to aluminum crystals at ≈100 GPa for temperatures up to melting (4000 K) using ab initio molecular dynamics (AIMD) simulation. In addition to the role of the thermostat, the effect of using either conventional or the recently introduced harmonically mapped averaging (HMA) method is considered. The effect of AIMD time-step size ∆t on the ensemble averages gauges accuracy, while for a given ∆t, the stochastic uncertainty (computed using block averaging) provides the metric for precision. We identify the rate of convergence of block averages (with respect to block size) as an important issue in this context, as it imposes a minimum simulation length required to achieve reliable statistics, and it differs considerably among the methods. We observe that HMA with a Langevin thermostat in an NVT simulation shows the best performance, from the point of view of accuracy, precision, and simulation length. In addition, we introduce a novel HMA-based ensemble average for the temperature. In application to NVE simulations, the new formulation exhibits much smaller fluctuations compared to the conventional kinetic-energy approach; however, it provides only marginal improvement in uncertainty due to strong negative correlations exhibited by the conventional form (which acts to reduce its uncertainty but also slows convergence of the block averages).

No abstract

No abstract