Molecular mechanics force-field parameterization procedures

A set of procedures and guidelines are presented for the estimation of bond length, bond angle, and torsional potential constants for molecular mechanics force fields. The force field constants are ultimately derived by “subtracting” nonbonded molecular mechanics energies from corresponding molecular orbital energies using a model compound containing the chemical structure to be parameterized. Case study examples of bond length, bond angle, and torsional rotation force field parameterizations are presented. A general discussion of molecular mechanics force field parameterization strategy is included for reference and completeness. Finally, a curve-fitting program to generate force field parameters from raw data is given in Appendix I.     

A CNDO/2 study of the OH torsion in phenol and its hydrogen bonded complex with pyridine

The CNDO/2 molecular orbital method has been applied to the study of the OH torsion, in phenol and phenol—pyridine hydrogen bonded complex. The calculated torsional barrier (13.58 kJ mol^?1) and force constant (5.4 × 10^?20 J rad^?2) of phenol agree well with the experimental quantities. The calculated force constant of the corresponding vibration in phenol—pyridine is increased sixfold, reproducing closely the rise in the torsional frequency observed when phenol is complexed to strong acceptors. It is shown that according to CNDO theory, most of the increase can be attributed to the influence of the intermolecular force field and not to a major change in the torsional force constant.     

Computation of pressure components due to Class II force fields

It has been shown that purely angle dependent terms in a force field do not contribute to the total pressure in a molecular simulation. However, this is not the case for the individual components of the pressure tensor, and is also untrue for crossterms in the force field including bond stretch. In this article, we show that virial contributions to the pressure tensor are easily computed in terms of bond distance vectors and atom forces for the bond topologies present in a Class II force field. Results from a simulation of a phospholipid biomembrane using the cff97 force field show that angle and torsional crossterms make a significant contribution to the pressure tensor.     

Ab initio calculations and molecular mechanics (MM3) force field development for ammonium and protonated aliphatic amines

The geometries and vibrational frequencies of 11 training molecules containing the ammonium ion moiety were calculated at the MP2/6-31+G* level of theory. Various torsional energy profiles were also calculated using this basis set. From those ab initio calculations, a molecular mechanics (MM3) force field was developed using our Parameter Analysis and Refinement Toolkit System (PARTS). Using this set of parameters, the MM3 force field was found to well reproduce the molecular geometries and vibrational spectra for the all training molecules. CPU time was reduced from days to seconds. The availability of this new force field dramatically increases the feasibility of the computer-assisted drug design involving ammonium and protonated amino groups. © 1997 by John Wiley & Sons, Inc. J Comput Chem 18 : 1371–1391, 1997     

An improved OPLS-AA force field for carbohydrates

This work describes an improved version of the original OPLS-all atom (OPLS-AA) force field for carbohydrates (Damm et al., J Comp Chem 1997, 18, 1955). The improvement is achieved by applying additional scaling factors for the electrostatic interactions between 1,5- and 1,6-interactions. This new model is tested first for improving the conformational energetics of 1,2-ethanediol, the smallest polyol. With a 1,5-scaling factor of 1.25 the force field calculated relative energies are in excellent agreement with the ab initio-derived data. Applying the new 1,5-scaling makes it also necessary to use a 1,6-scaling factor for the interactions between the C4 and C6 atoms in hexopyranoses. After torsional parameter fitting, this improves the conformational energetics in comparison to the OPLS-AA force field. The set of hexopyranoses included in the torsional parameter derivation consists of the two anomers of D-glucose, D-mannose, and D-galactose, as well as of the methyl-pyranosides of D-glucose, D-mannose. Rotational profiles for the rotation of the exocyclic group and of different hydroxyl groups are also compared for the two force fields and at the ab initio level of theory. The new force field reduces the overly high barriers calculated using the OPLS-AA force field. This leads to better sampling, which was shown to produce more realistic conformational behavior for hexopyranoses in liquid simulation. From 10-ns molecular dynamics (MD) simulations of alpha-D-glucose and alpha-D-galactose the ratios for the three different conformations of the hydroxymethylene group and the average (3)J(H,H) coupling constants are derived and compared to experimental values. The results obtained for OPLS-AA-SEI force field are in good agreement with experiment whereas the properties derived for the OPLS-AA force field suffer from sampling problems. The undertaken investigations show that the newly derived OPLS-AA-SEI force field will allow simulating larger carbohydrates or polysaccharides with improved sampling of the hydroxyl groups.     

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.     

