A valence bond model for aqueous Cu(II) and Zn(II) ions in the AMOEBA polarizable force field

A general molecular mechanics (MM) model for treating aqueous Cu²⁺ and Zn²⁺ ions was developed based on valence bond (VB) theory and incorporated into the atomic multipole optimized energetics for biomolecular applications (AMOEBA) polarizable force field. Parameters were obtained by fitting MM energies to that computed by ab initio methods for gas‐phase tetra‐ and hexa‐aqua metal complexes. Molecular dynamics (MD) simulations using the proposed AMOEBA‐VB model were performed for each transition metal ion in aqueous solution, and solvent coordination was evaluated. Results show that the AMOEBA‐VB model generates the correct square‐planar geometry for gas‐phase tetra‐aqua Cu²⁺ complex and improves the accuracy of MM model energetics for a number of ligation geometries when compared to quantum mechanical (QM) computations. On the other hand, both AMOEBA and AMOEBA‐VB generate results for Zn²⁺–water complexes in good agreement with QM calculations. Analyses of the MD trajectories revealed a six‐coordination first solvation shell for both Cu²⁺ and Zn²⁺ ions in aqueous solution, with ligation geometries falling in the range reported by previous studies. © 2012 Wiley Periodicals, Inc.     

Polarizable interaction potential for molecular dynamics simulations of actinoids(III) in liquid water

In this work, we have developed a polarizable classical interaction potential to study actinoids(III) in liquid water. This potential has the same analytical form as was recently used for lanthanoid(III) hydration [M. Duvail, P. Vitorge, and R. Spezia, J. Chem. Phys. 130, 104501 (2009)]. The hydration structure obtained with this potential is in good agreement with the experimentally measured ion-water distances and coordination numbers for the first half of the actinoid series. In particular, the almost linearly decreasing water-ion distance found experimentally is replicated within the calculations, in agreement with the actinoid contraction behavior. We also studied the hydration of the last part of the series, for which no structural experimental data are available, which allows us to provide some predictive insights on these ions. In particular we found that the ion-water distance decreases almost linearly across the series with a smooth decrease of coordination number from nine to eight at the end.     

Refinement of the NMR structures of alpha-conotoxin MI using molecular dynamics simulation with explicit solvent water and a full molecular force field

Three NMR structures of alpha-conotoxin MI, a potent antagonist of the nicotinic acetylcholine receptor, have been refined using molecular dynamics (MD) simulation with explicit water. Although the convergence of the NMR structures of alpha-conotoxin MI was not sufficient to provide detailed structural features, the average structures obtained from MD simulations converged to one conformation, providing structural characteristics. The resulting structure was also found to be consistent with the results of amide proton-exchange experiments. These results demonstrate that MD simulation with explicit solvent water is very useful in refining NMR structures.     

Three-site and five-site fixed-charge water models compatible with AMOEBA force field

In a typical biomolecular simulation using Atomic Multipole Optimized Energetics for Biomolecular Applications (AMOEBA) force field, the vast majority molecules in the simulation box consist of water, and these water molecules consume the most CPU power due to the explicit mutual induction effect. To improve the computational efficiency, we here develop two new nonpolarizable water models (with flexible bonds and fixed charges) that are compatible with AMOEBA solute: the 3-site AW3C and 5-site AW5C. To derive the force-field parameters for AW3C and AW5C, we fit to six experimental liquid thermodynamic properties: liquid density, enthalpy of vaporization, dielectric constant, isobaric heat capacity, isothermal compressibility and thermal expansion coefficient, at a broad range of temperatures from 261.15 to 353.15 K under 1.0 atm pressure. We further validate our AW3C and AW5C water models by showing that they can well reproduce the radial distribution function g(r), self-diffusion constant D, and hydration free energy from the AMOEBA03 water model and the experimental observations. Furthermore, we show that our AW3C and AW5C water models can greatly accelerate (>5 times) the bulk water as well as biomolecular simulations when compared to AMOEBA water. Specifically, we demonstrate that the applications of AW3C and AW5C water models to simulate a DNA duplex lead to a threefold acceleration, and in the meanwhile well maintain the structural properties as the fully polarizable AMOEBA water. We expect that our AW3C and AW5C water models hold great promise to be widely applied to simulate complex bio-molecules using the AMOEBA force field.     

New all-atom force field for molecular dynamics simulation of an AlPO(4)-34 molecular sieve

