Force field for copper clusters and nanoparticles

An atomic force field for simulating copper clusters and nanoparticles is developed. More than 2000 cluster configurations of varying size and shape are used to constrain the parametrization of the copper force field. Binding energies for these training clusters were computed using density functional theory. Extensive testing shows that the copper force field is fast and reliable for near‐equilibrium structures of clusters, ranging from only a few atoms to large nanoparticles that approach bulk structure. Nonequilibrium dissociation and compression structures that are included in the training set are also well described by the force field. Implications for molecular dynamics simulations and extensions to other metallic and covalent systems are discussed. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

Internal coordinate modeling of DNA: Force field comparisons

The Jumna internal coordinate program for modeling nucleic acids was extended to include the force field developed for the Amber program. This forms a bridge between internal and Cartesian coordinate modeling approaches. Using the extensive conformational mapping and substate search facilities available within Jumna, we rigorously compared the behavior of the different force fields and also of different continuum solvent models. The results, which help to explain trends seen in earlier minimization and molecular dynamics simulations, point to the superiority of the latest Amber parameterization (Parm94) and to a surprising degree of agreement with the Flex force field originally developed for Jumna. © 1997 John Wiley & Sons, Inc. J Comput Chem 18 :1043–1055, 1997     

RedMD—Reduced molecular dynamics package

We developed a software package (RedMD) to perform molecular dynamics simulations and normal mode analysis of reduced models of proteins, nucleic acids, and their complexes. With RedMD one can perform molecular dynamics simulations in a microcanonical ensemble, with Berendsen and Langevin thermostats, and with Brownian dynamics. We provide force field and topology generators which are based on the one‐bead per residue/nucleotide elastic network model and its extensions. The user can change the force field parameters with the command line options that are passed to generators. Also, the generators can be modified, for example, to add new potential energy functions. Normal mode analysis tool is available for elastic or anisotropic network models. The program is written in C and C++ languages and the structure/topology of a molecule is based on an XML format. OpenMP technology for shared‐memory architectures was used for code parallelization. The code is distributed under GNU public licence and available at http://bionano.icm.edu.pl/software/ . © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

Polarizable and flexible model for ethanol

A polarizable, flexible model for ethanol is obtained based on an extensive series of B3LYP/6-311++G(d,p) calculations and molecular dynamics simulations. The ethanol model includes electric-field dependence in both the atomic charges and the intramolecular degrees of freedom. Field-dependent intramolecular potentials have been attempted only once previously, for OH and HH stretches in water [P. Cicu et al., J. Chem. Phys. 112, 8267 (2000)]. The torsional potential involving the hydrogen-bonding hydrogen in ethanol is found to be particularly field sensitive. The methodology for developing field-dependent potentials can be readily generalized to other molecules and is discussed in detail. Molecular dynamics simulations of bulk ethanol are performed and the results are assessed based on comparisons with the self-diffusion coefficient [N. Karger et al., J. Chem. Phys. 93, 3437 (1990)], dielectric constant [J. T. Kindt and C. A. Schmuttenmaer, J. Phys. Chem. 100, 10373 (1996)], enthalpy of vaporization [R. C. Wilhoit and B. J. Zwolinski, J. Phys. Chem. Ref. Data, Suppl. 2, 2 (1973)], and experimental interatomic distributions [C. J. Benmore and Y. L. Loh, J. Chem. Phys. 112, 5877 (2000)]. The simultaneous variation of the atomic charges and the intramolecular potentials requires modified equations of motion and a multiple time step algorithm has been implemented to solve these equations. The article concludes with a discussion of the bulk structure and properties with an emphasis on the hydrogen bonding network.     

Molecular dynamics simulations of hexahydro-1,3,5-trinitro-1,3,5-s-triazine (RDX) using a combined Sorescu-Rice-Thompson AMBER force field

We present the results of molecular dynamics simulations of crystalline hexahydro-1,3,5-trinitro-1,3,5-s-triazine (RDX) using the SRT-AMBER force field (P. M. Agrawal et al., J. Phys. Chem. B 2006, 110, 5721), which combines the rigid-molecule force field developed by Sorescu-Rice-Thompson (D. C. Sorescu, B. M. Rice, and D. L. Thompson, J. Phys. Chem. B 1997, 101, 798) with the intramolecular interactions obtained from the Generalized AMBER Force Field (Wang et al., J. Comput. Chem. 2004, 25, 1157). The calculated crystal density at room conditions is about 10% lower than the measured value, while the lattice parameters and thermodynamic melting point are within about 5% at ambient pressure. The chair and inverted chair conformation, bond lengths, and bond angles of the RDX molecule are accurately predicted; however, there are some inaccuracies in the calculated orientations of the NO2 groups. The SRT-AMBER force field predicts overall reasonable results, but modifications, probably in the torsional parameters, are needed for a more accurate force field.     