Force-field parametrization and molecular dynamics simulations of Congo red

Congo red, a diazo dye widely used in medical diagnosis, is known to form supramolecular systems in solution. Such a supramolecular system may interact with various proteins. In order to examine the nature of such complexes empirical force field parameters for the Congo red molecule were developed. The parametrization of bonding terms closely followed the methodology used in the development of the charmm22 force field, except for the calculation of charges. Point charges were calculated from a fit to a quantum mechanically derived electrostatic potential using the CHELP-BOW method. Obtained parameters were tested in a series of molecular dynamics simulations of both a single molecule and a micelle composed of Congo red molecules. It is shown that newly developed parameters define a stable minimum on the hypersurface of the potential energy and crystal and ab initio geometries and rotational barriers are well reproduced. Furthermore, rotations around C-N bonds are similar to torsional vibrations observed in crystals of diphenyl-diazene, which confirms that the flexibility of the molecule is correct. Comparison of results obtained from micelles molecular dynamics simulations with experimental data shows that the thermal dependence of micelle creation is well reproduced.     

Force field parametrization and molecular dynamics simulation of flexible POSS-linked (NHC; phosphine) Ru catalytic complexes

In recent years, N-heterocyclic carbene (NHC) or phospine groups have been put forward as candidate catalysts ligands for olefin metathesis reactions to be performed using multistep methods. Some of these proposed ligands contain polyhedral oligomeric silsesquioxane (POSS) structures linked to NHC rings by means of alkyl chains. Some important properties for the prediction of catalytic activity, such as the theoretically defined buried volume, are related to the conformational characteristics of these complex ligands that can be studied through molecular dynamics simulations. However, the chemical structure of resulting catalytic complexes usually contains atoms or groups that are not included in the common forcefields used in simulations. In this work we focus on complexes formed by a catalytic metal center (Ru) with both phospine and POSS-linked NHC groups. The central part of the complexes contain atoms and groups that have bonds, bond angles, and torsional angles whose parameters have not been previously evaluated and included in existing force fields. We have performed basic ab initio quantum mechanical calculations based on the density functional theory to obtain energies for this central section. The force field parameters for bonds, bond angles, and torsional angles are then calculated from an analysis of energies calculated for the equilibrium and different locally deformed structures. Nonbonded interactions are also conveniently evaluated. From subsequent molecular dynamics simulations, we have obtained results that illustrate the conformational characteristics most closely connected with the catalytic activity.     

1,4-Interactions in molecular mechanics calculations of Ethers

Van der Waals and electrostatic interactions are found to be insufficient for the calculation of conformational energies of ethers by molecular mechanics. Low order torsional potential functions must be added for the potential about C-O bonds. A onefold term necessary for the CCOC-fragment is interpreted to be a substitute for gauche interactions present in CCCC-, but missing in CCOC-fragments. For the COCO fragment the anomeric effect must be included explicitly as another torsional energy term, but no such term is required to stabilize the gauche conformation for OCCO. With the resulting ether force field the geometries and energies of model compounds, many of them 1,3-dioxanes, are calculated with good accuracy.     

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.     

Controlling the shape and flexibility of arylamides: a combined ab initio, ab initio molecular dynamics, and classical molecular dynamics study

Using quantum chemistry plus ab initio molecular dynamics and classical molecular dynamics methods, we address the relationship between molecular conformation and the biomedical function of arylamide polymers. Specifically, we have developed new torsional parameters for a class of these polymers and applied them in a study of the interaction between a representative arylamide and one of its biomedical targets, the anticoagulant drug heparin. Our main finding is that the torsional barrier of a C(aromatic)-C(carbonyl) bond increases significantly upon addition of an o-OCH2CH2NH3+ substituent on the benzene ring. Our molecular dynamics studies that are based on the original general AMBER force field (GAFF) and GAFF modified to include our newly developed torsional parameters show that the binding mechanism between the arylamide and heparin is very sensitive to the choice of torsional potentials. Ab initio molecular dynamics simulation of the arylamide independently confirms the degree of flexibility we obtain by classical molecular dynamics when newly developed torsional potentials are used.     

Derivation of class II force fields. VIII. Derivation of a general quantum mechanical force field for organic compounds

