Force matching as a stepping stone to QM/MM CB[8] host/guest binding free energies: a SAMPL6 cautionary tale

Use of quantum mechanical/molecular mechanical (QM/MM) methods in binding free energy calculations, particularly in the SAMPL challenge, often fail to achieve improvement over standard additive (MM) force fields. Frequently, the implementation is through use of reference potentials, or the so-called “indirect approach”, and inherently relies on sufficient overlap existing between MM and QM/MM configurational spaces. This overlap is generally poor, particularly for the use of free energy perturbation to perform the MM to QM/MM free energy correction at the end states of interest (e.g., bound and unbound states). However, by utilizing MM parameters that best reproduce forces obtained at the desired QM level of theory, it is possible to lessen the configurational disparity between MM and QM/MM. To this end, we sought to use force matching to generate MM parameters for the SAMPL6 CB[8] host–guest binding challenge, classically compute binding free energies, and apply energetic end state corrections to obtain QM/MM binding free energy differences. For the standard set of 11 molecules and the bonus set (including three additional challenge molecules), error statistics, such as the root mean square deviation (RMSE) were moderately poor (5.5 and 5.4 kcal/mol). Correlation statistics, however, were in the top two for both standard and bonus set submissions (\(R^{2}\) of 0.42 and 0.26, \(\tau\) of 0.64 and 0.47 respectively). High RMSE and moderate correlation strongly indicated the presence of systematic error. Identifiable issues were ameliorated for two of the guest molecules, resulting in a reduction of error and pointing to strong prospects for the future use of this methodology.     

Importance of van der Waals Interactions in QM/MM Simulations

The importance of accurately treating van der Waals interactions between the quantum mechanical (QM) and molecular mechanical (MM) atoms in hybrid QM/MM simulations has been investigated systematically. First, a set of van der Waals (vdW) parameters was optimized for an approximate density functional method, the self-consistent charge-tight binding density functional (SCC-DFTB) approach, based on small hydrogen-bonding clusters. The sensitivity of condensed phase observables to the SCC-DFTB vdW parameters was then quantitatively investigated by SCC-DFTB/MM simulations of several model systems using the optimized set and two sets of extreme vdW parameters selected from the CHARMM22 forcefield. The model systems include a model FAD molecule in solution and a solvated enediolate, and the properties studied include the radial distribution functions of water molecules around the solute (model FAD and enediolate), the reduction potential of the model FAD and the potential of mean force for an intramolecular proton transfer in the enediolate. Although there are noticeable differences between parameter sets for gas-phase clusters and solvent structures around the solute, thermodynamic quantities in the condensed phase (e.g., reduction potential and potential of mean force) were found to be less sensitive to the numerical values of vdW parameters. The differences between SCC-DFTB/MM results with the three vdW parameter sets for SCC-DFTB atoms were explained in terms of the effects of the parameter set on solvation. The current study has made it clear that efforts in improving the reliability of QM/MM methods for energetical properties in the condensed phase should focus on components other than van der Waals interactions between QM and MM atoms.     

Parameter analysis and refinement toolkit system and its application in MM3 parameterization for phosphine and its derivatives

The multiparameter multistep relaxation (MPMSR) method, a routine within a new suite of parameterization programs entitled parameter analysis and refinement toolkit system (PARTS), was developed to assist in the development of molecular mechanics (MM3 and MM2) force field parameters and represents an ongoing effort in our laboratories to generate more accurate force fields in shorter times. In contrast to other computerized parameterization approaches, this method simulates intuition guided trial-and-error and has been used successfully within our laboratories to develop MM2 and MM3 force fields. The primary aim of this approach is to minimize human inspection time and effort, with simultaneous improvement in the efficiency and accuracy of the parameterization process. In an effort to validate the generality of the MPMSR method, a well parameterized data set of phosphine derivatives was reexamined. With the identical set of training molecules used in the original MM3 phosphine parameterization and with minimal human intervention, MPMSR shortened the process from several months to approximately five days. Although the previous phosphine force field is well parameterized, the newly generated MPMSR set of parameters has achieved an overall better fit to the experimentally observed data and ab initio calculations. © 1996 by John Wiley & Sons, Inc.     

Mixed ab initio QM/MM modeling using frozen orbitals and tests with alanine dipeptide and tetrapeptide

Methodology is discussed for mixed ab initio quantum mechanics/molecular mechanics modeling of systems where the quantum mechanics (QM) and molecular mechanics (MM) regions are within the same molecule. The ab initio QM calculations are at the restricted Hartree–Fock level using the pseudospectral method of the Jaguar program while the MM part is treated with the OPLS force fields implemented in the IMPACT program. The interface between the QM and MM regions, in particular, is elaborated upon, as it is dealt with by “breaking” bonds at the boundaries and using Boys-localized orbitals found from model molecules in place of the bonds. These orbitals are kept frozen during QM calculations. Results from tests of the method to find relative conformational energies and geometries of alanine dipeptides and alanine tetrapeptides are presented along with comparisons to pure QM and pure MM calculations. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1468–1494, 1999     

