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

Investigating the role of dispersion corrections and anharmonic effects on the phase transition in SrZrS3: A systematic analysis from AIMD free energy calculations

A thermally driven needle-like (NL) to distorted perovskite (DP) phase transition in SrZrS3 was investigated by means of ab initio free energy calculations accelerated by machine learning. As a first step, a systematic screening of the methods to include long-range interactions in semilocal density functional theory Perdew-Burke-Ernzerhof calculations was performed. Out of the ten correction schemes tested, the Tkatchenko-Scheffler method with iterative Hirshfeld partitioning method was found to yield the best match between calculated and experimental lattice geometries, while predicting the correct order of stability of NL and DP phases at zero temperature. This method was then used in free energy calculations, performed using several approaches, so as to determine the effect of various anharmonicity contributions, such as the anisotropic thermal lattice expansion or the thermally induced internal structure changes, on the phase transition temperature (TNP→DP). Accounting for the full anharmonicity by combining the NPT molecular dynamics data with thermodynamic integration with harmonic reference provided our best estimate of TNL→DP = 867 K. Although this result is ∼150 K lower than the experimental value, it still provides an improvement by nearly 300 K compared to the previous theoretical report by Koocher et al. [Inorg. Chem. 62, 11134-11141 (2023)].

Anharmonic Correction to Adsorption Free Energy from DFT-Based MD Using Thermodynamic Integration

Adsorption processes are often governed by weak interactions for which the estimation of entropy contributions by means of the harmonic approximation is prone to be inaccurate. Thermodynamic integration (TI) from the harmonic to the fully interacting system (λ-path integration) can be used to compute anharmonic corrections. Here, we combine TI with (curvilinear) internal coordinates in periodic systems to make the formalism available in computational studies. Our implementation of ab initio molecular dynamics in VASP is independent of the reaction path and can be thus applied to study adsorption processes relative to the gas phase and does hence provide a useful tool for computational catalysis. We discuss the application of the approach on three model systems for which exact semianalytical solutions exist and illustrate and quantify the importance of anharmonic vibrations, hindered rotations, and hindered translations (dissociation). Eventually, we apply the method to study the adsorption of small adsorbates in a zeolite (H-SSZ-13).

Implementation of harmonically mapped averaging in LAMMPS, and effect of potential truncation on anharmonic properties

Implementation of the harmonically mapped averaging (HMA) framework in the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) is presented for on-the-fly computations of the energy, pressure, and heat capacity of crystalline systems during canonical molecular dynamics simulations. HMA has a low central processing unit and storage requirements and is straightforward to use. As a case study, the properties of the Lennard-Jones and embedded-atom model (parameterized for nickel) crystals are computed. The results demonstrate the higher efficiency of the new class compared to the inbuilt LAMMPS classes for calculating these properties. However, HMA loses its effectiveness in systems where diffusion occurs in the crystal, and an example is presented to allow this behavior to be recognized. In addition to its improved precision, HMA is less affected by small errors introduced by having a larger time step in molecular dynamics simulations. We also present an analysis of the effect of potential truncation on anharmonic properties, and show that artifacts of truncation on the HMA averages can be eliminated simply by shifting the potential energy to zero at the truncation radius. Full properties can be obtained by adding easily computed values for the lattice and harmonic properties using the untruncated potential.

Free Energy of Classical Molecular Crystals by Thermodynamic Integration from a Harmonic Reference

We develop an algorithm for calculating the normal modes of vibration of mechanical systems with constraints, particularly of molecules with rigid bonds and models of rigid molecules, and use it to obtain the harmonic free energy of a crystal. The anharmonic correction is then calculated by the conventional thermodynamic integration over temperature in the NVT ensemble. Attention is paid to finite-size errors, tail corrections, thermostat choice, ergodicity, and other sources of inaccuracies. The calculated free energy of ice XIV modeled by the TIP4P/2005 potential agrees with the previously reported value and is by one order more accurate.

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

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.

Powder diffraction and crystal structure prediction identify four new coumarin polymorphs

Coumarin, a simple, commodity chemical isolated from beans in 1820, has, to date, only yielded one solid state structure. Here, we report a rich polymorphism of coumarin grown from the melt. Four new metastable forms were identified and their crystal structures were solved using a combination of computational crystal structure prediction algorithms and X-ray powder diffraction. With five crystal structures, coumarin has become one of the few rigid molecules showing extensive polymorphism at ambient conditions. We demonstrate the crucial role of advanced electronic structure calculations including many-body dispersion effects for accurate ranking of the stability of coumarin polymorphs and the need to account for anharmonic vibrational contributions to their free energy. As such, coumarin is a model system for studying weak intermolecular interactions, crystallization mechanisms, and kinetic effects.