Symmetrization of the AMBER and CHARMM force fields

The AMBER and CHARMM force fields are analyzed from the viewpoint of the permutational symmetry of the potential for feasible exchanges of identical atoms and chemical groups in amino and nucleic acids. In each case, we propose schemes for symmetrizing the potentials, which greatly facilitate the bookkeeping associated with constructing kinetic transition networks via geometry optimization. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

Hierarchical atom type definitions and extensible all‐atom force fields

The extensibility of force field is a key to solve the missing parameter problem commonly found in force field applications. The extensibility of conventional force fields is traditionally managed in the parameterization procedure, which becomes impractical as the coverage of the force field increases above a threshold. A hierarchical atom‐type definition (HAD) scheme is proposed to make extensible atom type definitions, which ensures that the force field developed based on the definitions are extensible. To demonstrate how HAD works and to prepare a foundation for future developments, two general force fields based on AMBER and DFF functional forms are parameterized for common organic molecules. The force field parameters are derived from the same set of quantum mechanical data and experimental liquid data using an automated parameterization tool, and validated by calculating molecular and liquid properties. The hydration free energies are calculated successfully by introducing a polarization scaling factor to the dispersion term between the solvent and solute molecules. © 2015 Wiley Periodicals, Inc.     

ForceFit: A code to fit classical force fields to quantum mechanical potential energy surfaces

The ForceFit program package has been developed for fitting classical force field parameters based upon a force matching algorithm to quantum mechanical gradients of configurations that span the potential energy surface of the system. The program, which runs under UNIX and is written in C++, is an easy‐to‐use, nonproprietary platform that enables gradient fitting of a wide variety of functional force field forms to quantum mechanical information obtained from an array of common electronic structure codes. All aspects of the fitting process are run from a graphical user interface, from the parsing of quantum mechanical data, assembling of a potential energy surface database, setting the force field, and variables to be optimized, choosing a molecular mechanics code for comparison to the reference data, and finally, the initiation of a least squares minimization algorithm. Furthermore, the code is based on a modular templated code design that enables the facile addition of new functionality to the program. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

Carbohydrate solution simulations: producing a force field with experimentally consistent primary alcohol rotational frequencies and populations

We present a CHARMM Carbohydrate Solution Force Field (CSFF) suitable for nanosecond molecular dynamics computer simulations. The force field was derived from a recently published sugar parameter set.1 Dihedral angle parameters for the primary alcohol as well as the secondary hydroxyl groups were adjusted. Free energy profiles of the hydroxymethyl group for two monosaccharides (beta-D-glucose and beta-D-galactose) were calculated using the new parameter set and compared with similar force fields. Equilibrium rotamer populations obtained from the CSFF are in excellent agreement with NMR data (glucose gg:gt:tg approximately 66:33:1 and galactose gg:gt:tg approximately 4:75:21). In addition, the primary alcohol rotational frequency is on the nanosecond time scale, which conforms to experimental observations. Equilibrium population distributions of the primary alcohol conformers for glucose and galactose are reached within 10 nanoseconds of molecular dynamics simulations. In addition, gas phase vibrational frequencies computed for beta-D-glucose using this force field compare well with experimental frequencies. Carbohydrate parameter sets that produce both conformational energies and rotational frequencies for the pyranose primary alcohol group that are consistent with experimental observations should allow for increased accuracy in modeling the flexibility of biologically important (1-6)-linked saccharides in solution.     

A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS force-field parameter sets 53A5 and 53A6

Successive parameterizations of the GROMOS force field have been used successfully to simulate biomolecular systems over a long period of time. The continuing expansion of computational power with time makes it possible to compute ever more properties for an increasing variety of molecular systems with greater precision. This has led to recurrent parameterizations of the GROMOS force field all aimed at achieving better agreement with experimental data. Here we report the results of the latest, extensive reparameterization of the GROMOS force field. In contrast to the parameterization of other biomolecular force fields, this parameterization of the GROMOS force field is based primarily on reproducing the free enthalpies of hydration and apolar solvation for a range of compounds. This approach was chosen because the relative free enthalpy of solvation between polar and apolar environments is a key property in many biomolecular processes of interest, such as protein folding, biomolecular association, membrane formation, and transport over membranes. The newest parameter sets, 53A5 and 53A6, were optimized by first fitting to reproduce the thermodynamic properties of pure liquids of a range of small polar molecules and the solvation free enthalpies of amino acid analogs in cyclohexane (53A5). The partial charges were then adjusted to reproduce the hydration free enthalpies in water (53A6). Both parameter sets are fully documented, and the differences between these and previous parameter sets are discussed.     

Revised CHARMM force field parameters for iron‐containing cofactors of photosystem II