A class II valence force field covering a broad range of organic molecules has been derived employing ab initio quantum mechanical "observables." The procedure includes selecting representative molecules and molecular structures, and systematically sampling their energy surfaces as described by energies and energy first and second derivatives with respect to molecular deformations. In this article the procedure for fitting the force field parameters to these energies and energy derivatives is briefly reviewed. The application of the methodology to the derivation of a class II quantum mechanical force field (QMFF) for 32 organic functional groups is then described. A training set of 400 molecules spanning the 32 functional groups was used to parameterize the force field. The molecular families comprising the functional groups and, within each family, the torsional angles used to sample different conformers, are described. The number of stationary points (equilibria and transition states) for these molecules is given for each functional group. This set contains 1324 stationary structures, with 718 minimum energy structures and 606 transition states. The quality of the fit to the quantum data is gauged based on the deviations between the ab initio and force field energies and energy derivatives. The accuracy with which the QMFF reproduces the ab initio molecular bond lengths, bond angles, torsional angles, vibrational frequencies, and conformational energies is then given for each functional group. Consistently good accuracy is found for these computed properties for the various types of molecules. This demonstrates that the methodology is broadly applicable for the derivation of force field parameters across widely differing types of molecular structures. Copyright 2001 John Wiley & Sons, Inc. J Comput Chem 22: 1782-1800, 2001     

A molecular dynamics simulation of the melting points and glass transition temperatures of myo- and neo-inositol

The heat of sublimation, density, melting point, and glass transition temperature are calculated for myo- and neo-inositol, using the condensed-phase optimized molecular potentials for atomistic simulation studies (COMPASS) force field and molecular dynamics techniques. Our results show that the calculated heats of sublimation and density are very close to the experimental values for both compounds. Furthermore, our simulated melting temperatures for myo- and neo-inositol also compare very well to the experimentally obtained data. The glass transition temperatures for myo- and neo-inositol have been calculated to be ca. 494 K and ca. 518 K, respectively, and the shape of the volume versus temperature plots produced are typical for a glass transition. As a result, it is our view that the COMPASS force field suitably describes these two compounds in molecular simulations and that molecular dynamics techniques, combined with this force field, can be used to simulate the melt and glass transitions for such molecules.     

A force field for allenes and for nonlinear acetylenes within the MM2 approximation

Geometries and heats of formation for allenes may now be calculated with an extended version of MM2. The torsional barriers about bonds in an acetylenic unit are much more complicated than originally thought, if the unit is nonlinear. Ab initio calculations have been used to show just what these torsional barriers are, and the MM2 force field has been adapted to this information.     

Balancing solvation and intramolecular interactions: toward a consistent generalized Born force field

The efficient and accurate characterization of solvent effects is a key element in the theoretical and computational study of biological problems. Implicit solvent models, particularly generalized Born (GB) continuum electrostatics, have emerged as an attractive tool to study the structure and dynamics of biomolecules in various environments. Despite recent advances in this methodology, there remain limitations in the parametrization of many of these models. In the present work, we demonstrate that it is possible to achieve a balanced implicit solvent force field by further optimizing the input atomic radii in combination with adjusting the protein backbone torsional energetics. This parameter optimization is guided by the potentials of mean force (PMFs) between amino acid polar groups, calculated from explicit solvent free energy simulations, and by conformational equilibria of short peptides, obtained from extensive folding and unfolding replica exchange molecular dynamics (REX-MD) simulations. Through the application of this protocol, the delicate balance between the competing solvation forces and intramolecular forces appears to be better captured, and correct conformational equilibria for a range of both helical and beta-hairpin peptides are obtained. The same optimized force field also successfully folds both beta-hairpin trpzip2 and mini-protein Trp-Cage, indicating that it is quite robust. Such a balanced, physics-based force field will be highly applicable to a range of biological problems including protein folding and protein structural dynamics.     

An extensible and systematic force field,ESFF, for molecular modeling of organic,inorganic, and organometallic systems

ESFF is a rule-based force field designed for modeling organic, inorganic, and organometallic systems. To cover this broad range of molecular systems, ESFF was developed in an extensible and systematic manner. Several unique features were introduced including pseudoangle and a dot product function representing torsion energy terms. The partial atomic charges that are topology-dependent are determined from ab initio (DFT) calculated electronegativity and hardness for valence orbitals. The van der Waals parameters are charge-dependent, and correlated with the ionization potential for atoms in various valence states. To obtain a set of well-defined and physically meaningful parameters, ESFF employs semiempirical rules to translate atomic-based parameters to parameters typically associated with a covalent valence force field. The atomic parameters depend not only on atom type, but also on internal type, thus resulting in a more accurate force field. This article presents the theory and the method used to develop the force field. The force field has been applied to molecular simulations of a wide variety of systems including nucleic acids, peptides, hydrocarbons, porphyrins, transition metal complexes, zeolites, and organometallic compounds. Agreement with the experimental results indicates that ESFF is a valuable tool in molecular simulations for understanding and predicting both crystal and gas phase molecular structures.     