A force field of the triclinic framework of AlPO(4)-34, important in methanol-hydrocarbon conversion reactions, was developed using an empirical potential function. Molecular dynamics simulation of an AlPO(4)-34 triclinic framework segment of 1216 atoms, containing the template molecules isopropylamine and water, was performed with explicit consideration of atomic charges. The average RMS difference between instantaneous positions of the framework atoms during 1 ns simulation and their positions in the structure determined from single crystal X-ray diffraction was calculated, and the average structure of the flexible framework was determined. The computed Debye-Waller factors and simulated FTIR spectra are in good agreement with the experimental data. The new force field permits detailed molecular dynamics simulations of flexible, charged aluminophosphate molecular sieves which should lead to a better understanding of the catalytic processes and the crucial role played by templating molecules.     

Predicting water uptake in poly(perfluorosulfonic acids) using force field simulation methods

Free energy perturbation methods were applied to predict water contents in hydrated poly(perfluorosulfonic acids) (PPFSA). The simulations were based on the TEAM force field which was derived from quantum mechanical data calculated for small molecules using density functional theory (DFT) and thermodynamic data of molecular liquids and crystal. The equilibrium water contents in three PPFSA polymers (Nafion-117, Nafion-115 and Hyflon) were predicted by evaluating excess chemical potentials of water in hydrated polymers and in pure water. High level of precision measured by average uncertainty of ca. 0.1 kcal mol(-1), and accuracy in terms of deviation from experimental data by ca. 0.2 kcal mol(-1) were obtained in the predicted excess chemical potentials. The predicted amounts of water uptake agree well with experimental values. In addition, the equilibrium and dynamic properties of hydrated Nafion-117 were calculated and the results agree well with the existing experimental and computational data. The entropy and enthalpy contributions in the calculated excess chemical potentials are analyzed and the results are consistent with intuition. A linear correlation between the entropies and enthalpies is identified for the systems studied, which indicates that just increasing the interaction energies between water and host materials does not guarantee enhancement of the water uptake.     

Effective force field for liquid hydrogen fluoride from ab initio molecular dynamics simulation using the force-matching method

A recently developed force-matching method for obtaining effective force fields for condensed matter systems from ab initio molecular dynamics (MD) simulations has been applied to fit a simple nonpolarizable two-site pairwise force field for liquid hydrogen fluoride. The ab initio MD in this case was a Car-Parrinello (CP) MD simulation of 64 HF molecules at nearly ambient conditions within the Becke-Lee-Yang-Parr approximation to the electronic density functional theory. The force-matching procedure included a fit of short-ranged nonbonded forces, bonded forces, and atomic partial charges. The performance of the force-match potential was examined for the gas-phase dimer and for the liquid phase at various temperatures. The model was able to reproduce correctly the bent structure and energetics of the gas-phase dimer, while the results for the structural properties, self-diffusion, vibrational spectra, density, and thermodynamic properties of liquid HF were compared to both experiment and the CP MD simulation. The force-matching model performs well in reproducing nearly all of the liquid properties as well as the aggregation behavior at different temperatures. The model is computationally cheap and compares favorably to many more computationally expensive potential energy functions for liquid HF.     

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.     

Simulation of phase separation in alcohol/water mixtures using two-body force field and standard molecular dynamics

Standard molecular dynamics simulations have been carried out on pure alcohols and alcohol/water mixtures. A simple atom-atom force field consisting of Lennard-Jones potentials plus coulombic terms over atomic point charges, but without explicit polarization terms, has been specifically fitted to reproduce several experimental properties of the pure alcohols, and has been used for mixtures by developing combination rules with the TIP3P water model. Densities, enthalpies of vaporization, radial distribution functions, self-diffusion coefficients, and rotational correlation functions of the pure alcohols are well reproduced and compare favorably with those from more sophisticated force fields. Some key aspects of the phase behaviour are correctly reproduced by the molecular dynamics simulation, showing a distinct demixing process for the n-butanol/water mixture as opposed to the stability of the t-butanol/water mixtures. The results demonstrate the ability of a molecular dynamics simulation, even in its standard form and with easily accessible time ranges, but with a carefully optimized force field, to simulate and, to a certain extent, predict the properties of binary mixtures.     

United atom force field for phospholipid membranes: Constant pressure molecular dynamics simulation of dipalmitoylphosphatidicholine/water system

We refined the united atom field for the simulations of phospholipid membranes. To validate this potential we performed 1000-ps constant pressure simulation of a dipalmitoylphosphatidicholine (DPPC) bilayer at T=50° C. The average area per head group (61.6±0.6) Ų obtained in our simulation agrees well with the measured one of (62.9±1.3) Ų. The calculated S_CD order parameters for the Sn-2 hydrocarbon tail also display a good agreement with the experiment. The conformations of head groups in our simulations of the liquid crystal phase are different than the ones observed in the crystal structure. ©1999 John Wiley & Sons, Inc. J Comput Chem 20, 531–545, 1999     

Polarizable coarse‐grained models for molecular dynamics simulation of liquid cyclohexane

