Nonpolarizable Force Fields Research Articles

Prediction of hydration free energies for the SAMPL4 data set with the AMOEBA polarizable force field

Hydration free energy calculations are often used to validate molecular simulation methodologies and molecular mechanics force fields. We use the free-energy perturbation method together with the AMOEBA polarizable force field and the Poltype parametrization protocol to predict the hydration free energies of 52 molecules as part of the SAMPL4 blind challenge. For comparison, similar calculations are performed using the non-polarizable General Amber force field. Against our expectations, the latter force field gives the better results compared to experiment. One possible explanation is the sensitivity of the AMOEBA results to the conformation used for parametrization.

How does the electronic continuum model perform in the prediction of the surface tension of salt solutions?

The electronic continuum (EC) model uses a scaling of the charges of the ions in order to model implicitly the polarization into nonpolarizable models. This scaling procedure is applied here to two standard nonpolarizable force fields to investigate the salt concentration dependence of the surface tension and density of NaCl aqueous solutions. The composition of the interface and the orientation of the water molecules at the water surface are reported for different combinations of force fields.

Molecular dynamics study of response of liquid N,N-dimethylformamide to externally applied electric field using a polarizable force field.

The behavior of Liquid N,N-dimethylformamide subjected to a wide range of externally applied electric fields (from 0.001 V/nm to 1 V/nm) has been investigated through molecular dynamics simulation. To approach the objective the AMOEBA polarizable force field was extended to include the interaction of the external electric field with atomic partial charges and the contribution to the atomic polarization. The simulation results were evaluated with quantum mechanical calculations. The results from the present force field for the liquid at normal conditions were compared with the experimental and molecular dynamics results with non-polarizable and other polarizable force fields. The uniform external electric fields of higher than 0.01 V/nm have a significant effect on the structure of the liquid, which exhibits a variation in numerous properties, including molecular polarization, local cluster structure, rotation, alignment, energetics, and bulk thermodynamic and structural properties.

QM/MM Simulations of Mg and Zn Solvation

Mg(II) and Zn(II) ions serve a variety of the biological functions, such as signaling, catalysis, and structure. Despite similar ionic radii (0.72 A for Mg(II) and 0.74 A for Zn(II)), these ions bind to different types of ligands; Mg(II) predominantly binds to phosphates, forming species like Mg-ATP while Zn(II) typically binds to cysteine and histidine side chains. Some of the differences in the biochemistry of Mg(II) and Zn(II) ions are apparent in the solvation energies of these ions in water; despite the similar ionic radii, the hydration free energy of Zn(II) is 30 kcal/mol more favorable than Mg(II). In this work, the hydration properties of Mg(II) and Zn(II) have been explored computationally. The performance in the prediction of the energetics and solvation structure of Mg(II) and Zn(II) are compared. The CHARMM non-polarizable force field incorrectly predicts that Mg(II) has the more favorable solvation free energy. The Drude polarizable model improves the description, but still underestimates the relative solvation energy by 17 kcal/mol. Only a QM/MM simulation using the CHARMM/TURBOMOLE code was able to predict the relative solvation energy and structure of Mg(II) and Zn(II) in good agreement with experiment [1]. ETS-NOCV analysis indicates that this difference is due to increased water-to-ion charge transfer interactions in Zn(II) compare to the Mg(II) due to the higher Lewis acidity of Zn(II).[1] Riahi, S., Roux, B., Rowley, C.N. Can. J. Chem. 2013, 91(7), 552558.

Development of an Effective Polarizable Bond Method for Biomolecular Simulation

An effective polarizable bond (EPB) model has been developed for computer simulation of proteins. In this partial polarizable approach, all polar groups of amino acids are treated as polarizable, and the relevant polarizable parameters were determined by fitting to quantum calculated electrostatic properties of these polar groups. Extensive numerical tests on a diverse set of proteins (including 1IEP, 1MWE, 1NLJ, 4COX, 1PGB, 1K4C, 1MHN, 1UBQ, 1IGD) showed that this EPB model is robust in MD simulation and can correctly describe the structure and dynamics of proteins (both soluble and membrane proteins). Comparison of the computed hydrogen bond properties and dynamics of proteins with experimental data and with results obtained from the nonpolarizable force field clearly demonstrated that EPB can produce results in much better agreement with experiment. The averaged deviation of the simulated backbone N-H order parameter of the B3 immunoglobobulin-binding domain of streptococcal protein G from experimental observation is 0.0811 and 0.0332 for Amber99SB and EPB, respectively. This new model inherited the effective character of the classic force field and the fluctuating feature of previous polarizable models. Different from other polarizable models, the polarization cost energy is implicitly included in the present method. As a result, the present method avoids the problem of over polarization and is numerically stable and efficient for dynamics simulation. Finally, compared to the traditional fixed AMBER charge model, the present method only adds about 5% additional computational time and is therefore highly efficient for practical applications.

