S/G-1: an ab initio force-field blending frozen Hermite Gaussian densities and distributed multipoles. Proof of concept and first applications to metal cations

Classical Exchange Polarization: An Anisotropic Variable Polarizability Model

Polarizability, or the tendency of the electron distribution to distort under an electric field, often depends on the local chemical environment. For example, the polarizability of a chloride ion is larger in gas phase compared to a chloride ion solvated in water. This effect is due to the restriction the Pauli exclusion principle places on the allowed electron states. Because no two electrons can occupy the same state, when a highly polarizable atom comes in close contact with other atoms or molecules, the space of allowed states can dramatically decrease. This constraint suggests that an accurate molecular mechanics polarizability model should depend on the radial distance between neighboring atoms. This paper introduces a variable polarizability model within the framework of the HIPPO (Hydrogen-like Intermolecular Polarizable Potential) force field, by damping the polarizability as a function of the orbital overlap of two atoms. This effectively captures the quantum mechanical exchange polarization effects, without explicit utilization of antisymmetrized wave functions. We show that the variable polarizability model remarkably improves the two-body polarization energies and three-body energies of ion-ion and ion-water systems. Under this model, no manual tuning of atomic polarizabilities for monatomic ions is required; the gas-phase polarizability can be used because an appropriate damping function is able to correct the polarizability at short range.

Development of the Quantum-Inspired SIBFA Many-Body Polarizable Force Field: Enabling Condensed-Phase Molecular Dynamics Simulations

We present the extension of the Sum of Interactions Between Fragments Ab initio Computed (SIBFA) many-body polarizable force field to condensed-phase molecular dynamics (MD) simulations. The quantum-inspired SIBFA procedure is grounded on simplified integrals obtained from localized molecular orbital theory and achieves full separability of its intermolecular potential. It embodies long-range multipolar electrostatics (up to quadrupole) coupled to a short-range penetration correction (up to charge-quadrupole), exchange repulsion, many-body polarization, many-body charge transfer/delocalization, exchange dispersion, and dispersion (up to C10). This enables the reproduction of all energy contributions of ab initio symmetry-adapted perturbation theory (SAPT(DFT)) gas-phase reference computations. The SIBFA approach has been integrated within the Tinker-HP massively parallel MD package. To do so, all SIBFA energy gradients have been derived and the approach has been extended to enable periodic boundary conditions simulations using smooth particle mesh Ewald. This novel implementation also notably includes a computationally tractable simplification of the many-body charge transfer/delocalization contribution. As a proof of concept, we perform a first computational experiment defining a water model fitted on a limited set of SAPT(DFT) data. SIBFA is shown to enable a satisfactory reproduction of both gas-phase energetic contributions and condensed-phase properties highlighting the importance of its physically motivated functional form.

Ab initio effective one-electron potential operators: Applications for charge-transfer energy in effective fragment potentials

The concept of effective one-electron potentials (EOPs) has proven to be extremely useful in efficient description of electronic structure of chemical systems, especially extended molecular aggregates such as interacting molecules in condensed phases. Here, a general method for EOP-based elimination of electron repulsion integrals is presented, that is tuned toward the fragment-based calculation methodologies such as the second generation of the effective fragment potentials (EFP2) method. Two general types of the EOP operator matrix elements are distinguished and treated either via the distributed multipole expansion or the extended density fitting (DF) schemes developed in this work. The EOP technique is then applied to reduce the high computational costs of the effective fragment charge-transfer (CT) terms being the bottleneck of EFP2 potentials. The alternative EOP-based CT energy model is proposed, derived within the framework of intermolecular perturbation theory with Hartree-Fock noninteracting reference wavefunctions, compatible with the original EFP2 formulation. It is found that the computational cost of the EFP2 total interaction energy calculation can be reduced by up to 38 times when using the EOP-based formulation of CT energy, as compared to the original EFP2 scheme, without compromising the accuracy for a wide range of weakly interacting neutral and ionic molecular fragments. The proposed model can thus be used routinely within the EFP2 framework.

Fast Evaluation of Two-Center Integrals over Gaussian Charge Distributions and Gaussian Orbitals with General Interaction Kernels