Force field parameters for polarizable coarse‐grained (CG) supra‐atomic models of liquid cyclohexane are proposed. Two different bead sizes were investigated, one representing two fine‐grained (FG) CH₂r united atoms of the cyclohexane ring, and one representing three FG CH₂r united atoms. Electronic polarizability is represented by a massless charge‐on‐spring particle connected to each CG bead. The model parameters were calibrated against the experimental density and heat of vaporization of liquid cyclohexane, and the free energy of cyclohexane hydration. Both models show good agreement with thermodynamic properties of cyclohexane, yet overestimate the self‐diffusion. The dielectric properties of the polarizable models agree very well with experiment. © 2015 Wiley Periodicals, Inc.     

Polarizable molecular mechanics studies of Cu(I)/Zn(II) superoxide dismutase: Bimetallic binding site and structured waters

The existence of a network of structured waters in the vicinity of the bimetallic site of Cu/Zn‐superoxide dismutase (SOD) has been inferred from high‐resolution X‐ray crystallography. Long‐duration molecular dynamics (MD) simulations could enable to quantify the lifetimes and possible interchanges of these waters between themselves as well as with a ligand diffusing toward the bimetallic site. The presence of several charged or polar ligands makes it necessary to resort to second‐generation polarizable potentials. As a first step toward such simulations, we benchmark in this article the accuracy of one such potential, sum of interactions between fragments Ab initio computed (SIBFA), by comparisons with quantum mechanics (QM) computations. We first consider the bimetallic binding site of a Cu/Zn‐SOD, in which three histidines and a water molecule are bound to Cu(I) and three histidines and one aspartate are bound to Zn(II). The comparisons are made for different His6 complexes with either one or both cations, and either with or without Asp and water. The total net charges vary from zero to three. We subsequently perform preliminary short‐duration MD simulations of 296 waters solvating Cu/Zn‐SOD. Six representative geometries are selected and energy‐minimized. Single‐point SIBFA and QM computations are then performed in parallel on model binding sites extracted from these six structures, each of which totals 301 atoms including the closest 28 waters from the Cu metal site. The ranking of their relative stabilities as given by SIBFA is identical to the QM one, and the relative energy differences by both approaches are fully consistent. In addition, the lowest‐energy structure, from SIBFA and QM, has a close overlap with the crystallographic one. The SIBFA calculations enable to quantify the impact of polarization and charge transfer in the ranking of the six structures. Five structural waters, which connect Arg141 and Glu131, are endowed with very high dipole moments (2.7–3.0 Debye), equal and larger than the one computed by SIBFA in ice‐like arrangements (2.7 D).     

Enthalpies of mixing predicted using molecular dynamics simulations and OPLS force field

Enthalpy of mixing (EOM) is one of the most basic thermodynamic properties of mixtures. To assess feasibility of predicting EOM using force field simulation methods, fifteen (15) representative binary mixtures were investigated using MD simulations based on OPLS and TIP4P force fields. The simulation conditions and errors were carefully examined. The precision level of 0.04 kJ/mol was obtained for calculated EOM data. However, the predictions, measured by deviations from experimental data, were only qualitatively correct. Among various factors influencing the accuracy of predictions, force field quality representing interactions among different molecules plays the most significant role. Using methanol/benzene and ethanol/benzene as examples, we demonstrated that non-additive interaction terms between polarizable atoms can be used to significantly improve the quality of predictions. In addition, it appears that charge-dependent LJ parameters are required in order to represent the polarization effects accurately.     

Automation of AMOEBA polarizable force field parameterization for small molecules

A protocol to generate parameters for the AMOEBA polarizable force field for small organic molecules has been established, and polarizable atomic typing utility, Poltype, which fully automates this process, has been implemented. For validation, we have compared with quantum mechanical calculations of molecular dipole moments, optimized geometry, electrostatic potential, and conformational energy for a variety of neutral and charged organic molecules, as well as dimer interaction energies of a set of amino acid side chain model compounds. Furthermore, parameters obtained in gas phase are substantiated in liquid-phase simulations. The hydration free energy (HFE) of neutral and charged molecules have been calculated and compared with experimental values. The RMS error for the HFE of neutral molecules is less than 1 kcal/mol. Meanwhile, the relative error in the predicted HFE of salts (cations and anions) is less than 3% with a correlation coefficient of 0.95. Overall, the performance of Poltype is satisfactory and provides a convenient utility for applications such as drug discovery. Further improvement can be achieved by the systematic study of various organic compounds, particularly ionic molecules, and refinement and expansion of the parameter database.     

Study of amino- and pyridyl-containing zinc(II) and cadmium(II) complexes by molecular mechanics using a force field based on the Gillespie-Kepert model

