Design Of Next Generation Force Fields From AB Initio Computations: Beyond Point Charges Electrostatics

Exchange Repulsion in Quantum Mechanical/Effective Fragment Potential Excitation Energies: Beyond Polarizable Embedding

Hybrid quantum mechanical and molecular mechanical (QM/MM) approaches facilitate computational modeling of large biological and materials systems. Typically, in QM/MM, a small region of the system is modeled with an accurate quantum mechanical method and its surroundings with a more efficient alternative, such as a classical force field or the effective fragment potential (EFP). The reliability of QM/MM calculations depends largely on the treatment of interactions between the two subregions, also known as embedding. The polarizable embedding, which allows mutual polarization between solvent and solute, is considered to be essential for describing electronic excitations in polar solvents. In this work, we employ the QM/EFP model and extend the polarizable embedding by incorporating two short-range terms-a charge penetration correction to the electrostatic term and the exchange-repulsion term-both of which are modeled with one-electron contributions to the quantum Hamiltonian. We evaluate the accuracy of these terms by computing excitation energies across 37 molecular clusters consisting of biologically relevant chromophores surrounded by polar solvent molecules. QM/EFP excitation energies are compared to the fully quantum mechanical calculations with the configuration interaction singles (CIS) method. We find that the charge penetration correction diminishes the accuracy of the QM/EFP calculations. On the other hand, while the effect of exchange-repulsion is negligible for most ππ* transitions, the exchange-repulsion significantly improves description of nπ* transitions with blue solvatochromic shifts. As a result, addition of the exchange-repulsion term improves the overall accuracy of QM/EFP. Performances of QM/EFP models remain similar when excitation energies are modeled with cc-pVDZ and aug-cc-pVDZ basis sets.

Dispersion Interactions in QM/EFP

The dispersion energy term between quantum-mechanical (QM) and classical (represented by effective fragment potentials, EFP) subsystems is developed and implemented. A new formulation is based on long-range perturbation theory and uses dynamic polarizability tensors of the effective fragments and electric field integrals and orbital energies of the quantum-mechanical subsystem. No parametrization is involved. The accuracy of the QM-EFP dispersion energy is tested on a number of model systems; the average mean unsigned error is 0.8 kcal/mol or 13% with respect to the symmetry adapted perturbation theory on the S22 data set of noncovalent interactions. The computational cost of the dispersion energy computation is low compared to the self-consistent field calculation of the QM subsystem. The dispersion energy is sensitive to the level of theory employed for the QM part and to the electrostatic interactions in the system. The latter means that the dispersion interactions in the QM/EFP method are not purely two-body but have more complex many-body behavior.

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

We demonstrate as a proof of principle the capabilities of a novel hybrid MM'/MM polarizable force field to integrate short-range quantum effects in molecular mechanics (MM) through the use of Gaussian electrostatics. This lead to a further gain in accuracy in the representation of the first coordination shell of metal ions. It uses advanced electrostatics and couples two point dipole polarizable force fields, namely, the Gaussian electrostatic model (GEM), a model based on density fitting, which uses fitted electronic densities to evaluate nonbonded interactions, and SIBFA (sum of interactions between fragments ab initio computed), which resorts to distributed multipoles. To understand the benefits of the use of Gaussian electrostatics, we evaluate first the accuracy of GEM, which is a pure density-based Gaussian electrostatics model on a test Ca(II)-H2O complex. GEM is shown to further improve the agreement of MM polarization with ab initio reference results. Indeed, GEM introduces nonclassical effects by modeling the short-range quantum behavior of electric fields and therefore enables a straightforward (and selective) inclusion of the sole overlap-dependent exchange-polarization repulsive contribution by means of a Gaussian damping function acting on the GEM fields. The S/G-1 scheme is then introduced. Upon limiting the use of Gaussian electrostatics to metal centers only, it is shown to be able to capture the dominant quantum effects at play on the metal coordination sphere. S/G-1 is able to accurately reproduce ab initio total interaction energies within closed-shell metal complexes regarding each individual contribution including the separate contributions of induction, polarization, and charge-transfer. Applications of the method are provided for various systems including the HIV-1 NCp7-Zn(II) metalloprotein. S/G-1 is then extended to heavy metal complexes. Tested on Hg(II) water complexes, S/G-1 is shown to accurately model polarization up to quadrupolar response level. This opens up the possibility of embodying explicit scalar relativistic effects in molecular mechanics thanks to the direct transferability of ab initio pseudopotentials. Therefore, incorporating GEM-like electron density for a metal cation enable the introduction of nonambiguous short-range quantum effects within any point-dipole based polarizable force field without the need of an extensive parametrization.