Computational studies of water exchange around aqueous Li+ with polarizable potential models

To enhance our understanding of the mechanism of water-exchange around aqueous Li(+), we carried out a systematic study on this system using molecular dynamics simulations with polarizable potential models. The mechanistic properties associated with the water-exchange process, such as potentials of mean force, time dependent transmission coefficients, and the corresponding rate constants, were examined using transition rate theory, the reactive flux method, and Grote-Hynes treatments of the dynamic response of the solvent. We compared the computed rate theory results with results from previous corresponding studies in which classical non-polarizable force fields were used. Our computed barrier heights for water exchange are significantly larger than those obtained using classical non-polarizable force fields. We also studied the effect of pressure on water-exchange rates and the corresponding activation volume. Our computed rate results for water exchange increase with pressure; therefore, a small negative activation volume is observed.

A polarizable force field of dipalmitoylphosphatidylcholine based on the classical Drude model for molecular dynamics simulations of lipids.

A polarizable force field of saturated phosphatidylcholine-containing lipids based on the classical Drude oscillator model is optimized and used in molecular dynamics simulations of bilayer and monolayer membranes. The hierarchical parametrization strategy involves the optimization of parameters for small molecules representative of lipid functional groups, followed by their application in larger model compounds and full lipids. The polar headgroup is based on molecular ions tetramethyl ammonium and dimethyl phosphate, the esterified glycerol backbone is based on methyl acetate, and the aliphatic lipid hydrocarbon tails are based on linear alkanes. Parameters, optimized to best represent a collection of gas and liquid properties for these compounds, are assembled into a complete model of dipalmitoylphosphatidylcholine (DPPC) lipids that is tested against the experimental properties of bilayer and monolayer membranes. The polarizable model yields average structural properties that are in broad accord with experimental data. The area per lipid of the model is 60 Å(2), slightly smaller than the experimental value of 63 Å(2). The order parameters from nuclear magnetic resonance deuterium quadrupolar splitting measures, the electron density profile, and the monolayer dipole potential are in reasonable agreement with experimental data, and with the nonpolarizable CHARMM C36 lipid force field.

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.

Molecular dynamics simulations of n-hexane at 1-butyl-3-methylimidazolium bis(trifluoromethylsulfonyl) imide interface

Molecular dynamics simulations of n-hexane adsorbed onto the interface of 1-butyl-3-methylimidazolium bis(trifluoromethylsulfonyl) imide ([bmim][Tf2N]) are performed at three n-hexane surface densities, ranged from 0.7 to 2.3 μmol/m(2) at 300 K. For [bmim][Tf2N] room-temperature ionic liquid, we use a non-polarizable all-atom force field with the partial atomic charges based on ab initio calculations for the isolated ion pair. The net charges of the ions are ±0.89e, which mimics the anion to cation charge transfer and polarization effects. The OPLS-AA force field is employed for modeling of n-hexane. The surface tension is computed using the mechanical route and its value decreases with increase of the n-hexane surface density. The [bmim][Tf2N]/n-hexane interface is analyzed using the intrinsic method, and the structural and dynamic properties of the interfacial, sub-interfacial, and central layers are computed. We determine the surface roughness, global and intrinsic density profiles, and orientation ordering of the molecules to describe the structure of the interface. We further compute the survival probability, normal and lateral self-diffusion coefficients, and re-orientation correlation functions to elucidate the effects of n-hexane on dynamics of the cations and anions in the layers.

QM/MM molecular dynamics simulations of the hydration of Mg(II) and Zn(II) ions

