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.     

Implementation of extended Lagrangian dynamics in GROMACS for polarizable simulations using the classical Drude oscillator model

Explicit treatment of electronic polarization in empirical force fields used for molecular dynamics simulations represents an important advancement in simulation methodology. A straightforward means of treating electronic polarization in these simulations is the inclusion of Drude oscillators, which are auxiliary, charge‐carrying particles bonded to the cores of atoms in the system. The additional degrees of freedom make these simulations more computationally expensive relative to simulations using traditional fixed‐charge (additive) force fields. Thus, efficient tools are needed for conducting these simulations. Here, we present the implementation of highly scalable algorithms in the GROMACS simulation package that allow for the simulation of polarizable systems using extended Lagrangian dynamics with a dual Nosé–Hoover thermostat as well as simulations using a full self‐consistent field treatment of polarization. The performance of systems of varying size is evaluated, showing that the present code parallelizes efficiently and is the fastest implementation of the extended Lagrangian methods currently available for simulations using the Drude polarizable force field.     

Toward accurate solvation dynamics of lanthanides and actinides in water using polarizable force fields: from gas-phase energetics to hydration free energies

In this contribution, we focused on the use of polarizable force fields to model the structural, energetic, and thermodynamical properties of lanthanides and actinides in water. In a first part, we chose the particular case of the Th(IV) cation to demonstrate the capabilities of the AMOEBA polarizable force field to reproduce both reference ab initio gas-phase energetics and experimental data including coordination numbers and radial distribution functions. Using such model, we predicted the first polarizable force field estimate of Th(IV) solvation free energy, which accounts for −1,638 kcal/mol. In addition, we proposed in a second part of this work a full extension of the SIBFA (Sum of Interaction Between Fragments Ab initio computed) polarizable potential to lanthanides (La(III) and Lu(III)) and to actinides (Th(IV)) in water. We demonstrate its capabilities to reproduce all ab initio contributions as extracted from energy decomposition analysis computations, including many-body charge transfer and discussed its applicability to extended molecular dynamics and its parametrization on high-level post-Hartree–Fock data.     

Comparing polarizable force fields to ab initio calculations reveals nonclassical effects in condensed phases

In a recent work [Giese and York J. Chem. Phys. 120, 9903 (2004)] showed that many-body force field models based solely on pairwise Coulomb screening cannot simultaneously reproduce both gas-phase and condensed-phase polarizability limits. In particular, polarizable force fields applied to bifurcated water chains have been demonstrated to be overpolarized with respect to ab initio methods. This behavior was ascribed to the neglect of coupling between many-body exchange and polarization. In the present article we reproduce those results using different ab initio levels of theory and a polarizable model based on the chemical-potential equalization principle. Moreover we show that, when hydrogen-bond (H-bond) forming systems are considered, an additional nonclassical effect, i.e., intermolecular charge transfer, must be taken into account. Such effect, contrarily to that of coupling between many-body exchange and polarization, makes classical polarizable force fields underpolarized. In the case of water at standard conditions, being H-bonded geometries much more probable than the bifurcated ones, intermolecular charge transfer is the dominant effect.     

Variational calculation of quantum mechanical/molecular mechanical free energy with electronic polarization of solvent

Quantum mechanical/molecular mechanical (QM/MM) free energy calculation presents a significant challenge due to an excessive number of QM calculations. A useful approach for reducing the computational cost is that based on the mean field approximation to the QM subsystem. Here, we describe such a mean-field QM/MM theory for electronically polarizable systems by starting from the Hartree product ansatz for the total system and invoking a variational principle of free energy. The MM part is then recast to a classical polarizable model by introducing the charge response kernel. Numerical test shows that the potential of mean force (PMF) thus obtained agrees quantitatively with that obtained from a direct QM/MM calculation, indicating the utility of self-consistent mean-field approximation. Next, we apply the obtained method to prototypical reactions in several qualitatively different solvents and make a systematic comparison of polarization effects. The results show that in aqueous solution the PMF does not depend very much on the water models employed, while in nonaqueous solutions the PMF is significantly affected by explicit polarization. For example, the free energy barrier for a phosphoryl dissociation reaction in acetone and cyclohexane is found to increase by more than 10 kcal/mol when switching the solvent model from an empirical to explicitly polarizable one. The reason for this is discussed based on the parametrization of empirical nonpolarizable models.     