Photosystem II is a complex protein–cofactor machinery that splits water molecules into molecular oxygen, protons, and electrons. All‐atom molecular dynamics simulations have the potential to contribute to our general understanding of how photosystem II works. To perform reliable all‐atom simulations, we need accurate force field parameters for the cofactor molecules. We present here CHARMM bonded and non‐bonded parameters for the iron‐containing cofactors of photosystem II that include a six‐coordinated heme moiety coordinated by two histidine groups, and a non‐heme iron complex coordinated by bicarbonate and four histidines. The force field parameters presented here give water interaction energies and geometries in good agreement with the quantum mechanical target data. © 2017 Wiley Periodicals, Inc.     

Detailed considerations for a balanced and broadly applicable force field: a study of substituted benzenes modeled with OPLS-AA

Modern classical force fields have been traditionally parameterized by attempting to maximize agreement to any number of experimental and/or quantum mechanical target properties. As these force fields are pushed towards obtaining quantitative estimates of often subtle energetic differences, stringent and consistent parameterization criteria, particularly in regard to charge distributions, are required to ensure that systematic errors cancel, that parameters are transferable between molecules, and that performance does not significantly deteriorate when using more approximate methods, such as with continuum solvent models. Relative free energies of hydration are presented here for 40 mono- and disubstituted benzenes modeled with the OPLS-AA force field; heats of vaporization and pure liquid densities at standard conditions are presented when experimental data is available. Overall agreement between OPLS-AA and experiment is remarkable (average error = 0.5 kcal/mol for DeltaDeltaG(hydration), 1.0 kcal/mol for DeltaH(vap) (0), 0.02 g/mL for densities), yet several functional groups are identified as having consistent and correctable errors (alkyl-, nitro-, and thiobenzenes). Relative free energies of hydration obtained with rigorous free energy perturbations using explicit solvent are also compared with energies from minimizations using a generalized Born model (GB). There is high correlation between these estimates (R = 0.99), and as demonstrated here, reparameterization of the aforementioned groups can be guided with rapid GB calculations.     

A comparison of methods to compute the potential of mean force.

Most processes occurring in a system are determined by the relative free energy between two or more states because the free energy is a measure of the probability of finding the system in a given state. When the two states of interest are connected by a pathway, usually called reaction coordinate, along which the free-energy profile is determined, this profile or potential of mean force (PMF) will also yield the relative free energy of the two states. Twelve different methods to compute a PMF are reviewed and compared, with regard to their precision, for a system consisting of a pair of methane molecules in aqueous solution. We analyze all combinations of the type of sampling (unbiased, umbrella-biased or constraint-biased), how to compute free energies (from density of states or force averaging) and the type of coordinate system (internal or Cartesian) used for the PMF degree of freedom. The method of choice is constraint-bias simulation combined with force averaging for either an internal or a Cartesian PMF degree of freedom.     

Calculation of the water-cyclohexane transfer free energies of neutral amino acid side-chain analogs using the OPLS all-atom force field

We calculated the free energy of solvation of the neutral analogs of 18 amino acid side-chains (not including glycine and proline) using the OPLS all-atom force field in TIP4P water, SPC water, and cyclohexane by molecular dynamics simulation and thermodynamic integration. The average unsigned errors in the free energies of solvation in TIP4P, SPC, and cyclohexane are 4.4, 4.9, and 2.1 kJ/mol respectively. Most of the calculated hydration free energies are not favorable enough compared to experiment. The largest errors are found for tryptophan, histidine, glutamic acid, and glutamine. The average unsigned errors in the free energy of transfer from TIP4P to cyclohexane and from SPC to cyclohexane are 4.0 and 4.1 kJ/mol, respectively. The largest errors, of more than 7.5 kJ/mol, are found for histidine, glutamine, and glutamatic acid.     

Modeling the effect of ion-induced shock waves and DNA breakage with the reactive CHARMM force field

Ion-induced DNA damage is an important effect underlying ion beam cancer therapy. This article introduces the methodology of modeling DNA damage induced by a shock wave caused by a projectile ion. Specifically it is demonstrated how single- and double strand breaks in a DNA molecule could be described by the reactive CHARMM (rCHARMM) force field implemented in the program MBN Explorer. The entire workflow of performing the shock wave simulations, including obtaining the crucial simulation parameters, is described in seven steps. Two exemplary analyses are provided for a case study simulation serving to: (a) quantify the shock wave propagation and (b) describe the dynamics of formation of DNA breaks. The article concludes by discussing the computational cost of the simulations and revealing the possible maximal computational time for different simulation set-ups.     

Structure and Thermodynamics of Mg:Phosphate Interactions in Water: A Simulation Study

The association of Mg²⁺ and H₂PO₄^? in water can give insights into Mg:phosphate interactions in general, which are very widespread, but for which experimental data is surprisingly sparse. It is studied through molecular dynamics simulations (>100 ns) by using the polarizable AMOEBA force field, and the association free energy is computed for the first time. Explicit consideration of outer‐sphere and two types of inner‐sphere association provides considerable insight into the dynamics and thermodynamics of ion pairing. After careful assessment of the computational approximations, the agreement with experimental values indicates that the methodology can be extended to other inorganic and biological Mg:phosphate interactions in solution.     