Molecular mechanics (MM4) study of amines

The MM4 force field has been extended to include aliphatic amines. About 20 amines have been examined to obtain a set of useful molecular mechanics parameters for this class. The vibrational spectra of seven amines (172 frequencies) calculated by MM4 have an overall rms error of 27 cm(-1), compared with corresponding MM4 value of 24 cm(-1) for alkanes. The rms and signed average errors of the moments of inertia of nine simple amines compared with the experimental data were 0.18% and -0.004%, respectively. The heats of formation of 30 amines were also studied. The MM4 weighted standard deviation is 0.41 kcal/mol, compared with experiment. Electronegativity effects occur in the hydrocarbon portion of an amine from the nitrogen, and are accounted for by including electronegativity induced changes in bond lengths and angles, and induced dipole-dipole interactions in the molecule. Negative hyperconjugation results from the presence of the lone pair of electrons on nitrogen, and leads to the Bohlmann bands in the infrared, and also to strong and unusual geometric changes in the molecules (Bohlmann effect), all of which are fairly well accounted for. The conformational energies in amines appear to be less straightforward than those for most other classes of molecules, apparently because of the Bohlmann effect, and these are probably not yet completely understood. In general, the agreement between the MM4 calculated results and the available data is reasonably good.     

Alcohols,ethers, carbohydrates,and related compounds. I. The MM4 force field for simple compounds

Simple alcohols and ethers have been studied with the MM4 force field. The structures of 13 molecules have been well fit using the MM4 force field. Moments of inertia have been fit with rms percentage errors as indicated: 18 moments for ethers, 0.28%; 21 moments for alcohols, 0.22%. Rotational barriers and conformational equilibria have also been examined, and the experimental and ab initio results are reproduced substantially better with MM4 than they were with MM3. Much of the improvement comes from the use of additional interaction terms in the force constant matrix, of which the torsion-bend and torsion-torsion are particularly important. Induced dipoles are included in the calculation, and dipole moments are reasonably well fit. It has been possible for the first time to fit conformational energetic data for both open chain and cyclic alcohols (e.g., propanol and cyclohexanol) with the same parameter set. For vibrational spectra, over a total of 82 frequencies, the rms error is 27 cm(-1), as opposed to 38 cm(-1) with MM3. Both the alpha and beta bond shortening resulting from the presence of the electronegative oxygen atom in the molecule are well reproduced. The electronegativity of the oxygen is sufficient that one must also include not only the alpha and beta electronegativity effects on bond lengths, but also on angle distortions, if structures are to be well reproduced. The heats of formation of 32 alcohols and ethers were fit overall to within experimental error (weighted standard deviation error 0.26 kcal/mol).     

Partial hessian fitting for determining force constant parameters in molecular mechanics

We present a new protocol for deriving force constant parameters that are used in molecular mechanics (MM) force fields to describe the bond‐stretching, angle‐bending, and dihedral terms. A 3 × 3 partial matrix is chosen from the MM Hessian matrix in Cartesian coordinates according to a simple rule and made as close as possible to the corresponding partial Hessian matrix computed using quantum mechanics (QM). This partial Hessian fitting (PHF) is done analytically and thus rapidly in a least‐squares sense, yielding force constant parameters as the output. We herein apply this approach to derive force constant parameters for the AMBER‐type energy expression. Test calculations on several different molecules show good performance of the PHF parameter sets in terms of how well they can reproduce QM‐calculated frequencies. When soft bonds are involved in the target molecule as in the case of secondary building units of metal‐organic frameworks, the MM‐optimized geometry sometimes deviates significantly from the QM‐optimized one. We show that this problem is rectified effectively by use of a simple procedure called Katachi that modifies the equilibrium bond distances and angles in bond‐stretching and angle‐bending terms. © 2016 Wiley Periodicals, Inc.     

Efficient treatment of out-of-plane bend and improper torsion interactions in MM2, MM3, and MM4 molecular mechanics calculations

Simple and very efficient formulas are presented for four-body out-of-plane bend (used in MM2 and MM3 force fields) and improper torsion (used in the MM4 force field) internal coordinates and their first and second derivatives. The use of a small set of bend and stretch intermediates allows for order of magnitude decreases in calculation time for potential energies and their first and second derivatives, which are required in molecular mechanics calculations. The formulas are eminently suitable for use in molecular simulations of systems with complicated bond networks. © 1997 John Wiley & Sons, Inc. J Comput Chem 18 : 1804–1811, 1997     

Paramfit: Automated optimization of force field parameters for molecular dynamics simulations