Polarizable simulations with second order interaction model (POSSIM) force field: Developing parameters for protein side‐chain analogues

A previously introduced polarizable simulations with second‐order interaction model (POSSIM) force field has been extended to include parameters for small molecules serving as models for peptide and protein side‐chains. Parameters have been fitted to permit reproducing many‐body energies, gas‐phase dimerization energies, and geometries and liquid‐phase heats of vaporization and densities. Quantum mechanical and experimental data have been used as the target for the fitting. The POSSIM framework combines accuracy of a polarizable force field and computational efficiency of the second‐order approximation of the full‐scale induced point dipole polarization formalism. The resulting parameters can be used for simulations of the parameterized molecules themselves or their analogues. In addition to this, these force field parameters are currently being used in further development of the POSSIM fast polarizable force field for proteins. © 2012 Wiley Periodicals, Inc.     

A polarizable force field for water using an artificial neural network

A force field for liquid water including polarization effects has been constructed using an artificial neural network (ANN). It is essential to include a many-body polarization effect explicitly into a potential energy function in order to treat liquid water which is dense and highly polar. The new potential energy function is a combination of empirical and nonempirical potentials. The TIP4P model was used for the empirical part of the potential. For the nonempirical part, an ANN with a back-propagation of error algorithm (BPNN) was introduced to reproduce the complicated many-body interaction energy surface from ab initio quantum mechanical calculations. BPNN, described in terms of a matrix, provides enough flexibility to describe the complex potential energy surface (PES). The structural and thermodynamic properties, calculated by isobaric-isothermal (constant-NPT) Monte Carlo simulations with the new polarizable force field for water, are compatible with experimental results. Thus, the simulation establishes the validity of using our estimated PES with a polarization effect for accurate predictions of liquid state properties. Applications of this approach are simple and systematic so that it can easily be applied to the development of other force fields besides the water-water system.     

Polarizable empirical force field for alkanes based on the classical Drude oscillator model

Recent extensions of potential energy functions used in empirical force field calculations have involved the inclusion of electronic polarizability. To properly include this extension into a potential energy function it is necessary to systematically and rigorously optimize the associated parameters based on model compounds for which extensive experimental data are available. In the present work, optimization of parameters for alkanes in a polarizable empirical force field based on a classical Drude oscillator is presented. Emphasis is placed on the development of parameters for CH3, CH2, and CH moieties that are directly transferable to long chain alkanes, as required for lipids and other biomolecules. It is shown that a variety of quantum mechanical and experimental target data are reproduced by the polarizable model. Notable is the proper treatment of the dielectric constant of pure alkanes by the polarizable force field, a property essential for the accurate treatment of, for example, hydrophobic solvation in lipid bilayers. The present alkane force field will act as the basis for the aliphatic moieties in an extensive empirical force field for biomolecules that includes the explicit treatment of electronic polarizability.     

Behavior of polarizable models in presence of strong electric fields. I. Origin of nonlinear effects in water point-charge systems

In the current opinion, the inclusion of polarization response in classical computer simulations is considered as one of the most important and urgent improvements to be implemented in modern empirical potential models. In this work we focus on the capability of polarizable models, based on the pairwise Coulomb interactions, to model systems where strong electric fields enter into play. As shown by Masia, Probst, and Rey (MPR) [in J. Chem. Phys. 121, 7362 (2004)], when a molecule interacts with point charges, polarizable models show underpolarization with respect to ab initio methods. We prove that this underpolarization, clearly related to nonlinear polarization effects, cannot be simply ascribed to the lack of hyperpolarization in the polarizable models, as suggested by MPR. Analysis of the electron-density rearrangement induced on a water molecule by a point charge reveals a twofold level of polarization response. One level involves intramolecular charge transfer on the whole molecular volume, with the related polarization exhibiting a seemingly linear behavior with the external electric field. The other nonlinear polarization level occurs only at strong electric fields and is found to be strictly correlated to the quantum-mechanical nature of the water molecule. The latter type of polarization has a local character, being limited to the space region of the water lone pairs.     