Protegrin‐1 orientation and physicochemical properties in membrane bilayers studied by potential of mean force calculations

Protegrin‐1 (PG‐1) belongs to the family of antimicrobial peptides. It interacts specifically with the membrane of a pathogen and kills the pathogen by releasing its cellular contents. To fully understand the energetics governing the orientation of PG‐1 in different membrane environments and its effects on the physicochemical properties of the peptide and membrane bilayers, we have performed the potential of mean force (PMF) calculations as a function of its tilt angle at four distinct rotation angles in explicit membranes composed of either DLPC (1,2‐dilauroylphosphatidylcholine) or POPC (1‐palmitoyl‐2‐oleoylphosphatidylcholine) lipid molecules. The resulting PMFs in explicit lipid bilayers were then used to search for the optimal hydrophobic thickness of the EEF1/IMM1 implicit membrane model in which a two‐dimensional PMF in the tilt and rotation space was calculated. The PMFs in explicit membrane systems clearly reveal that the energetically favorable tilt angle is affected by both the membrane hydrophobic thickness and the PG‐1 rotation angle. Local thinning of the membrane around PG‐1 is observed upon PG‐1 tilting. The thinning is caused by both hydrophobic mismatch and arginine‐lipid head group interactions. The two‐dimensional PMF in the implicit membrane is in good accordance with those from the explicit membrane simulations. The ensemble‐averaged Val16 ¹⁵N and ¹³CO chemical shifts weighted by the two‐dimensional PMF agree fairly well with the experimental values, suggesting the importance of peptide dynamics in calculating such ensemble properties for direct comparison with experimental observables. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

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.     

Evaluating nonpolarizable nucleic acid force fields: A systematic comparison of the nucleobases hydration free energies and chloroform‐to‐water partition coefficients

Nucleic acid force fields have been shown to reproduce structural properties of DNA and RNA very well, but comparative studies with respect to thermodynamic properties are rare. As a test for thermodynamic properties, we have computed hydration free energies and chloroform‐to‐water partition coefficients of nucleobases using the AMBER‐99, AMBER‐gaff, CHARMM‐27, GROMOS‐45a4/53a6 and OPLS‐AA force fields. A mutual force field comparison showed a very large spread in the calculated thermodynamic properties, demonstrating that some of the parameter sets require further optimization. The choice of solvent model used in the simulation does not have a significant effect on the results. Comparing the hydration free energies obtained by the various force fields to the adenine and thymine experimental values showed a very large deviation for the GROMOS and AMBER parameter sets. Validation against experimental partition coefficients showed good agreement for the CHARMM‐27 parameter set. In view of mutation studies, differences in partition coefficient between two bases were also compared, and good agreement between experiments and calculations was found for the AMBER‐99 parameter set. Overall, the CHARMM‐27 parameter set performs best with respect to the thermodynamic properties tested here. © 2012 Wiley Periodicals, Inc.     

Improved precision and efficiency of free energy calculations for small systems using lambda-scaled atomic masses and separating conformational and transformational sampling

We present results showing the importance of appropriate treatment of atomic masses in molecular dynamics (MD)-based single topology free-energy perturbations (FEPs) on small molecule systems. The reversibility of gas phase simulations is significantly improved by scaling the atomic mass of mutated atoms with the lambda variable normally used for the scaling of energy terms. Because this effect is less pronounced for solvated systems, it will not cancel in estimates of the relative hydration free energy difference. The advantage of mass scaling is demonstrated by a null mutation of ethane to ethane and the calculation of the relative hydration free energy difference between ethane and n-propane. Furthermore, it is found that the simulation time necessary for converged MD/FEPs is prohibitively large for relative hydration free energy calculations on cyclic alkanes. Therefore, we explore an alternative free energy pathway including strongly constrained conformations to improve convergence in FEP simulations of flexible molecules.     

Peptide free‐energy profile is strongly dependent on the force field: Comparison of C96 and AMBER95

The C96 and AMBER95 force fields were compared with small model peptides Ac‐(Ala)_n‐NMe (Ac = CH₃CO, NMe = NHCH₃, n=2 and 3) in vacuo and in TIP3P water by computing the free‐energy profiles using multicanonical molecular dynamics method. The C96 force field is a modified version of the AMBER95 force field, which was adjusted to reproduce the energy difference between extended β‐ and constrained α‐helical energies for the alanine tetrapeptide, obtained by the high level ab initio MO method. The slight modification resulted in a large difference in the free energy profiles. The C96 force field prefers relatively extended conformers, whereas the AMBER95 force field favors turn conformations. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 748–762, 2000