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

Anharmonic Correction to Free Energy Barriers from DFT-Based Molecular Dynamics Using Constrained Thermodynamic Integration

For the calculation of anharmonic contributions to free energy barriers, constrained thermodynamic λ-path integration (λ-TI) from a harmonic reference force field to density functional theory is presented as an alternative to the established Blue Moon ensemble method (ξ-TI), in which free energy gradients along the reaction coordinate ξ are integrated. With good agreement in all cases, the λ-TI method is benchmarked against the ξ-TI method for several reactions, including the internal CH₃ group rotation in ethane, a nucleophilic substitution of CH₃Cl, a retro-Diels-Alder reaction, and a proton transfer in zeolite H-SSZ-13. An advantage of λ-TI is that one can use virtually any reference state to compute anharmonic contributions to reaction free energies or free energy barriers. This is particularly relevant for catalysis, where it is now possible to compute anharmonic corrections to the free energy of a transition state relative to any reference, for example, the most stable state of the active site and the reactants in the gas phase. This is in contrast to ξ-TI, where free energy barriers can only be computed relative to an initial state with all reactants coadsorbed. Finally, the Bennett acceptance ratio method combined with λ-TI is demonstrated to reduce the number of required integration grid points with tolerable accuracy, favoring thus λ-TI over ξ-TI in terms of computational efficiency.

In Silico Drug Repurposing Approach: Investigation of Mycobacterium tuberculosis FadD32 Targeted by FDA-Approved Drugs

Background: Despite the enormous efforts made towards combating tuberculosis (TB), the disease remains a major global threat. Hence, new drugs with novel mechanisms against TB are urgently needed. Fatty acid degradation protein D32 (FadD32) has been identified as a promising drug target against TB, the protein is required for the biosynthesis of mycolic acids, hence, essential for the growth and multiplication of the mycobacterium. However, the FadD32 mechanism upon the binding of FDA-approved drugs is not well established. Herein, we applied virtual screening (VS), molecular docking, and molecular dynamic (MD) simulation to identify potential FDA-approved drugs against FadD32. Methodology/Results: VS technique was found promising to identify four FDA-approved drugs (accolate, sorafenib, mefloquine, and loperamide) with higher molecular docking scores, ranging from -8.0 to -10.0 kcal/mol. Post-MD analysis showed that the accolate hit displayed the highest total binding energy of -45.13 kcal/mol. Results also showed that the accolate hit formed more interactions with FadD32 active site residues and all active site residues displayed an increase in total binding contribution. RMSD, RMSF, Rg, and DCCM analysis further supported that the presence of accolate exhibited more structural stability, lower bimolecular flexibility, and more compactness into the FadD32 protein. Conclusions: Our study revealed accolate as the best potential drug against FadD32, hence a prospective anti-TB drug in TB therapy. In addition, we believe that the approach presented in the current study will serve as a cornerstone to identifying new potential inhibitors against a wide range of biological targets.

Shadow Hamiltonian in classical NVE molecular dynamics simulations involving Coulomb interactions

Microcanonical ensemble (NVE) Molecular Dynamics (MD) computer simulations are performed with negligible energy drift for systems incorporating Coulomb interactions and complex constraint schemes. In principle, such systems can now be simulated in the NVE ensemble for millisecond time scales, with no requirement for system thermostatting. Numerical tools for assessing drift in MD simulations are outlined, and drift rates of 10^-6 K/μs are demonstrated for molten salts, polar liquids, and room temperature ionic liquids. Such drift rates are six orders of magnitude smaller than those typically quoted in the literature. To achieve this, the standard Ewald method is slightly modified so the first four derivatives of the real space terms go smoothly to zero at the truncation distance, r_c. New methods for determining standard Ewald errors and the new perturbation errors introduced by the smoothing procedure are developed and applied, these taking charge correlation effects explicitly into account. The shadow Hamiltonian, E_s, is shown to be the strictly conserved quantity in these systems, and standard errors in the mean of one part in 10¹⁰ are routinely calculated. Expressions for the shadow Hamiltonian are improved over previous work by accounting for O(h⁴) terms, where h is the MD time step. These improvements are demonstrated by means of extreme out-of-equilibrium simulations. Using the new methodology, the very low diffusion coefficients of room temperature 1-hexyl-3-methyl-imidazolium chloride are determined from long NVE trajectories in which the equations of motion are known to be integrated correctly, with negligible drift.

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

Absolute chemical potentials for complex molecules in fluid phases: A centroid reference for predicting phase equilibria

Solid-fluid phase equilibria are difficult to predict in simulations because bound degrees of freedom in the crystal phase must be converted to free translations and rotations in the fluid phase. Here, we avoid the solid-to-fluid transformation step by starting with chemical potentials for two reference systems, one for the fluid phase and one for the solid phase. For the solid, we start from the Einstein crystal and transform to the fully interacting molecular crystal. For the fluid phase, we introduce a new reference system, the "centroid," and then transform to gas phase molecules. We illustrate the new calculations by predicting the sublimation vapor pressure of succinic acid in the temperature range of 300 K-350 K.

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.