The hydration of Mg2+ and Zn2+ is examined using molecular dynamics simulations using 3 computational approaches of increasing complexity: the CHARMM nonpolarizable force field based on the TIP3P water model, the Drude polarizable force field based on the SWM4-NDP water model, and a combined QM/MM approach in which the inner coordination sphere is represented using a high-quality density functional theory (DFT) model (PBE/def2-TZVPP), and the remainder of the bulk water solvent is represented using the polarizable SWM4-NDP water model. The characteristic structural distribution functions (radial, angular, and tilt) are comparedand show very good agreement between the polarizable force field and QM/MM approaches. They predict an average Mg–O distance of 2.11 Å and an Zn–O distance of 2.13 Å, in good agreement with the available experimental neutron scattering and EXAFS data, while the Mg–O distances calculated using the nonpolarizable force field are 0.1 Å too short. Mg2+ (aq) and Zn2+ (aq) both have a coordination number of 6 and have a remarkably similar octahedral coordination mode, despite the chemical differences between these ions. Thermodynamic integration was used to calculate the relative hydration free energies of these ions (ΔΔGhydr). The nonpolarizable model is in error by 60 kcal mol– 1 and incorrectly predicts that Mg2+ has the more negative hydration energy. The Drude polarizable model predicts a ΔΔGhydr of only –13.2 kcal kcal mol– 1, an improvement over the results of the nonpolarizable force field, but still signficantly different than the experimental value of –30.1 kcal mol–1. The combined QM/MM approach performs much better, predicting a ΔΔGhydr of –34.8 kcal mol–1 in excellent agreement with experiment. These calculations support the experimental observation that Zn2+ has more favourable solvation free energy than Mg2+ despite having a very similar solvation structure.

Interchain hydrophobic clustering promotes rigidity in HIV-1 protease flap dynamics: new insights from Molecular Dynamics

The dynamics of HIV-1 protease (HIV-pr), a drug target for HIV infection, has been studied extensively by both computational and experimental methods. The flap dynamics of HIV-pr is considered to be more important for better ligand binding and enzymatic actions. Moreover, it has been demonstrated that the drug-induced mutations can change the flap dynamics of HIV-pr affecting the binding affinity of the ligands. Therefore, detailed understanding of flap dynamics is essential for designing better inhibitors. Previous computational investigations observed significant variation in the flap opening in nanosecond time scale indicating that the dynamics is highly sensitive to the simulation protocols. To understand the sensitivity of the flap dynamics on the force field and simulation protocol, molecular dynamics simulations of HIV-pr have been performed with two different AMBER force fields, ff99 and ff02. Two different trajectories (20 ns each) were obtained using the ff99 and ff02 force field. The results showed polarizable force field (ff02) make the flap tighter than the nonpolarizable force field (ff99). Some polar interactions and hydrogen bonds involving flap residues were found to be stronger with ff02 force field. The formation of interchain hydrophobic cluster (between flap tip of one chain and active site wall of another chain) was found to be dominant in the semi-open structures obtained from the simulations irrespective of the force field. It is proposed that an inhibitor, which will promote this interchain hydrophobic clustering, may make the flaps more rigid, and presumably the effect of mutation would be small on ligand binding.

Towards molecular dynamics simulations of chiral room-temperature ionic liquids

We developed an all-atom non-polarizable force field for simulations of two chiral room-temperature ionic liquids (RTILs) derived from 1-n-butyl-3-methylimidazolium bromide ([bmim][Br]); namely, (R)-1-butyl-3-(3-hydroxy-2-methylpropyl)imidazolium bromide (hydroxypropyl) and 1-butyl-3-[(1R)-nopyl]imidazolium bromide (nopyl). The force field adopts the CHARMM parameters for intramolecular and repulsion–dispersion interactions, and it employs reduced partial atomic charges of the ions which we derived by quantum-mechanical calculations. The atomic charges of a cation are first assigned by the ChelpG method to yield the integer net charge +e. The ChelpG charges of the cation are then uniformly scaled down by a factor equal to the absolute value of a charge on the bromine anion obtained by the QTAIM method for the cation–anion pair. The net charges of the ions are around ±0.8e, which mimics the anion to cation charge transfer and many-body effects. Molecular dynamics (MD) simulations are used to estimate the bulk density and heat of vaporization of the RTILs at 300K and 1bar. MD simulations in the slab geometry combined with the intrinsic interface analysis are then employed to compare the vapor–liquid interface of the chiral RTILs with the interface of [bmim][Br] at 300K. We focus on the structural properties of the interfaces by computing the density profiles and surface roughness. The chiral RTILs are currently tested for separation of racemic mixtures of α-pinenes and the simulations complement experimental characterization of the chiral RTILs.