A molecular mechanics model for imatinib and imatinib:kinase binding

Imatinib is an important anticancer drug, which binds specifically to the Abl kinase and blocks its signalling activity. To model imatinib:protein interactions, we have developed a molecular mechanics force field for imatinib and four close analogues, which is consistent with the CHARMM force field for proteins and nucleic acids. Atomic charges and Lennard‐Jones parameters were derived from a supermolecule ab initio approach. We considered the ab initio energies and geometries of a probe water molecule interacting with imatinib fragments at 32 different positions. We considered both a neutral and a protonated imatinib. The final RMS deviation between the ab initio and force field energies, averaged over both forms, was 0.2 kcal/mol. The model also reproduces the ab initio geometry and flexibility of imatinib. To apply the force field to imatinib:Abl simulations, it is also necessary to determine the most likely imatinib protonation state when it binds to Abl. This was done using molecular dynamics free energy simulations, where imatinib is reversibly protonated during a series of MD simulations, both in solution and in complex with Abl. The simulations indicate that imatinib binds to Abl in its protonated, positively‐charged form. To help test the force field and the protonation prediction, we did MD free energy simulations that compare the Abl binding affinities of two imatinib analogs, obtaining good agreement with experiment. Finally, two new imatinib variants were considered, one of which is predicted to have improved Abl binding. This variant could be of interest as a potential drug. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

A new soft‐core potential function for molecular dynamics applied to the prediction of protein loop conformations

We have developed a new soft‐core potential function for the conformational search of complex systems with molecular dynamics. The potential function was designed to maintain the main equilibrium properties of the original force field, which means that the soft‐core potential gives physically realistic performance also without additional restraints, different from most of the previous soft‐core potential functions. The performance of the method was demonstrated by applying it to the problem of finding native conformations for protein loops. Short loops from neocarzinostatin and parvalbumin were used as the first test cases. The use of the new soft‐core potential function was shown to improve significantly the performance of molecular dynamics in the search of the native conformation of protein loops. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 388–397, 2000     

CHARMM general force field: A force field for drug‐like molecules compatible with the CHARMM all‐atom additive biological force fields

The widely used CHARMM additive all‐atom force field includes parameters for proteins, nucleic acids, lipids, and carbohydrates. In the present article, an extension of the CHARMM force field to drug‐like molecules is presented. The resulting CHARMM General Force Field (CGenFF) covers a wide range of chemical groups present in biomolecules and drug‐like molecules, including a large number of heterocyclic scaffolds. The parametrization philosophy behind the force field focuses on quality at the expense of transferability, with the implementation concentrating on an extensible force field. Statistics related to the quality of the parametrization with a focus on experimental validation are presented. Additionally, the parametrization procedure, described fully in the present article in the context of the model systems, pyrrolidine, and 3‐phenoxymethylpyrrolidine will allow users to readily extend the force field to chemical groups that are not explicitly covered in the force field as well as add functional groups to and link together molecules already available in the force field. CGenFF thus makes it possible to perform “all‐CHARMM” simulations on drug‐target interactions thereby extending the utility of CHARMM force fields to medicinally relevant systems. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

On the accuracy of force fields for predicting the physical properties of dimethylnitramine

The accuracy of three force fields for predicting the physical properties of dimethylnitramine (DMNA) has been investigated by using molecular dynamics simulations. The Sorescu, Rice, and Thompson (SRT) (J. Phys. Chem. B 1997, 101, 798) rigid-molecule, flexible generalized AMBER (J. Comput. Chem. 2004, 25, 1157), and Smith et al. flexible force fields (J. Phys. Chem. B 1999, 103, 705) were tested. The density, lattice parameters, isotherm, and melting point of DMNA are calculated using classical molecular dynamics. Except for the melting point, the predictions of the three force fields are in reasonable agreement with experimental values. The calculated thermodynamic melting points (Tmp) for the SRT, AMBER, and Smith et al. force fields are 380, 360, and 260 K, respectively. The experimental value is 331 K. Modifications of the torsional barriers in the AMBER force field resulted in Tmp = 346 K, in good agreement with the experimental value of 331 K. The calculated lattice parameters and bulk modulus are also improved with the modifications of the AMBER potential. The results indicate that, although not sufficiently accurate without modifications, the general force fields such as AMBER provide the basis for developing force fields that correctly predict the physical properties of nitramines.     

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     

Force field parameters for the simulation of modified histone tails