The generation of bond, angle, and torsion parameters for classical molecular dynamics force fields typically requires fitting parameters such that classical properties such as energies and gradients match precalculated quantum data for structures that scan the value of interest. We present a program, Paramfit, distributed as part of the AmberTools software package that automates and extends this fitting process, allowing for simplified parameter generation for applications ranging from single molecules to entire force fields. Paramfit implements a novel combination of a genetic and simplex algorithm to find the optimal set of parameters that replicate either quantum energy or force data. The program allows for the derivation of multiple parameters simultaneously using significantly fewer quantum calculations than previous methods, and can also fit parameters across multiple molecules with applications to force field development. Paramfit has been applied successfully to systems with a sparse number of structures, and has already proven crucial in the development of the Assisted Model Building with Energy Refinement Lipid14 force field. © 2014 Wiley Periodicals, Inc.     

A new strategy for the evaluation of force parameters from quantum mechanical computations

A new strategy for the determination of force parameters is presented. The equilibrium values appearing in the force field equations representing the “stretching” and “bending” of bonds are directly determined from quantum mechanical calculations without geometrical restrictions. The determination of the force parameters is carried out by means of a rigorous fitting between the quantum mechanic and the molecular mechanical energy variations arising from the perturbation of the geometric variables. The strategy presented here has been incorporated into a computer program named PAPQMD, which was developed in order to provide nonquantum mechanical experts with a powerful tool for the determination of approximate force parameters. The program was developed upon the assumption that force parameters are not universal, but they strongly depend on the molecular environment. This implies that the parametrization procedure should be done in a molecular model close to the molecule or molecules to be studied by means of molecular mechanical or dynamic methods, and consequently, it is no longer supposed that the variation of one geometrical parameter does not affect the rest of the molecular geometry. PAPQMD performs the fitting between molecular mechanics and quantum mechanical energies considering all the perturbations that the modification in one geometric variable causes in all the others, enabling the parametrization even of large molecules. The ability of our method to reproduce experimentally derived force parameters is discussed and compared with the widely used Hopfinger's strategy. The study of the behavior of PAPQMD and Hopfinger's strategies for reproducing the force parameters of two complex molecules demonstrates the superiority of the methodology presented here.     

The determination of force constants from isotopic frequencies using parameter representation

The problem of fixing the molecular force field using isotopic frequencies is approached from a new angle with the help of the parameter technique. A parameter space is defined with the parameters for the molecule and its isotopic substituent as the coordinates. Making use of the constraint that the F elements are invariant under isotopic substitution in the ideal case, it is shown that unique points must exist in the parameter space enabling the determination of the force field. The method is applied to harmonic vibrations of molecules with 2 × 2 secular equations. The displays in the parameter space indicate that two such points exist for all molecules. Thus there can be two sets of force fields capable of reproducing the observed frequencies for the molecule and its isotopic substituents. Checks with mean amplitudes of vibration favour only one set. The calculated values of the force constants compare well with the values from the literature.     

Improving the Q2MM method for transition state force field modeling

The Quantum‐to‐molecular mechanics method (Q2MM) for converting quantum mechanical transition states (TSs) to molecular mechanical minima has been modified to allow a fit to the “natural” reaction mode eigenvalue. The resulting force field gives an improved representation of the energy curvature at the TS, but can potentially give false responses to steric interactions. Ways to address this problem while staying close to the “natural” TS force field are discussed. © 2014 Wiley Periodicals, Inc.     

QM/MM vibrational mode tracking

Vibrational spectroscopy is a powerful tool to investigate the structure and dynamics of biomolecules. When small subsystems of large molecules such as active centers of enzymes are studied, quantum chemical calculations based on quantum mechanics/molecular mechanics (QM/MM) coupling schemes are a valuable means to interpret the spectra. The goal of this work is a methodological pilot study on how to selectively and thus efficiently extract certain vibrational information for extended molecular systems described by QM/MM methods. This is achieved by an extension of the mode tracking algorithm and a comparison with the partial Hessian diagonalization approach. After validating the methodology for the CO stretching vibration of 2-butanone and a delocalized CO stretch in acetylacetone, the stretching and bending modes of the CO ligand in CO myoglobin are tracked. Such systems represent an ideal application for mode tracking, because only a few strongly localized vibrations are sought for, while the large remainder of the molecule is of interest only as far as it affects these local vibrations. This influence is treated exactly by mode tracking.     

Chlorosilanes: Development and application of MM2 force field parameters

The geometries, relative conformational energies, and dipole moments of mono and polychlorosilanes have been calculated using ab initio molecular orbital (MO) theory. Calculations at the HF/3–21G(*) level, with the exception of dipole moments, give reasonable agreement with experimental data. A new MM2 force field for chlorosilanes, which includes terms for bond length shortening and bond angle compression due to the attachment of electronegative Cl atoms, has been developed on the basis of experimental and ab initio results. The new force field is generally successful in predicting structural parameters, but is unable to reproduce the dipole moments of several model systems. While dipole moment predictions are not the authors' main interest, this failure defines a shortcoming in the MM2 method. The new parameters have been applied to problems in the prediction of stereochemistries of cyclic systems, and compared with experimental results where data are available.     

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     