A new approach combining the molecular mechanics (MM) method and the Gillespie-Kepert model was applied to calculate the geometry and strain energy of zinc(II) and cadmium(II) complexes with amino- and pyridyl-containing ligands. High accuracy of calculations of the geometry was demonstrated for more than 20 complexes of these metals. Typical r.m.s. deviations between the calculated and experimental values (X-ray diffraction analysis) were 0.02 Å for bond lengths, 2° for bond angles, and 4° for torsion angles. The size-match selectivity of several macrocycles and polydentate open-chain ligands was studied. Correlations between the calculated strain energies of metal complexes and the experimental values of their stability constants and enthalpies of formation are discussed.     

CHARMM force field and molecular dynamics simulations of protonated polyethylenimine

As a gene delivery vector, polyethylenimine (PEI) shows one of the highest transfection efficiencies, while effectively protecting DNA from enzyme degradation. The distinctive charge pattern of protonated PEI is widely considered responsible for fundamental process such as DNA condensation into PEI/DNA polyplexes (which are able to enter cells via endocytosis), proton sponge effect (which triggers the release of polyplexes from endosome), and release of DNA from polyplexes (to be further processed inside the nucleus). Our investigations are largely motivated by the crucial need for a realistic molecular mechanics force field (FF) for PEI, and, accordingly, we focus on two major issues: (1) development of a new atomistic (CHARMM) FF for PEI in different protonation states, rigorously derived from high‐quality ab initio calculations performed on model polymers, and (2) molecular dynamics investigations of solvated PEI, providing a detailed picture of the dynamic structuring thereof in dependence on their size and protonation state. The modeled PEI chains are essentially described in terms of gyration radius, end‐to‐end distance, persistence length, radial distribution functions, coordination numbers, and diffusion coefficients. They turn out to be more rigid than in other computational studies and we find diffusion coefficients in fair agreement with experimental data. The developed atomistic FF proves adequate for the realistic modeling of the size and protonation behavior of linear PEI, either as individual chains or composing polyplexes. © 2017 Wiley Periodicals, Inc.     

Development of force field parameters for molecular simulation of polylactide

Polylactide is a biodegradable polymer that is widely used for biomedical applications, and it is a replacement for some petroleum based polymers in applications that range from packaging to carpeting. Efforts to characterize and further enhance polylactide based systems using molecular simulations have to this point been hindered by the lack of accurate atomistic models for the polymer. Thus, we present force field parameters specifically suited for molecular modeling of PLA. The model, which we refer to as PLAFF3, is based on a combination of the OPLS and CHARMM force fields, with modifications to bonded and nonbonded parameters. Dihedral angle parameters were adjusted to reproduce DFT data using newly developed CMAP dihedral cross terms, and the model was further adjusted to reproduce experimentally resolved crystal structure conformations, melt density, volume expansivity, and the glass transition temperature of PLA. We recommend the use of PLAFF3 in modeling PLA in its crystalline or amorphous states and have provided the necessary input files required for the publicly available molecular dynamics code GROMACS.     

Evaluation of solvation free energies for small molecules with the AMOEBA polarizable force field

The effects of electronic polarization in biomolecular interactions will differ depending on the local dielectric constant of the environment, such as in solvent, DNA, proteins, and membranes. Here the performance of the AMOEBA polarizable force field is evaluated under nonaqueous conditions by calculating the solvation free energies of small molecules in four common organic solvents. Results are compared with experimental data and equivalent simulations performed with the GAFF pairwise‐additive force field. Although AMOEBA results give mean errors close to “chemical accuracy,” GAFF performs surprisingly well, with statistically significantly more accurate results than AMOEBA in some solvents. However, for both models, free energies calculated in chloroform show worst agreement to experiment and individual solutes are consistently poor performers, suggesting non‐potential‐specific errors also contribute to inaccuracy. Scope for the improvement of both potentials remains limited by the lack of high quality experimental data across multiple solvents, particularly those of high dielectric constant. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Comparative study of Zn(II) and Cd(II) ions adsorption on charged carbon nano tubes: molecular dynamics approach

A comparative study was performed on adsorption process of Zn(II) and Cd(II) ions using charged carbon nano tubes (CNTs). The results show that the adsorption behavior of Zn(II) and Cd(II) ions have different patterns. To understand the reason for this difference, the ion–water complex properties and the reduced density of water around CNT were studied. The results reveal that the water molecules make high-, low-, and bulk regions around CNTs, which can affect the ion movement toward CNT. The ion–water complex characteristics during the adsorption demonstrated that the Cd(II)–water complex has a large number of exchange events between 4 and 8 coordinate numbers where Zn(II) has only 6 coordinate numbers during the adsorption process. The Cd(II)–water complex was found less stable and therefore has faster translational dynamics. These properties help cadmium ion to be adsorbed faster into CNT despite its heavier weight.