We present efficient algorithms for computing two-center integrals and integral derivatives, with general interaction kernels K(r₁₂), over Gaussian charge distributions of general angular momenta l. While formulated in terms of traditional ab initio integration techniques, full derivations and required secondary information, as well as a reference implementation, are provided to make the content accessible to other fields. Concretely, the presented algorithms are based on an adaption of the McMurchie-Davidson Recurrence Relation (MDRR) combined with analytical properties of the solid harmonic transformation; this obviates all intermediate recurrences except the adapted MDRR itself, and allows it to be applied to fully contracted auxiliary kernel integrals. The technique is particularly well-suited for semiempirical molecular orbital methods, where it can serve as a more general and efficient replacement of Slater-Koster tables, and for first-principles quantum chemistry methods employing density fitting. But the formalism's high efficiency and ability of handling general interaction kernels K(r₁₂) and multipolar Gaussian charge distributions may also be of interest for modeling electrostatic interactions and short-range exchange and charge penetration effects in classical force fields and model potentials. With the presented technique, a 4894 × 4894 univ-JKFIT Coulomb matrix J_AB = (A|1/r₁₂|B) (183 MiB) can be computed in 50 ms on a Q2'2018 notebook CPU, without any screening or approximations.

Molecular Dynamics Simulations of Ionic Liquids and Electrolytes Using Polarizable Force Fields

Many applications in chemistry, biology, and energy storage/conversion research rely on molecular simulations to provide fundamental insight into structural and transport properties of materials with high ionic concentrations. Whether the system is comprised entirely of ions, like ionic liquids, or is a mixture of a polar solvent with a salt, e.g., liquid electrolytes for battery applications, the presence of ions in these materials results in strong local electric fields polarizing solvent molecules and large ions. To predict properties of such systems from molecular simulations often requires either explicit or mean-field inclusion of the influence of polarization on electrostatic interactions. In this manuscript, we review the pros and cons of different treatments of polarization ranging from the mean-field approaches to the most popular explicit polarization models in molecular dynamics simulations of ionic materials. For each method, we discuss their advantages and disadvantages and emphasize key assumptions as well as their adjustable parameters. Strategies for the development of polarizable models are presented with a specific focus on extracting atomic polarizabilities. Finally, we compare simulations using polarizable and nonpolarizable models for several classes of ionic systems, discussing the underlying physics that each approach includes or ignores, implications for implementation and computational efficiency, and the accuracy of properties predicted by these methods compared to experiments.

Steric "attraction": not by dispersion alone

Non-covalent interactions between neutral, sterically hindered organic molecules generally involve a strong stabilizing contribution from dispersion forces that in many systems turns the 'steric repulsion' into a 'steric attraction'. In addition to London dispersion, such systems benefit from electrostatic stabilization, which arises from a short-range effect of charge penetration and gets bigger with increasing steric bulk. In the present work, we quantify this contribution for a diverse set of molecular cores, ranging from unsubstituted benzene and cyclohexane to their derivatives carrying tert-butyl, phenyl, cyclohexyl and adamantyl substituents. While the importance of electrostatic interactions in the dimers of sp2-rich (e.g., π-conjugated) cores is well appreciated, less polarizable assemblies of sp3-rich systems with multiple short-range CH···HC contacts between the bulky cyclohexyl and adamantyl moieties are also significantly influenced by electrostatics. Charge penetration is drastically larger in absolute terms for the sp2-rich cores, but still has a non-negligible effect on the sp3-rich dimers, investigated herein, both in terms of their energetics and equilibrium interaction distances. These results emphasize the importance of this electrostatic effect, which has so far been less recognized in aliphatic systems compared to London dispersion, and are therefore likely to have implications for the development of force fields and methods for crystal structure prediction.

Computational Simulations of DNA Polymerases: Detailed Insights on Structure/Function/Mechanism from Native Proteins to Cancer Variants

Genetic information is vital in the cell cycle of DNA-based organisms. DNA polymerases (DNA Pols) are crucial players in transactions dealing with these processes. Therefore, the detailed understanding of the structure, function, and mechanism of these proteins has been the focus of significant effort. Computational simulations have been applied to investigate various facets of DNA polymerase structure and function. These simulations have provided significant insights over the years. This perspective presents the results of various computational studies that have been employed to research different aspects of DNA polymerases including detailed reaction mechanism investigation, mutagenicity of different metal cations, possible factors for fidelity synthesis, and discovery/functional characterization of cancer-related mutations on DNA polymerases.

Capturing Many-Body Interactions with Classical Dipole Induction Models