Direct folding simulation of a long helix in explicit water

A recently proposed Polarizable Hydrogen Bond (PHB) method has been employed to simulate the folding of a 53 amino acid helix (PDB ID 2KHK) in explicit water. Under PHB simulation, starting from a fully extended structure, the peptide folds into the native state as confirmed by measured time evolutions of radius of gyration, root mean square deviation (RMSD), and native hydrogen bond. Free energy and cluster analysis show that the folded helix is thermally stable under the PHB model. Comparison of simulation results under, respectively, PHB and standard nonpolarizable force field demonstrates that polarization is critical for stable folding of this long α-helix.

Ion Pairing in Aqueous Lithium Salt Solutions with Monovalent and Divalent Counter-Anions

Molecular dynamics simulations of concentrated aqueous solutions of LiCl and Li2SO4 were conducted in order to provide molecular insight into recent neutron scattering data. The structures predicted from the molecular dynamics simulations using standard nonpolarizable force fields provided a very poor fit to the experiment; therefore, refinement was needed. The electronic polarizability of the medium was effectively accounted for by implementing the electronic continuum correction, which practically means rescaling the ionic charges. Consistent with previous studies, we found that this approach in each case provided a significantly improved fit to the experimental data, which was further enhanced by slightly adjusting the radius of the lithium ion. The polarization effect was particularly pronounced in the Li2SO4 solution where the ions in the nonpolarizable simulations tended to cluster unphysically. With the above alterations, the employed force field displayed an excellent fit to the neutron scattering data and provided a useful interpretative framework for the experimental measurements. At the same time, the present study underlines the importance of solvent polarization effects in hydration of ions with high charge density.

A Polarizable Force-Field for Cholesterol and Sphingomyelin.

A polarizable force-field, based on the Drude oscillator model, has been developed for cholesterol and the sphingomyelin class of lipids crucial to membrane raft formation, with testing performed on several 100 ns simulations. These have been validated against experimental observables as well as previous, nonpolarizable force-fields. Membrane bilayer properties, such as area-per-lipid and membrane thickness, produce results of comparable quantitative accuracy to those from the nonpolarizable force-field, while the membrane dipole potential is computed to be approximately 260 mV for a pure sphingomyelin bilayer, showing quantitative agreement with experimental results.

Predicting Mutation-Induced Stark Shifts in the Active Site of a Protein with a Polarized Force Field

The electric field inside a protein has a significant effect on the protein structure, function, and dynamics. Recent experimental developments have offered a direct approach to measure the electric field by utilizing a nitrile-containing inhibitor as a probe that can deliver a unique vibration to the specific site of interest in the protein. The observed frequency shift of the nitrile stretching vibration exhibits a linear dependence on the electric field at the nitrile site, thus providing a direct measurement of the relative electric field. In the present work, molecular dynamics simulations were carried out to compute the electric field shift in human aldose reductase (hALR2) using a polarized protein-specific charge (PPC) model derived from fragment-based quantum-chemistry calculations in implicit solvent. Calculated changes of electric field in the active site of hALR2 between the wild type and mutants were directly compared with measured vibrational frequency shifts (Stark shifts). Our study demonstrates that the Stark shifts calculated using the PPC model are in much better agreement with the experimental data than widely used nonpolarizable force fields, indicating that the electronic polarization effect is important for the accurate prediction of changes in the electric field inside proteins.

Thermodynamics of iodide adsorption at the instantaneous air-water interface

We performed molecular dynamics simulations using both polarizable and non-polarizable force fields to study the adsorption of iodide to the air-water interface. A novel aspect of our analysis is that the progress of ion adsorption is measured as the distance from the instantaneous interface, which is defined by a coarse-graining scheme proposed recently by Willard and Chandler ["Instantaneous liquid interfaces," J. Phys. Chem. B 114, 1954-1958 (2010)]. Referring structural and thermodynamic quantities to the instantaneous interface unmasks molecular-scale details that are obscured by thermal fluctuations when the same quantities are referred to an average measure of the position of the interface, such as the Gibbs dividing surface. Our results suggest that an ion adsorbed at the interface resides primarily in the topmost water layer, and the interfacial location of the ion is favored by enthalpy and opposed by entropy.

