Efficient single-run implementation of generalized Einstein relation to compute transport coefficients: A binary-based time sampling.

A robust and simple implementation of the generalized Einstein formulation using single equilibrium molecular dynamics simulation is introduced to compute diffusion and shear viscosity. The unique features underlying this framework are as follows: (1) The use of a simple binary-based method to sample time-dependent transport coefficients results in a uniform distribution of data on a logarithmic time scale. Although we sample "on-the-fly," the algorithm is readily applicable for post-processing analysis. Overlapping same-length segments are not sampled as they indicate strong correlations. (2) Transport coefficients are estimated using a power law fitting function, a generalization of the standard linear relation, that accurately describes the long-time plateau. (3) The use of a generalized least squares (GLS) fitting estimator to explicitly consider correlations between fitted data points results in a reliable estimate of the statistical uncertainties in a single run. (4) The covariance matrix for the GLS method is estimated analytically using the Wiener process statistics and computed variances. (5) We provide a Python script to perform the fits and automate the procedure to determine the optimal fitting domain. The framework is applied to two fluids, binary hard sphere and a Lennard-Jones near the triple point, and the validity of the single-run estimates is verified against multiple independent runs. The approach should be applicable to other transport coefficients since the diffusive limit is universal to all of them. Given its rigor and simplicity, this methodology can be readily incorporated into standard molecular dynamics packages using on-the-fly or post-processing analysis.

Effects of thermostatting in molecular dynamics on anharmonic properties of crystals: Application to fcc Al at high pressure and temperature.

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).

Harmonic Oscillator Staging Coordinates for Efficient Path Integral Simulations of Quantum Oscillators and Crystals.

Imaginary-time path integral (PI) is a rigorous quantum mechanical tool to compute static properties at finite temperatures. However, the stiff nature of the internal PI modes poses a sampling challenge. This is commonly tackled using staging coordinates, in which the free particle (FP) contribution of the PI action is diagonalized. We introduce novel and simple staging coordinates that diagonalize the entire action of the harmonic oscillator (HO) model, rendering it efficiently applicable to (exclusively) systems with harmonic character, such as quantum oscillators and crystals. The method is not applicable to fluids or systems with imaginary modes. Unlike FP staging, the HO staging provides a unique treatment of the centroid mode. We provide implementation schemes for PIMC and PIMD simulations in NVT ensemble. Sampling efficiency is assessed in terms of the precision and accuracy of estimating the energy and heat capacity of a one-dimensional HO and an asymmetric anharmonic oscillator (AO). In PIMC, the HO coordinates propose collective moves that perfectly sample the HO contribution, then (for AO) the residual anharmonic term is sampled using standard Metropolis method. This results in a high acceptance rate and, hence, high precision, in comparison to the FP staging. In PIMD, the HO coordinates naturally prescribe definitions for the fictitious masses, yielding equal frequencies of all modes when applied to the HO model. This allows for a substantially larger time step sizes relative to standard staging, without affecting accuracy or integrator stability. For completeness, we also present results using normal mode (NM) coordinates, based on both HO and FP models. While staging and NM coordinates show similar performance (for FP or HO), staging is computationally preferable due to its cheaper scaling with the number of beads. The simplicity and the enhanced sampling gained by the HO coordinates open avenues for efficient estimation of nuclear quantum effects in more complex systems with harmonic character, such as real molecular bonds and quantum crystals.

A comparative study of methods to compute the free energy of an ordered assembly by molecular simulation.

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.

No system-size anomalies in entropy of bcc iron at Earth's inner-core conditions.

New molecular modeling data show that the entropy of bcc iron exhibits no system-size anomalies, implying that it should be feasible to compute accurate free energies of this system using first-principles methods without requiring a prohibitively large number of atoms. Conclusions are based on rigorous calculations of size-dependent free energies for a Sutton-Chen model of iron previously fit to ab initio calculations, and refute statements recently appearing in the literature indicating that the size of the simulation cell is critical for stabilization of the bcc phase.

Harmonically Assisted Methods for Computing the Free Energy of Classical Crystals by Molecular Simulation: A Comparative Study.

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.

Reformulation of Ensemble Averages via Coordinate Mapping.

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.

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.