The nonadditive many-body interactions are significant for structural and thermodynamic properties of condensed phase systems. In this work we examined the many-body interaction energy of a large number of common organic/biochemical molecular clusters, which consist of 18 chemical species and cover nine common organic elements, using the Møller-Plesset perturbation theory to the second order (MP2) [ Møller et al. Phys. Rev. 1934 , 46 , 618 . ]. We evaluated the capability of Thole-based dipole induction models to capture the many-body interaction energy. Three models were compared: the original model and parameters used by the AMOEBA force field, a variation of this original model where the damping parameters have been reoptimized to MP2 data, and a third model where the damping function form applied to the permanent electric field is modified. Overall, we find the simple classical atomic dipole models are able to capture the 3- and 4-body interaction energy across a wide variety of organic molecules in various intermolecular configurations. With modified Thole models, it is possible to further improve the agreement with MP2 results. These models were also tested on systems containing metal/halogen ions to examine the accuracy and transferability. This work suggests that the form of damping function applied to the permanent electrostatic field strongly affects the distance dependence of polarization energy at short intermolecular separations.

Charge-dependent many-body exchange and dispersion interactions in combined QM/MM simulations

Accurate modeling of the molecular environment is critical in condensed phase simulations of chemical reactions. Conventional quantum mechanical/molecular mechanical (QM/MM) simulations traditionally model non-electrostatic non-bonded interactions through an empirical Lennard-Jones (LJ) potential which, in violation of intuitive chemical principles, is bereft of any explicit coupling to an atom's local electronic structure. This oversight results in a model whereby short-ranged exchange-repulsion and long-ranged dispersion interactions are invariant to changes in the local atomic charge, leading to accuracy limitations for chemical reactions where significant atomic charge transfer can occur along the reaction coordinate. The present work presents a variational, charge-dependent exchange-repulsion and dispersion model, referred to as the charge-dependent exchange and dispersion (QXD) model, for hybrid QM/MM simulations. Analytic expressions for the energy and gradients are provided, as well as a description of the integration of the model into existing QM/MM frameworks, allowing QXD to replace traditional LJ interactions in simulations of reactive condensed phase systems. After initial validation against QM data, the method is demonstrated by capturing the solvation free energies of a series of small, chlorine-containing compounds that have varying charge on the chlorine atom. The model is further tested on the SN2 attack of a chloride anion on methylchloride. Results suggest that the QXD model, unlike the traditional LJ model, is able to simultaneously obtain accurate solvation free energies for a range of compounds while at the same time closely reproducing the experimental reaction free energy barrier. The QXD interaction model allows explicit coupling of atomic charge with many-body exchange and dispersion interactions that are related to atomic size and provides a more accurate and robust representation of non-electrostatic non-bonded QM/MM interactions.

Modeling Molecular Interactions in Water: From Pairwise to Many-Body Potential Energy Functions

Almost 50 years have passed from the first computer simulations of water, and a large number of molecular models have been proposed since then to elucidate the unique behavior of water across different phases. In this article, we review the recent progress in the development of analytical potential energy functions that aim at correctly representing many-body effects. Starting from the many-body expansion of the interaction energy, specific focus is on different classes of potential energy functions built upon a hierarchy of approximations and on their ability to accurately reproduce reference data obtained from state-of-the-art electronic structure calculations and experimental measurements. We show that most recent potential energy functions, which include explicit short-range representations of two-body and three-body effects along with a physically correct description of many-body effects at all distances, predict the properties of water from the gas to the condensed phase with unprecedented accuracy, thus opening the door to the long-sought "universal model" capable of describing the behavior of water under different conditions and in different environments.

Beyond Born-Mayer: Improved Models for Short-Range Repulsion in ab Initio Force Fields

Short-range repulsion within intermolecular force fields is conventionally described by either Lennard-Jones (A/r(12)) or Born-Mayer (A exp(-Br)) forms. Despite their widespread use, these simple functional forms are often unable to describe the interaction energy accurately over a broad range of intermolecular distances, thus creating challenges in the development of ab initio force fields and potentially leading to decreased accuracy and transferability. Herein, we derive a novel short-range functional form based on a simple Slater-like model of overlapping atomic densities and an iterated stockholder atom (ISA) partitioning of the molecular electron density. We demonstrate that this Slater-ISA methodology yields a more accurate, transferable, and robust description of the short-range interactions at minimal additional computational cost compared to standard Lennard-Jones or Born-Mayer approaches. Finally, we show how this methodology can be adapted to yield the standard Born-Mayer functional form while still retaining many of the advantages of the Slater-ISA approach.