Multipole-Based Force Fields from ab Initio Interaction Energies and the Need for Jointly Refitting All Intermolecular Parameters.

Distributed atomic multipole (MTP) moments promise significant improvements over point charges (PCs) in molecular force fields, as they (a) more realistically reproduce the ab initio electrostatic potential (ESP) and (b) allow to capture anisotropic atomic properties such as lone pairs, conjugated systems, and σ holes. The present work focuses on the question of whether multipolar electrostatics instead of PCs in standard force fields leads to quantitative improvements over point charges in reproducing intermolecular interactions. To this end, the interaction energies of two model systems, benzonitrile (BZN) and formamide (FAM) homodimers, are characterized over a wide range of dimer conformations. It is found that although with MTPs the monomer ab initio ESP can be captured better by about an order of magnitude compared to point charges (PCs), this does not directly translate into better describing ab initio interaction energies compared to PCs. Neither ESP-fitted MTPs nor refitted Lennard-Jones (LJ) parameters alone demonstrate a clear superiority of atomic MTPs. We show that only if both electrostatic and LJ parameters are jointly optimized in standard, nonpolarizable force fields, atomic are MTPs clearly beneficial for reproducing ab initio dimerization energies. After an exhaustive exponent scan, we find that for both BZN and FAM, atomic MTPs and a 9-6 LJ potential can reproduce ab initio interaction energies with ∼30% (RMSD 0.13 vs 0.18 kcal/mol) less error than point charges (PCs) and a 12-6 LJ potential. We also find that the improvement due to using MTPs with a 9-6 LJ potential is considerably more pronounced than with a 12-6 LJ potential (≈ 10%; RMSD 0.19 versus 0.21 kcal/mol).

Exploring Adsorption of Water and Ions on Carbon Surfaces Using a Polarizable Force Field

Graphene, carbon nanotubes, and fullerenes are of great interest due to their unique properties and diverse applications in biology, molecular electronics, and materials science. Therefore, there is demand for methods that can accurately model the interface between carbon surfaces and their environment. In this letter we compare results for complexes of water, potassium ion, and chloride ion with graphene, carbon nanotube, and fullerene surfaces using a standard non-polarizable force field (OPLS-AA), a polarizable force field (OPLS-AAP), DFT, and ab initio theory. For interactions with water, OPLS-AA with the TIP3P or TIP4P water models describes the interactions with benzene (C(6)H(6)) and coronene (C(24)H(12)) well; however, for acenes larger than circumcoronene (C(54)H(18)) and especially for C(60), the interaction energies are somewhat too weak and polarization is needed. For ions interacting with carbon surfaces, inclusion of polarization is essential, and OPLS-AAP is found to perform well in comparison to the highest-level quantum mechanical methods. Overall, OPLS-AAP provides an accurate and computationally efficient force field for modeling condensed-phase systems featuring carbon surfaces.

Polarizable Site Charge Model at Liquid/Solid Interfaces for Describing Surface Polarity: Application to Structure and Molecular Dynamics of Water/Rutile TiO2(110) Interface

We present a novel scheme to construct a polarizable force field for liquid/solid interfaces, which takes into account the effect of the surface polarity induced by liquid-solid interactions explicitly. We extend the charge response kernel (CRK) method for molecules to solid surfaces by introducing the surface CRK. The CRK parameters are systematically determined by the first-principles calculations in the slab model with the dipole-correction method. Our methodology is applied to the water/clean rutile TiO2(110) interface. Structures and induced charges of a single water molecule attached to the TiO2 surface optimized by our polarizable force field show good agreement with those predicted by the first-principles calculations. Further, we carried out MD simulations for the liquid water/TiO2 interface and found three stable structures of water attached to the TiO2 surface. Two of them are predicted by both the polarizable and the nonpolarizable force fields, while the polarizable force field model predicts a structure of water with the hydrogen and oxygen atoms interacting with the oxygen atom of the surface TiO2 and the hydrogen atom of the other water molecule, respectively, which was reported by the previous first-principles MD simulation. This indicates that the dipole moments of water and TiO2 induced by the water-TiO2 interactions have significant impact on molecular conformations of the water/TiO2 interface.