A theoretical study of the vibrational spectra,geometry and force field of thiourea

A theoretical force field for the molecular vibrations of thiourea has been determined from ab initio calculations at the Hartree-Fock level using the 3-21G* basis set. The reliability of the force field is analyzed by calculating the vibrational frequencies for the deuterated and ¹⁵N isotopomers. Frequencies calculated from the force field are utilized to critically examine the experimental assignments for thiourea and deuterated thiourea. Theoretical geometry, the calculated IR and Raman band intensifies are analyzed.     

An ab initio study of the structure, torsional potential energy function, and electric properties of disilane, ethane, and their deuterated isotopomers

Highly accurate ab initio computations of the molecular structure and properties, torsional potential energy function, and harmonic force field of disilane and ethane have been carried out. Equilibrium parameters as well as vibrational corrections have been evaluated. In addition, for these systems a vibrational averaging procedure has been employed for calculating the dipole moment of molecules which have no permanent dipole moment, i.e., SiH(3)SiD(3) and CH(3)CD(3). The molecular and spectroscopic properties calculated for ethane and its isotopomers provide a calibration against known experimental data, allowing us to estimate the reliability of our computed results for disilane for which there is much less experimental data. The goal of the present study is to predict the molecular parameters, with estimated uncertainties, that determine the microwave spectrum of SiH(3)SiD(3).     

A new force field (ECEPP-05) for peptides, proteins, and organic molecules

Parametrization and testing of a new all-atom force field for organic molecules and peptides with fixed bond lengths and bond angles are described. The van der Waals parameters for both the organic molecules and the peptides were taken from J. Phys. Chem. B 2003, 107, 7143 and J. Phys. Chem. B 2004, 108, 12181. First, the values of the 1-4 nonbonded and electrostatic scale factors appropriate to the new force field were determined by computing the conformational energies of six model molecules, namely, ethanol, ethylamine, propanol, propylamine, 1,2-ethanediol, and 1,3-propanediol with different values of these factors. The partial atomic charges of these molecules were obtained by fitting to the electrostatic potentials calculated with the HF/6-31G quantum-mechanical method. Two different charge models (single- and multiple-conformation-derived) were also considered. We demonstrated that the charge model has a stronger effect on the conformational energies than the 1-4 scaling. The choice of a charge model affected the conformational energies of even the smallest molecules considered, whereas the effect of the 1-4 electrostatic or nonbonded scaling was apparent only for 1,3-propanediol. The best agreement with high-level ab initio data was obtained with the multiple-conformation-derived charges and with no scaling of the 1-4 nonbonded or electrostatic interactions (scale factors of 1.0). Next, the torsional parameters of a large number of neutral and charged organic molecules, assumed to be models of the side chains of the 20 naturally occurring amino acids, were computed by fitting to rotational energy profiles obtained from ab initio MP2/6-31G calculations. The quality of the fits was high with average errors for torsional profiles of less than 0.2 kcal/mol. To derive the torsional parameters for the peptide backbone, the partial atomic charges of the 20 neutral and charged amino acids were obtained by fitting to the electrostatic potentials of terminally blocked amino acids using the HF/6-31G quantum-mechanical method. Then, the phi-psi energy maps of Ac-Ala-NMe and Ac-Gly-NMe were computed using MP2/6-31G//HF/6-31G quantum-mechanical methods. The phi-psi energy map of Ac-Ala-NMe was used for refinement of the nonbonded parameters for the backbone nitrogen and hydrogen bonded to it. Subsequently, the main-chain torsional parameters were obtained by fitting the molecular mechanics energies to the phi-psi energy maps of Ac-Ala-NMe and Ac-Gly-NMe. The transferability of the entire force field was demonstrated by reproducing the main energy minima of terminally blocked Ala3 from the literature. The performance of the force field was also evaluated by simulating crystal structures of small peptides. By comparison of simulated and experimental data, examination of the torsional-angle and atom-positional root-mean-square deviations of the energy-minimized crystal structures from the corresponding X-ray model structures demonstrated high accuracy of the force field.