We describe the development of force field parameters for methylated lysines and arginines, and acetylated lysine for the CHARMM all‐atom force field. We also describe a CHARMM united‐atom force field for modified sidechains suitable for use with fragment‐based docking methods. The development of these parameters is based on results of ab initio quantum mechanics calculations of model compounds with subsequent refinement and validation by molecular mechanics and molecular dynamics simulations. The united‐atom parameters are tested by fragment docking to target proteins using the MCSS procedure. The all‐atom force field is validated by molecular dynamics simulations of multiple experimental structures. In both sets of calculations, the computational predictions using the force field were compared to the corresponding experimental structures. We show that the parameters yield an accurate reproduction of experimental structures. Together with the existing CHARMM force field, these parameters will enable the general modeling of post‐translational modifications of histone tails. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

Alchemical prediction of hydration free energies for SAMPL

Hydration free energy calculations have become important tests of force fields. Alchemical free energy calculations based on molecular dynamics simulations provide a rigorous way to calculate these free energies for a particular force field, given sufficient sampling. Here, we report results of alchemical hydration free energy calculations for the set of small molecules comprising the 2011 Statistical Assessment of Modeling of Proteins and Ligands challenge. Our calculations are largely based on the Generalized Amber Force Field with several different charge models, and we achieved RMS errors in the 1.4-2.2 kcal/mol range depending on charge model, marginally higher than what we typically observed in previous studies (Mobley et al. in J Phys Chem B 111(9):2242-2254, 2007, J Chem Theory Comput 5(2):350-358, 2009, J Phys Chem B 115:1329-1332, 2011; Nicholls et al. in J Med Chem 51:769-779, 2008; Klimovich and Mobley in J Comput Aided Mol Design 24(4):307-316, 2010). The test set consists of ethane, biphenyl, and a dibenzyl dioxin, as well as a series of chlorinated derivatives of each. We found that, for this set, using high-quality partial charges from MP2/cc-PVTZ SCRF RESP fits provided marginally improved agreement with experiment over using AM1-BCC partial charges as we have more typically done, in keeping with our recent findings (Mobley et al. in J Phys Chem B 115:1329-1332, 2011). Switching to OPLS Lennard-Jones parameters with AM1-BCC charges also improves agreement with experiment. We also find a number of chemical trends within each molecular series which we can explain, but there are also some surprises, including some that are captured by the calculations and some that are not.     

Improved treatment of cyclic β‐amino acids and successful prediction of β‐peptide secondary structure using a modified force field: AMBER*C

We added parameters to the AMBER* force field to model cyclic β‐amino acid derivatives more accurately within the commonly used MacroModel program. In an effort to generate an improved treatment of cyclohexane and cyclopentane conformational preferences, carbon–carbon torsional parameters were modified and incorporated into a force field we call AMBER*C. Simulation of trans‐2‐aminocyclohexanecarboxylic acid (trans‐ACHC) and trans‐2‐aminocyclopentanecarboxylic acid (trans‐ACPC) derivatives using AMBER*C produces more realistic energy differences between (pseudo)diaxial and (pseudo)diequatorial conformations than does simulation using AMBER*. AMBER*C molecular dynamics simulations more accurately reproduce the experimental hydrogen‐bonding tendencies of simple diamide derivatives of trans‐ACHC and trans‐ACPC than do simulations using the AMBER* force field. More importantly, this modified force field allows accurate qualitative prediction of the helical secondary structures adopted by β‐amino acid homo‐oligomers. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 763–773, 2000     

Using one‐step perturbation to predict the effect of changing force‐field parameters on the simulated folding equilibrium of a β‐peptide in solution

Computer simulation using molecular dynamics is increasingly used to simulate the folding equilibria of peptides and small proteins. Yet, the quality of the obtained results depends largely on the quality of the force field used. This comprises the solute as well as the solvent model and their energetic and entropic compatibility. It is, however, computational very expensive to perform test simulations for each combination of force‐field parameters. Here, we use the one‐step perturbation technique to predict the change of the free enthalpy of folding of a β‐peptide in methanol solution due to changing a variety of force‐field parameters. The results show that changing the solute backbone partial charges affects the folding equilibrium, whereas this is relatively insensitive to changes in the force constants of the torsional energy terms of the force field. Extending the cut‐off distance for nonbonded interactions beyond 1.4 nm does not affect the folding equilibrium. The same result is found for a change of the reaction‐field permittivity for methanol from 17.7 to 30. The results are not sensitive to the criterion, e.g., atom‐positional RMSD or number of hydrogen bonds, that is used to distinguish folded and unfolded conformations. Control simulations with perturbed Hamiltonians followed by backward one‐step perturbation indicated that quite large perturbations still yield reliable results. Yet, perturbing all solvent molecules showed where the limitations of the one‐step perturbation technique are met. The evaluated methodology constitutes an efficient tool in force‐field development for molecular simulation by reducing the number of required separate simulations by orders of magnitude. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