Effect of polarization on the opsin shift in rhodopsins. 2. Empirical polarization models for proteins

The explicit treatment of polarization as a many-body interaction in condensed-phase systems represents a current problem in empirical force-field development. Although a variety of efficient models for molecular polarization have been suggested, polarizable force fields are still far from common use nowadays. In this work, we consider interactive polarization models employing Thole's short-range damping scheme and assess them for application on polypeptides. Despite the simplicity of the model, we find mean polarizabilities and anisotropies of amino acid side chains in excellent agreement with MP2/cc-pVQZ benchmark calculations. Combined with restrained electrostatic potential (RESP) derived atomic charges, the models are applied in a quantum-mechanical/molecular-mechanical (QM/MM) approach. An iterative scheme is used to establish a self-consistent mutual polarization between the QM and MM moieties. This ansatz is employed to study the influence of the protein polarizability on calculated optical properties of the protonated Schiff base of retinal in rhodopsin (Rh), bacterio-rhodopsin (bR), and pharaonis sensory rhodopsin II (psRII). The shifts of the excitation energy due to the instantaneous polarization response of the protein to the charge transfer on the retinal chromophore are quantified using the high level ab initio multireference spectroscopy-oriented configuration interaction (SORCI) method. The results are compared with those of previously published QM1/QM2/MM models for bR and psRII.     

Combining the polarizable Drude force field with a continuum electrostatic Poisson–Boltzmann implicit solvation model

In this work, we have combined the polarizable force field based on the classical Drude oscillator with a continuum Poisson–Boltzmann/solvent‐accessible surface area (PB/SASA) model. In practice, the positions of the Drude particles experiencing the solvent reaction field arising from the fixed charges and induced polarization of the solute must be optimized in a self‐consistent manner. Here, we parameterized the model to reproduce experimental solvation free energies of a set of small molecules. The model reproduces well‐experimental solvation free energies of 70 molecules, yielding a root mean square difference of 0.8 kcal/mol versus 2.5 kcal/mol for the CHARMM36 additive force field. The polarization work associated with the solute transfer from the gas‐phase to the polar solvent, a term neglected in the framework of additive force fields, was found to make a large contribution to the total solvation free energy, comparable to the polar solute–solvent solvation contribution. The Drude PB/SASA also reproduces well the electronic polarization from the explicit solvent simulations of a small protein, BPTI. Model validation was based on comparisons with the experimental relative binding free energies of 371 single alanine mutations. With the Drude PB/SASA model the root mean square deviation between the predicted and experimental relative binding free energies is 3.35 kcal/mol, lower than 5.11 kcal/mol computed with the CHARMM36 additive force field. Overall, the results indicate that the main limitation of the Drude PB/SASA model is the inability of the SASA term to accurately capture non‐polar solvation effects. © 2018 Wiley Periodicals, Inc.     

Effective polarization energy of the naphthalene molecular crystal: a study on the polarizable force field

polarization energy of the localized charge in organic solids consists of electronic polarization energy, permanent electrostatic interactions, and inter/intra molecular relaxation energies. The effective electronic polarization energies for an electron/hole carrier were successfully estimated by AMOEBA polarizable force field in naphthalene molecular crystals. Both electronic polarization energy and permanent electrostatic interaction were in agreement with the preview experimental values. In addition, the influence of the multipoles from different distributed mutipole analysis (DMA) fitting options on the electrostatic interactions are discussed in this paper. We found that the multipoles obtained from Gauss-Hermite quadrature without diffuse function or grid-based quadrature with 0.325 Å H atomic radius will give reasonable electronic polarization energies and permanent interactions for electron and hole carriers.     