LICHEM: A QM/MM program for simulations with multipolar and polarizable force fields

We introduce an initial implementation of the LICHEM software package. LICHEM can interface with Gaussian, PSI4, NWChem, TINKER, and TINKER-HP to enable QM/MM calculations using multipolar/polarizable force fields. LICHEM extracts forces and energies from unmodified QM and MM software packages to perform geometry optimizations, single-point energy calculations, or Monte Carlo simulations. When the QM and MM regions are connected by covalent bonds, the pseudo-bond approach is employed to smoothly transition between the QM region and the polarizable force field. A series of water clusters and small peptides have been employed to test our initial implementation. The results obtained from these test systems show the capabilities of the new software and highlight the importance of including explicit polarization. © 2016 Wiley Periodicals, Inc.

Scalable improvement of SPME multipolar electrostatics in anisotropic polarizable molecular mechanics using a general short-range penetration correction up to quadrupoles

We propose a general coupling of the Smooth Particle Mesh Ewald SPME approach for distributed multipoles to a short-range charge penetration correction modifying the charge-charge, charge-dipole and charge-quadrupole energies. Such an approach significantly improves electrostatics when compared to ab initio values and has been calibrated on Symmetry-Adapted Perturbation Theory reference data. Various neutral molecular dimers have been tested and results on the complexes of mono- and divalent cations with a water ligand are also provided. Transferability of the correction is adressed in the context of the implementation of the AMOEBA and SIBFA polarizable force fields in the TINKER-HP software. As the choices of the multipolar distribution are discussed, conclusions are drawn for the future penetration-corrected polarizable force fields highlighting the mandatory need of non-spurious procedures for the obtention of well balanced and physically meaningful distributed moments. Finally, scalability and parallelism of the short-range corrected SPME approach are addressed, demonstrating that the damping function is computationally affordable and accurate for molecular dynamics simulations of complex bio- or bioinorganic systems in periodic boundary conditions.

Force Field for Mg(2+), Mn(2+), Zn(2+), and Cd(2+) Ions That Have Balanced Interactions with Nucleic Acids

Divalent metal ions are of fundamental importance to the function and folding of nucleic acids. Divalent metal ion-nucleic acid interactions are complex in nature and include both territorial and site specific binding. Commonly employed nonbonded divalent ion models, however, are often parametrized against bulk ion properties and are subsequently utilized in biomolecular simulations without considering any data related to interactions at specific nucleic acid sites. Previously, we assessed the ability of 17 different nonbonded Mg(2+) ion models to reproduce different properties of Mg(2+) in aqueous solution including radial distribution functions, solvation free energies, water exchange rates, and translational diffusion coefficients. In the present work, we depart from the recently developed 12-6-4 potential models for divalent metal ions developed by Li and Merz and tune the pairwise parameters for Mg(2+), Mn(2+), Zn(2+), and Cd(2+) binding dimethyl phosphate, adenosine, and guanosine in order to reproduce experimental site specific binding free energies derived from potentiometric pH titration data. We further apply these parameters to investigate a metal ion migration previously proposed to occur during the catalytic reaction of the hammerhead ribozyme. The new parameters are shown to be accurate and balanced for nucleic acid binding in comparison with available experimental data and provide an important tool for molecular dynamics and free energy simulations of nucleic acids where these ions may exhibit different binding modes.

Multipolar Ewald methods, 1: theory, accuracy, and performance

The Ewald, Particle Mesh Ewald (PME), and Fast Fourier–Poisson (FFP) methods are developed for systems composed of spherical multipole moment expansions. A unified set of equations is derived that takes advantage of a spherical tensor gradient operator formalism in both real space and reciprocal space to allow extension to arbitrary multipole order. The implementation of these methods into a novel linear-scaling modified “divide-and-conquer” (mDC) quantum mechanical force field is discussed. The evaluation times and relative force errors are compared between the three methods, as a function of multipole expansion order. Timings and errors are also compared within the context of the quantum mechanical force field, which encounters primary errors related to the quality of reproducing electrostatic forces for a given density matrix and secondary errors resulting from the propagation of the approximate electrostatics into the self-consistent field procedure, which yields a converged, variational, but nonetheless approximate density matrix. Condensed-phase simulations of an mDC water model are performed with the multipolar PME method and compared to an electrostatic cutoff method, which is shown to artificially increase the density of water and heat of vaporization relative to full electrostatic treatment.