Molecular dynamics simulations and free energy calculations on the enzyme 4‐hydroxyphenylpyruvate dioxygenase

4‐Hydroxyphenylpyruvate dioxygenase is a relevant target in both pharmaceutical and agricultural research. We report on molecular dynamics simulations and free energy calculations on this enzyme, in complex with 12 inhibitors for which experimental affinities were determined. We applied the thermodynamic integration approach and the more efficient one‐step perturbation. Even though simulations seem well converged and both methods show excellent agreement between them, the correlation with the experimental values remains poor. We investigate the effect of slight modifications on the charge distribution of these highly conjugated systems and find that accurate models can be obtained when using improved force field parameters. This study gives insight into the applicability of free energy methods and current limitations in force field parameterization. © 2011 Wiley Periodicals, Inc. J Comput Chem 2011     

Molecular dynamics simulations of surface-initiated melting of nitromethane

The melting of nitromethane initiated at solid-vacuum interfaces has been investigated using molecular dynamics nvt simulations with a realistic force field [D. C. Sorescu et al., J. Phys. Chem. B 104, 8406 (2000)]. The calculated melting point (251+/-5 K) is in good agreement with experiment (244.73 K) and values obtained previously (approximately 255.5 and 266.5+/-8 K) using other simulation methods [P. M. Agrawal et al., J. Chem. Phys. 119, 9617 (2003)]. Analyses of the molecular orientations and diffusion during the simulations as functions of the distance from the exposed surfaces show that the melting is a direct crystal-to-liquid transition, in which the molecules first gain rotational freedom, then mobility. There is a slight dependence of the melting temperature on the exposed crystallographic face.     

The QM‐MM interface for CHARMM‐deMon

We present a new QM/MM interface for fast and efficient simulations of organic and biological molecules. The CHARMM/deMon interface has been developed and tested to perform minimization and atomistic simulations for multi‐particle systems. The current features of this QM/MM interface include readability for molecular dynamics, tested compatibility with Free Energy Perturbation simulations (FEP) using the dual topology/single coordinate method. The current coupling scheme uses link atoms, but further extensions of the code to incorporate other available schemes are planned. We report the performance of different levels of theory for the treatment of the QM region, while the MM region was represented by a classical force‐field (CHARMM27) or a polarizable force‐field based on a simple Drude model. The current QM/MM implementation can be coupled to the dual‐thermostat method and the VV2 integrator to run molecular dynamics simulations. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Hybrid Monte Carlo with adaptive temperature in mixed-canonical ensemble: Efficient conformational analysis of RNA

A hybrid Monte Carlo method with adaptive temperature choice is presented that exactly generates the distribution of a mixed-canonical ensemble composed of two canonical ensembles at low and high temperature. The analysis of resulting Markov chains with the reweighting technique shows an efficient sampling of the canonical distribution at low temperature whereas the high temperature component facilitates conformational transitions, which allows shorter simulation times. The algorithm is tested by comparing analytical and numerical results for the small n-butane molecule before simulations are performed for a triribonucleotide. Sampling the complex multiminima energy landscape of this small RNA segment, we observe enforced crossing of energy barriers. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1689–1697, 1998     

Reversible peptide folding: Dependence on molecular force field used

Temperature‐dependent nuclear magnetic resonance (NMR) and CD spectra of methanol solutions of a β‐heptapeptide have been interpreted in such a way that the secondary structure, a 3₁₄‐helix, is assumed to be stable in a temperature range of between 298 and 393 K. This is in contrast to the results of a 50‐ns molecular dynamics simulation using the GROMOS 96 force field, which found a melting temperature of about 340 K. This discrepancy is addressed by further computational studies using the OPLS‐AA force field. The conformational energetics of N‐formyl‐3‐aminobutanamide in vacuo are obtained using ab initio and density functional quantum‐mechanical calculations at the HF/6‐31G*, B3LYP/6‐31G*, and B3LYP/6‐311+G* levels of theory. The results permit development of torsional parameters for the OPLS‐AA force field that reproduce the conformational energetics of the monomer. By varying the development procedure, three parameter sets are obtained that focus on reproducing either low‐energy or high‐energy conformations. These parameter sets are tested by simulating the reversible folding of the β‐heptapeptide in methanol. The melting temperature of the helix formed (>360 K) is found to be higher than the one obtained from simulations using the GROMOS 96 force field (∼340 K). Differences in the potential energy functions of the latter two force fields are evaluated and point to the origins of the difference in stability. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 774–787, 2000