A polarizable empirical force field for molecular dynamics simulation of liquid hydrocarbons

Electronic polarizability is usually treated implicitly in molecular simulations, which may lead to imprecise or even erroneous molecular behavior in spatially electronically inhomogeneous regions of systems such as proteins, membranes, interfaces between compounds, or mixtures of solvents. The majority of available molecular force fields and molecular dynamics simulation software packages does not account explicitly for electronic polarization. Even the simplest charge‐on‐spring (COS) models have only been developed for few types of molecules. In this work, we report a polarizable COS model for cyclohexane, as this molecule is a widely used solvent, and for linear alkanes, which are also used as solvents, and are the precursors of lipids, amino acid side chains, carbohydrates, or nucleic acid backbones. The model is an extension of a nonpolarizable united‐atom model for alkanes that had been calibrated against experimental values of the density, the heat of vaporization and the Gibbs free energy of hydration for each alkane. The latter quantity was used to calibrate the parameters governing the interaction of the polarizable alkanes with water. Subsequently, the model was tested for other structural, thermodynamic, dielectric, and dynamic properties such as trans/gauche ratios, excess free energy, static dielectric permittivity, and self‐diffusion. A good agreement with the experimental data for a large set of properties for each considered system was obtained, resulting in a transferable set of polarizable force‐field parameters for CH₂, CH₃, and CH₄ moieties. © 2014 Wiley Periodicals, Inc.     

Continuum polarizable force field within the Poisson-Boltzmann framework

We have developed and tested a complete set of nonbonded parameters for a continuum polarizable force field. Our analysis shows that the new continuum polarizable model is consistent with B3LYP/cc-pVTZ in modeling electronic response upon variation of dielectric environment. Comparison with experiment also shows that the new continuum polarizable model is reasonable, with accuracy similar to that of B3LYP/cc-pVTZ in reproduction of dipole moments of selected organic molecules in the gas phase. We have further tested the validity to interchange the Amber van der Waals parameters between the explicit and continuum polarizable force fields with a series of dimers. It can be found that the continuum polarizable model agrees well with MP2/cc-pVTZ, with deviations in dimer binding energies less than 0.9 kcal/mol in the aqueous dielectric environment. Finally, we have optimized atomic cavity radii with respect to experimental solvation free energies of 177 training molecules. To validate the optimized cavity radii, we have tested these parameters against 176 test molecules. It is found that the optimized Poisson-Boltzmann atomic cavity radii transfer well from the training set to the test set, with an overall root-mean-square deviation of 1.30 kcal/mol, an unsigned average error of 1.07 kcal/mol, and a correlation coefficient of 92% for all 353 molecules in both the training and test sets. Given the development documented here, the next natural step is the construction of a full protein/nucleic acid force field within the new continuum polarization framework.     

A polarizable ellipsoidal force field for halogen bonds

The anisotropic effects and short‐range quantum effects are essential characters in the formation of halogen bonds. Since there are an array of applications of halogen bonds and much difficulty in modeling them in classical force fields, the current research reports solely the polarizable ellipsoidal force field (PEff) for halogen bonds. The anisotropic charge distribution was represented with the combination of a negative charged sphere and a positively charged ellipsoid. The polarization energy was incorporated by the induced dipole model. The resulting force field is “physically motivated,” which includes separate, explicit terms to account for the electrostatic, repulsion/dispersion, and polarization interaction. Furthermore, it is largely compatible with existing, standard simulation packages. The fitted parameters are transferable and compatible with the general AMBER force field. This PEff model could correctly reproduces the potential energy surface of halogen bonds at MP2 level. Finally, the prediction of the halogen bond properties of human Cathepsin L (hcatL) has been found to be in excellent qualitative agreement with the cocrystal structures. © 2013 Wiley Periodicals, Inc.     

Further insights in the ability of classical nonadditive potentials to model actinide ion–water interactions