Understanding the structure and electronic properties of Th4+-water complexes

We present a systematic quantum chemistry study of the [Th(H2O)n]4+ (n = 1 to 10) complexes to gain insight into their electronic structure and properties: the effect of the ligand distribution on the valence shells of the thorium(IV) ion is studied by means of the electron localization function (ELF) topological analysis. Particular care is given to the study of the mono-aqua complex both at its equilibrium geometry, using various tools such as energy decomposition analyses (EDA), and along its dissociation pathway. Indeed, as several electronic states cross the Th4 +-H2O0 ground state along the minimum energy path, we demonstrate that the diabatic representation implemented in MOLPRO is able to generate reference potential energy surfaces that will lead to the evaluation of diabatic dissociation curves. The calculated diabatic interaction energy curve will allow for a consistent parameterization of new generation force fields dedicated to heavy metals based on quantum chemistry.

Efficient optimization of van der Waals parameters from bulk properties

Due to the computational cost involved, when developing a force field for new compounds, one often avoids fitting van der Waals (vdW) terms, instead relying on a general force field based on the atom type. Here, we provide a novel approach to efficiently optimize vdW terms, based on both ab initio dimer energies and condensed phase properties. The approach avoids the computational challenges of searching the parameter space by using an extrapolation method to obtain a reliable difference quotient for the parameter derivatives based on the central difference. The derivatives are then used in an active-space optimization method which convergences quadratically. This method is applicable to polarizable and nonpolarizable force fields, although we focus on the parameterization of the AMBER force field. The scaling of the electrostatic potential (ESP) of the compounds is also studied. The algorithm is tested on 12 compounds, reducing the root mean squared error (RMSE) of the density from 0.061 g/cm(3) with GAFF parameters to 0.004 g/cm(3) , and the heat of vaporization from 1.13 to 0.05 kcal/mol. This is done with only four iterations of molecular dynamic runs. Using the optimized vdW parameters, the RMSE of the self-diffusion (Dself ) was reduced from 1.22 × 10(-9) to 0.78 × 10(-9) m(2) s(-1) and the RMSE of the hydration free energies (ΔGsolv ) was reduced from 0.30 to 0.26 kcal/mol. Scaling the ESP to improve dimer energies resulted in the Dself  RMSE improving to 0.77× 10(-9) m(2) s(-1) , but the ΔGsolv  worsened to 0.33 kcal/mol.

A variational linear-scaling framework to build practical, efficient next-generation orbital-based quantum force fields

We introduce a new hybrid molecular orbital/density-functional modified divide-and-conquer (mDC) approach that allows the linear-scaling calculation of very large quantum systems. The method provides a powerful framework from which linear-scaling force fields for molecular simulations can be developed. The method is variational in the energy, and has simple, analytic gradients and essentially no break-even point with respect to the corresponding full electronic structure calculation. Furthermore, the new approach allows intermolecular forces to be properly balanced such that non-bonded interactions can be treated, in some cases, to much higher accuracy than the full calculation. The approach is illustrated using the second-order self-consistent charge density-functional tight-binding model (DFTB2). Using this model as a base Hamiltonian, the new mDC approach is applied to a series of water systems, where results show that geometries and interaction energies between water molecules are greatly improved relative to full DFTB2. In order to achieve substantial improvement in the accuracy of intermolecular binding energies and hydrogen bonded cluster geometries, it was necessary to extend the DFTB2 model to higher-order atom-centered multipoles for the second-order self-consistent intermolecular electrostatic term. Using generalized, linear-scaling electrostatic methods, timings demonstrate that the method is able to calculate a water system of 3000 atoms in less than half of a second, and systems of up to one million atoms in only a few minutes using a conventional desktop workstation.