Accurate assessment of the strain energy in a protein-bound drug using QM/MM X-ray refinement and converged quantum chemistry

An ongoing question regarding the energetics of protein‐ligand binding has been; what is the strain energy that a ligand pays (if any) when binding to its protein target? The traditional method to estimate strain energy uses force fields to calculate the energy difference between the ligand bound conformation and its nearest local minimum/global minimum on the gas‐phase or aqueous phase potential energy surface. This makes the implicit assumption that the underlying force field as well as the reference crystal structure is accurate. Herein, we use ibuprofen as a test case and compare MMFF and ab initio QM methods to identify the local and global minimum conformations. Nine low energy conformations were identified with HF/6‐31G* geometry optimization in vacuo. We also obtained highly accurate relative energies for ibuprofen's conformational energy surface based on M06/aug‐cc‐pVXZ (X = D and T), MP2/aug‐cc‐pVXZ (X = D and T) and the MP2/CBS method (with and without solvent corrections). Moreover, we curate and re‐refine the ibuprofen‐protein complex (PDB 2BXG) using QM/MM X‐ray refinement approaches (HF/6‐31G* was the QM method and the MM model was the AMBER force field ff99sb), which were compared with the low energy conformers to calculate the strain energy. The result indicates that there was an 88% reduction in ibuprofen conformation strain using the QM/MM refined structure versus the original PDB ibuprofen conformations. Furthermore, our results indicate that, due to its inherent limitations in estimating electrostatic interactions, force fields are not suitable to gauge strain energy for charged drug molecules like ibuprofen. The present work offers a carefully validated conformational potential energy surface for a drug molecule as well as a reliable QM/MM re‐refined X‐ray structure that can be used to test current structure‐based drug design approaches. © 2011 Wiley Periodicals, Inc. J Comput Chem, 2011     

Molecular mechanics calculations (MM2 and MM3) on enamines and aniline derivatives

The MM2 and MM3 force fields have been extended to cover this class of compounds. Structures, vibrational spectra, and other data for 13 compounds were examined and can be reproduced satisfactorily by MM3. Except for the spectra, the other data can be reproduced somewhat less well by MM2. © 1994 by John Wiley & Sons, Inc.     

Coupled semiempirical quantum mechanics and molecular mechanics (QM/MM) calculations on the aqueous solvation free energies of ionized molecules

The aqueous solvation free energies of ionized molecules were computed using a coupled quantum mechanical and molecular mechanical (QM/MM) model based on the AM1, MNDO, and PM3 semiempirical molecular orbital methods for the solute molecule and the TIP3P molecular mechanics model for liquid water. The present work is an extension of our model for neutral solutes where we assumed that the total free energy is the sum of components derived from the electrostatic/polarization terms in the Hamiltonian plus an empirical “nonpolar” term. The electrostatic/polarization contributions to the solvation free energies were computed using molecular dynamics (MD) simulation and thermodynamic integration techniques, while the nonpolar contributions were taken from the literature. The contribution to the electrostatic/polarization component of the free energy due to nonbonded interactions outside the cutoff radii used in the MD simulations was approximated by a Born solvation term. The experimental free energies were reproduced satisfactorily using variational parameters from the vdW terms as in the original model, in addition to a parameter from the one-electron integral terms. The new one-electron parameter was required to account for the short-range effects of overlapping atomic charge densities. The radial distribution functions obtained from the MD simulations showed the expected H-bonded structures between the ionized solute molecule and solvent molecules. We also obtained satisfactory results by neglecting both the empirical nonpolar term and the electronic polarization of the solute, i.e., by implementing a nonpolarization model. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1028–1038, 1999     

Molecular mechanics (MM3) calculations on lithium amide compounds

The MM3 force field has been extended to deal with the lithium amide molecules that are widely used as efficient catalysts for stereoselective asymmetric synthesis. The MM3 force field parameters have been determined on the basis of the ab initio MP2/6-31G* and/or DFT (B3LYP/6-31G*, B3-PW91/6-31G*) geometry optimization calculations. To evaluate the electronic interactions specific to the lithium amides derived from the diamine molecules properly, the Lewis bonding potential term for the interaction between the lithium atom and the nonbonded adjacent electronegative atom such as nitrogen was introduced into the MM3 force field. The bond dipoles were evaluated correctly from the electronic charges on the atoms calculated by fitting to the electrostatic potential at points selected. The MM3 results on the molecular structures, conformational energies, and vibrational spectra show good agreement with those from the quantum mechanical calculations.