Pursuing our efforts on the development of accurate classical models to simulate radionuclides in complex environments (Réal et al., J. Phys. Chem. A 2010, 114, 15913; Trumm et al. J. Chem. Phys. 2012, 136, 044509), this article places a large emphasis on the discussion of the influence of models/parameters uncertainties on the computed structural, dynamical, and temporal properties. Two actinide test cases, trivalent curium and tetravalent thorium, have been studied with three different potential energy functions, which allow us to account for the polarization and charge‐transfer effects occurring in hydrated actinide ion systems. The first type of models considers only an additive energy term for modeling ion/water charge‐transfer effects, whereas the other two treat cooperative charge‐transfer interactions with two different analytical expressions. Model parameters are assigned to reproduce high‐level ab initio data concerning only hydrated ion species in gas phase. For the two types of cooperative charge‐transfer models, we define two sets of parameters allowing or not to cancel out possible errors inherent to the force field used to model water/water interactions at the ion vicinity. We define thus five different models to characterize the solvation of each ion. For both ions, our cooperative charge‐transfer models lead to close results in terms of structure in solution: the coordination number is included within 8 and 9, and the mean ion/water oxygen distances are 2.45 and 2.49 Å, respectively, for Th(IV) and Cm(III). © 2012 Wiley Periodicals, Inc.     

Structure,properties, and nature of the BrF‐HX complexes: An Ab initio study

Structure and properties of complexes (energies and charge transfer) of complexes BrF‐HX (X = F, Cl, Br, I) have been investigated at the MP2/aug‐cc‐pVDZ (aug‐cc‐pVDZ‐pp basis sets for I) level. Two types of geometries (hydrogen‐bonded and halogen‐bonded) are observed. The calculated interaction energies show that the halogen bonded structures are more stable than the corresponding hydrogen‐bonded structures. To study the nature of the intermolecular interactions, symmetry‐adapted perturbation theory (SAPT) energy decomposition analysis reveals that the BrF‐HX complexes are dominantly electrostatic in nature. © 2011 Wiley Periodicals, Inc. Int J Quantum Chem, 2011     

Developing multisite empirical force field models for Pt(II) and cisplatin

We have developed empirical force field parameters for Pt(II) and cisplatin. Two force field frameworks were used—modified OPLS‐AA and our second‐order polarizable POSSIM. A seven‐site model was used for the Pt(II) ion. The goal was to create transferable parameter sets compatible with the force field models for proteins and general organic compounds. A number of properties of the Pt(II) ion and its coordination compounds have been considered, including geometries and energies of the complexes, hydration free energy, and radial distribution functions in water. Comparison has been made with experimental and quantum mechanical results. We have demonstrated that both versions are generally capable of reproducing key properties of the system, but the second‐order polarizable POSSIM formalism permits more accurate quantitative results to be obtained. For example, the energy of formation of cisplatin as calculated with the modified OPLS‐AA exhibited an error of 9.9%, while the POSSIM error for the same quantity was 6.2%. The produced parameter sets are transferable and suitable to be used in protein‐metal binding simulations in which position or even coordination of the ion does not have to be constrained using preexisting knowledge. © 2016 Wiley Periodicals, Inc.     

Recent applications and developments of charge equilibration force fields for modeling dynamical charges in classical molecular dynamics simulations

With the continuing advances in computational hardware and novel force fields constructed using quantum mechanics, the outlook for non-additive force fields is promising. Our work in the past several years has demonstrated the utility of polarizable force fields, in our hands those based on the charge equilibration formalism, for a broad range of physical and biophysical systems. We have constructed and applied polarizable force fields for small molecules, proteins, lipids, and lipid bilayers and recently have begun work on carbohydrate force fields. The latter area has been relatively untouched by force field developers with particular focus on polarizable, non-additive interaction potential models. In this review of our recent work, we discuss the formalism we have adopted for implementing the charge equilibration method for phase-dependent polarizable force fields, lipid molecules, and small-molecule carbohydrates. We discuss the methodology, related issues, and briefly discuss results from recent applications of such force fields.