Many-body exchange-repulsion in polarizable molecular mechanics. I. Orbital-based approximations and applications to hydrated metal cation complexes

We have quantified the extent of the nonadditivity of the short-range exchange-repulsion energy, E(exch-rep), in several polycoordinated complexes of alkali, alkaline-earth, transition, and metal cations. This was done by performing ab initio energy decomposition analyses of interaction energies in these complexes. The magnitude of E(exch-rep(n-body, n > 2)) was found to be strongly cation-dependent, ranging from close to zero for some alkali metal complexes to about 6 kcal/mol for the hexahydrated Zn(2+) complex. In all cases, the cation-water molecules, E(exch-rep(three-body)), has been found to be the dominant contribution to many-body exchange-repulsion effects, higher order terms being negligible. As the physical basis of this effect is discussed, a three-center exponential term was introduced in the SIBFA (Sum of Interactions Between Fragments Ab initio computed) polarizable molecular mechanics procedure to model such effects. The three-body correction is added to the two-center (two-body) overlap-like formulation of the short-range repulsion contribution, E(rep), which is grounded on simplified integrals obtained from localized molecular orbital theory. The present term is computed on using mostly precomputed two-body terms and, therefore, does not increase significantly the computational cost of the method. It was shown to match closely E(three-body) in a series of test cases bearing on the complexes of Ca(2+), Zn(2+), and Hg(2+). For example, its introduction enabled to restore the correct tetrahedral versus square planar preference found from quantum chemistry calculations on the tetrahydrate of Hg(2+) and [Hg(H(2)O)(4)](2+).

Connecting the virtual world of computers to the real world of medicinal chemistry

Drug discovery involves the simultaneous optimization of chemical and biological properties, usually in a single small molecule, which modulates one of nature's most complex systems: the balance between human health and disease. The increased use of computer-aided methods is having a significant impact on all aspects of the drug-discovery and development process and with improved methods and ever faster computers, computer-aided molecular design will be ever more central to the discovery process.

Polarizable molecular dynamics simulation of Zn(II) in water using the AMOEBA force field

The hydration free energy, structure, and dynamics of the zinc divalent cation are studied using a polarizable force field in molecular dynamics simulations. Parameters for the Zn(2+) are derived from gas-phase ab initio calculation of Zn(2+)-water dimer. The Thole-based dipole polarization is adjusted based on the Constrained Space Orbital Variations (CSOV) calculation while the Symmetry Adapted Perturbation Theory (SAPT) approach is also discussed. The vdW parameters of Zn(2+) have been obtained by comparing the AMOEBA Zn(2+)-water dimerization energy with results from several theory levels and basis sets over a range of distances. Molecular dynamics simulations of Zn(2+) solvation in bulk water are subsequently performed with the polarizable force field. The calculated first-shell water coordination number, water residence time and free energy of hydration are consistent with experimental and previous theoretical values. The study is supplemented with extensive Reduced Variational Space (RVS) and Electron Localization Function (ELF) computations in order to unravel the nature of the bonding in Zn(2+)(H(2)O)(n) (n=1,6) complexes and to analyze the charge transfer contribution to the complexes. Results show that the importance of charge transfer decreases as the size of Zn-water cluster grows due to anticooperativity and to changes in the nature of the metal-ligand bonds. Induction could be dominated by polarization when the system approaches condensed-phase and the covelant effects are eliminated from the Zn(II)-water interaction. To construct an "effective" classical polarizable potential for Zn(2+) in bulk water, one should therefore avoid over-fitting to the ab initio charge transfer energy of Zn(2+)-water dimer. Indeed, in order to avoid overestimation of condensed-phase many-body effects, which is crucial to the transferability of polarizable molecular dynamics, charge transfer should not be included within the classical polarization contribution and should preferably be either incorporated in to the pairwise van der Waals contribution or treated explicitly.