MM2 Force Field Research Articles

Mechanism of proteolysis in matrix metalloproteinase-2 revealed by QM/MM modeling.

The mechanism of enzymatic peptide hydrolysis in matrix metalloproteinase-2 (MMP-2) was studied at atomic resolution through quantum mechanics/molecular mechanics (QM/MM) simulations. An all-atom three-dimensional molecular model was constructed on the basis of a crystal structure from the Protein Data Bank (ID: 1QIB), and the oligopeptide Ace-Gln-Gly∼Ile-Ala-Gly-Nme was considered as the substrate. Two QM/MM software packages and several computational protocols were employed to calculate QM/MM energy profiles for a four-step mechanism involving an initial nucleophilic attack followed by hydrogen bond rearrangement, proton transfer, and C-N bond cleavage. These QM/MM calculations consistently yield rather low overall barriers for the chemical steps, in the range of 5-10 kcal/mol, for diverse QM treatments (PBE0, B3LYP, and BB1K density functionals as well as local coupled cluster treatments) and two MM force fields (CHARMM and AMBER). It, thus, seems likely that product release is the rate-limiting step in MMP-2 catalysis. This is supported by an exploration of various release channels through QM/MM reaction path calculations and steered molecular dynamics simulations.

Hybrid Quantum Mechanics/Molecular Mechanics/Coarse Grained Modeling: A Triple-Resolution Approach for Biomolecular Systems.

We present a hybrid quantum mechanics/molecular mechanics/coarse-grained (QM/MM/CG) multiresolution approach for solvated biomolecular systems. The chemically important active-site region is treated at the QM level. The biomolecular environment is described by an atomistic MM force field, and the solvent is modeled with the CG Martini force field using standard or polarizable (pol-CG) water. Interactions within the QM, MM, and CG regions, and between the QM and MM regions, are treated in the usual manner, whereas the CG-MM and CG-QM interactions are evaluated using the virtual sites approach. The accuracy and efficiency of our implementation is tested for two enzymes, chorismate mutase (CM) and p-hydroxybenzoate hydroxylase (PHBH). In CM, the QM/MM/CG potential energy scans along the reaction coordinate yield reaction energies that are too large, both for the standard and polarizable Martini CG water models, which can be attributed to adverse effects of using large CG water beads. The inclusion of an atomistic MM water layer (10 Å for uncharged CG water and 5 Å for polarizable CG water) around the QM region improves the energy profiles compared to the reference QM/MM calculations. In analogous QM/MM/CG calculations on PHBH, the use of the pol-CG description for the outer water does not affect the stabilization of the highly charged FADHOOH-pOHB transition state compared to the fully atomistic QM/MM calculations. Detailed performance analysis in a glycine-water model system indicates that computation times for QM energy and gradient evaluations at the density functional level are typically reduced by 40-70% for QM/MM/CG relative to fully atomistic QM/MM calculations.

Application of the molecular mechanics method to simulation of buckling of single-walled carbon nanotubes

The molecular mechanics (MM) method based on the Modified Morse force field is used to determine the critical buckling parameters and post-critical deformation modes of single-walled carbon nanotubes (SWCNTs) twisted at the edges. Computer simulation of the buckling and post-critical deformation of nanotubes is performed using two versions of the MM method: the standard MM method and a mixed method combining the molecular mechanics and molecular structural mechanics approaches (MM/MSM method). Computer simulation shows that the MM/MSM method gives acceptable critical twisting angles, buckling modes, and post-critical deformed configurations of nanotubes, comparable to the critical twisting angles, modes, and configurations obtained using the standard MM method. More precisely, for the (10,10) armchair SCWNT 12.2919nm long, the solutions obtained by both methods are similar up to the beginning of fracture of the SCWNT (at a twisting angle up to 240° for each edge), and for the (10,0) zigzag SCWNT 16.8980nm long, the solutions obtained by both methods are similar up to the final twisting angle of 360° for each edge. The critical linear twisting angle for the (10,10) armchair SWCNT obtained in the present study is compared with the critical linear twisting angles obtained by other authors using the REBO, DREIDING, and MM3 force fields. It is concluded that the critical linear twisting angle obtained in the present study using the Modified Morse force field is unrealistically underestimated compared to the critical linear twisting angles obtained by other authors using the REBO, DREIDING, and MM3 force fields.

Recent advances toward a general purpose linear-scaling quantum force field.

ConspectusThere is need in the molecular simulation community to develop new quantum mechanical (QM) methods that can be routinely applied to the simulation of large molecular systems in complex, heterogeneous condensed phase environments. Although conventional methods, such as the hybrid quantum mechanical/molecular mechanical (QM/MM) method, are adequate for many problems, there remain other applications that demand a fully quantum mechanical approach. QM methods are generally required in applications that involve changes in electronic structure, such as when chemical bond formation or cleavage occurs, when molecules respond to one another through polarization or charge transfer, or when matter interacts with electromagnetic fields. A full QM treatment, rather than QM/MM, is necessary when these features present themselves over a wide spatial range that, in some cases, may span the entire system. Specific examples include the study of catalytic events that involve delocalized changes in chemical bonds, charge transfer, or extensive polarization of the macromolecular environment; drug discovery applications, where the wide range of nonstandard residues and protonation states are challenging to model with purely empirical MM force fields; and the interpretation of spectroscopic observables. Unfortunately, the enormous computational cost of conventional QM methods limit their practical application to small systems. Linear-scaling electronic structure methods (LSQMs) make possible the calculation of large systems but are still too computationally intensive to be applied with the degree of configurational sampling often required to make meaningful comparison with experiment.In this work, we present advances in the development of a quantum mechanical force field (QMFF) suitable for application to biological macromolecules and condensed phase simulations. QMFFs leverage the benefits provided by the LSQM and QM/MM approaches to produce a fully QM method that is able to simultaneously achieve very high accuracy and efficiency. The efficiency of the QMFF is made possible by partitioning the system into fragments and self-consistently solving for the fragment-localized molecular orbitals in the presence of the other fragment’s electron densities. Unlike a LSQM, the QMFF introduces empirical parameters that are tuned to obtain very accurate intermolecular forces. The speed and accuracy of our QMFF is demonstrated through a series of examples ranging from small molecule clusters to condensed phase simulation, and applications to drug docking and protein–protein interactions. In these examples, comparisons are made to conventional molecular mechanical models, semiempirical methods, ab initio Hamiltonians, and a hybrid QM/MM method. The comparisons demonstrate the superior accuracy of our QMFF relative to the other models; nonetheless, we stress that the overarching role of QMFFs is not to supplant these established computational methods for problems where their use is appropriate. The role of QMFFs within the toolbox of multiscale modeling methods is to extend the range of applications to include problems that demand a fully quantum mechanical treatment of a large system with extensive configurational sampling.

Applications of density functional theory to iron-containing molecules of bioinorganic interest.

The past decades have seen an explosive growth in the application of density functional theory (DFT) methods to molecular systems that are of interest in a variety of scientific fields. Owing to its balanced accuracy and efficiency, DFT plays particularly useful roles in the theoretical investigation of large molecules. Even for biological molecules such as proteins, DFT finds application in the form of, e.g., hybrid quantum mechanics and molecular mechanics (QM/MM), in which DFT may be used as a QM method to describe a higher prioritized region in the system, while a MM force field may be used to describe remaining atoms. Iron-containing molecules are particularly important targets of DFT calculations. From the viewpoint of chemistry, this is mainly because iron is abundant on earth, iron plays powerful (and often enigmatic) roles in enzyme catalysis, and iron thus has the great potential for biomimetic catalysis of chemically difficult transformations. In this paper, we present a brief overview of several recent applications of DFT to iron-containing non-heme synthetic complexes, heme-type cytochrome P450 enzymes, and non-heme iron enzymes, all of which are of particular interest in the field of bioinorganic chemistry. Emphasis will be placed on our own work.

The effects of conformation and intermolecular hydrogen bonding on the structural and vibrational spectral data of naproxen molecule

The structural and vibrational properties of naproxen, an inhibitor of cyclooxygenase (COX) enzyme, were investigated by molecular modeling and experimental IR and Raman spectroscopic techniques. Possible conformers of the molecule were searched via a molecular dynamics simulation carried out with MM2 force field. The total energies, equilibrium geometries, force fields, IR and Raman spectral data of the found stable conformers were determined by means of geometry optimization and harmonic frequency calculations carried out using the B3LYP method and Pople-style basis sets of different size. The stability order obtained for the lowest-energy conformers was confirmed by high-accuracy thermochemistry calculations performed with G3MP2B3 composite method. Some electronic structure parameters of naproxen and the anharmonicity characters of its vibrational modes were determined by means of natural population analysis (NPA) and anharmonic frequency calculations at B3LYP/6-31++G(d,p) and B3LYP/6-311++G(d,p) levels of theory. A part of these calculations carried out for free naproxen molecule were repeated also for its energetically most favored dimer forms. Two different scaling procedures ((1) “SQM-FF methodology” and (2) “Dual scale factors”) were independently applied to the obtained harmonic vibrational spectral data to fit them to the corresponding experimental data. In the light of the obtained calculation results, which confirm the remarkable effects of conformation and intermolecular hydrogen bonding on the structural and vibrational spectral data, in particular, on those associated with the functional groups in the propanoic acid chain, a reliable assignment of the fundamental bands observed in the experimental IR and Raman spectra of the molecule was achieved.

QuanPol: A full spectrum and seamless QM/MM program

The quantum chemistry polarizable force field program (QuanPol) is implemented to perform combined quantum mechanical and molecular mechanical (QM/MM) calculations with induced dipole polarizable force fields and induced surface charge continuum solvation models. The QM methods include Hartree-Fock method, density functional theory method (DFT), generalized valence bond theory method, multiconfiguration self-consistent field method, Møller-Plesset perturbation theory method, and time-dependent DFT method. The induced dipoles of the MM atoms and the induced surface charges of the continuum solvation model are self-consistently and variationally determined together with the QM wavefunction. The MM force field methods can be user specified, or a standard force field such as MMFF94, Chemistry at Harvard Molecular Mechanics (CHARMM), Assisted Model Building with Energy Refinement (AMBER), and Optimized Potentials for Liquid Simulations-All Atom (OPLS-AA). Analytic gradients for all of these methods are implemented so geometry optimization and molecular dynamics (MD) simulation can be performed. MD free energy perturbation and umbrella sampling methods are also implemented.

Electron Momentum Spectroscopy of 1-Butene: A Theoretical Analysis Using Molecular Dynamics and Molecular Quantum Similarity

The results of experimental studies of the valence electronic structure of 1-butene by means of electron momentum spectroscopy (EMS) have been reinterpreted on the basis of molecular dynamical simulations in conjunction with the classical MM3 force field. The computed atomic trajectories demonstrate the importance of thermally induced nuclear dynamics in the electronic neutral ground state, in the form of significant deviations from stationary points on the potential energy surface and considerable variations of the C-C-C-C dihedral angle. These motions are found to have a considerable influence on the computed spectral bands and outer-valence electron momentum distributions. Euclidean distances between spherically averaged electron momentum densities confirm that thermally induced nuclear motions need to be fully taken into account for a consistent interpretation of the results of EMS experiments on conformationally flexible molecules.

Nature of Noncovalent Interactions in Catenane Supramolecular Complexes: Calibrating the MM3 Force Field with ab Initio, DFT, and SAPT Methods

The design and assembly of mechanically interlocked molecules, such as catenanes and rotaxanes, are dictated by various types of noncovalent interactions. In particular, [C-H···O] hydrogen-bonding and π-π stacking interactions in these supramolecular complexes have been identified as important noncovalent interactions. With this in mind, we examined the [3]catenane 2·4PF6 using molecular mechanics (MM3), ab initio methods (HF, MP2), several versions of density functional theory (DFT) (B3LYP, M0X), and the dispersion-corrected method DFT-D3. Symmetry adapted perturbation theory (DFT-SAPT) provides the highest level of theory considered, and we use the DFT-SAPT results both to calibrate the other electronic structure methods, and the empirical potential MM3 force field that is often used to describe larger catenane and rotaxane structures where [C-H···O] hydrogen-bonding and π-π stacking interactions play a role. Our results indicate that the MM3 calculated complexation energies agree qualitatively with the energetic ordering from DFT-SAPT calculations with an aug-cc-pVTZ basis, both for structures dominated by [C-H···O] hydrogen-bonding and π-π stacking interactions. When the DFT-SAPT energies are decomposed into components, we find that electrostatic interactions dominate the [C-H···O] hydrogen-bonding interactions, while dispersion makes a significant contribution to π-π stacking. Another important conclusion is that DFT-D3 based on M06 or M06-2X provides interaction energies that are in near-quantitative agreement with DFT-SAPT. DFT results without the D3 correction have important differences compared to DFT-SAPT, while HF and even MP2 results are in poor agreement with DFT-SAPT.

The unique structure of complete lipopolysaccharide isolated from semi-rough Plesiomonas shigelloides O37 (strain CNCTC 39/89) containing (2S)-O-(4-oxopentanoic acid)-α-d-Glcp (α-d-Lenose)

The complete structure of semi-rough lipopolysaccharide (SR-LPS) of Plesiomonas shigelloides CNCTC 39/89 (serotype O37) has been investigated by 1H and 13C NMR spectroscopy, matrix-assisted laser-desorption/ionization time-of-flight MS, and chemical methods. The following structure of the single unit of the O-antigen has been established: [Display omitted] in which α-d-Lenp is (2S)-O-(4-oxopentanoic acid)-α-d-Glcp residue which has not been found in nature. The absolute configuration of oxopentanoic acid moiety in α-d-Lenose residue was determined by NOESY experiment combined with molecular modeling (MM2 force field). The decasaccharide core is substituted at C-4 of the β-d-Glcp residue with a single pentasaccharide unit. Lipid A is built of a β-d-GlcpN4P-(1→6)-α-d-GlcpN1P disaccharide asymmetrically substituted with fatty acids. It was concluded that the core oligosaccharide and the lipid A are identical with those in P. shigelloides CNCTC 113/92 Niedziela et al. (2002)9 and Lukasiewicz et al. (2006).10

Photoelectron and Electron Momentum Spectroscopy of Tetrahydrofuran from a Molecular Dynamical Perspective

The results of experimental studies of the valence electronic structure of tetrahydrofuran employing He I photoelectron spectroscopy as well as Electron Momentum Spectroscopy (EMS) have been reinterpreted on the basis of Molecular Dynamical simulations employing the classical MM3 force field and large-scale quantum mechanical simulations employing Born-Oppenheimer Molecular Dynamics in conjunction with the dispersion corrected ωB97XD exchange-correlation functional. Analysis of the produced atomic trajectories demonstrates the importance of thermal deviations from the lowest energy path for pseudorotation, in the form of considerable variations of the ring-puckering amplitude. These deviations are found to have a significant influence on several outer-valence electron momentum distributions, as well as on the He I photoelectron spectrum.

Modeling Molecular Crystals by QM/MM: Self-Consistent Electrostatic Embedding for Geometry Optimizations and Molecular Property Calculations in the Solid.

We present an approach to model molecular crystals using an adaptive quantum mechanics/molecular mechanics (QM/MM) based protocol. The molecule of interest (or a larger cluster thereof) is described at an appropriate QM level and is embedded in a large array of MM atoms built up from crystal structure information. The nonbonded MM force field consists of atom-centered point charges and Lennard-Jones potentials using van der Waals parameters from the UFF force field. The point charges are initially derived from a single molecule DFT calculation and are then updated self-consistently in the field of point charges. Additional charges are fitted around the MM cluster to correct for missing long-range electrostatic effects. The geometry of the central complex can then be relaxed by quantum chemical calculations in the surrounding MM reaction field, hence capturing solid-state effects on the geometry. We demonstrate the accuracy of this approach for geometry optimization by successful modeling of the huge gas-to-solid bond contraction of HCN-BF3, the ability to reproduce periodic-DFT quality local geometries of solid VOCl3, and the geometry of [Ru(η(5)-Cp*)(η(3)-CH2CHCHC6H5)(NCCH3)2](2+), a difficult ruthenium allyl complex in the solid state. We further show that this protocol is well suited for subsequent molecular property calculations in the solid state (where accurate relaxed geometries are often required) as exemplified by transition metal NMR and EFG calculations of VOCl3 and a vanadium catechol complex in the solid state.

Equilibrium and Growth Morphology of Oligoacenes: Periodic Bond Chains Analysis of Naphthalene, Anthracene, and Pentacene Crystals

The athermal equilibrium and growth crystal shapes of a series of oligoacenes, namely, naphthalene, anthracene, and pentacene, were simulated, in a vacuum, by using three different sets of empirical potentials (UNI, UFF, and MM3 force fields). By applying the Hartman–Perdok method of the periodic bond chains (PBC), the surface profiles were obtained, providing the specific surface and attachment energy values, both for ideal and relaxed surfaces. Very good agreement among the three force fields employed was observed. From calculations, it ensues that surface relaxation only weakly affects surface and, even more, attachment energies of these semiconductor molecular crystals. It is noteworthy to point out that both equilibrium and growth shapes of these crystals are quite similar when concerning phases belonging to the same point group.

Assessment of Weak Intermolecular Interactions Across QM/MM Noncovalent Boundaries

An assessment of a number of quantum mechanical/molecular mechanical (QM/MM) combinations was performed for weak intermolecular interactions across noncovalent QM/MM 'boundaries'. The popular S22 data set, comprising of a number of weak hydrogen-bonded, dispersion-bound and complexes with mixed interactions was used for the assessment. A range of QM methods was combined with a number of popular MM force fields. The single-point interaction energies, at reference geometries, are presented as deviations to accurate CCSD(T)/CBS reference values. This investigation employed both additive and subtractive QM/MM schemes. The density functional has only a negligible effect; the choice of basis set was also negligible in terms of accuracy. The importance of selecting the most appropriate MM force field for accurately describing interactions at the noncovalent 'boundary' region has a dramatic effect on the accuracy.

On the effect of a variation of the force field, spatial boundary condition and size of the QM region in QM/MM MD simulations

During the past years, the use of combined quantum-classical, QM/MM, methods for the study of complex biomolecular processes, such as enzymatic reactions and photocycles, has increased considerably. The quality of the results obtained from QM/MM calculations is largely dependent on five aspects to be considered when setting up a molecular model: the QM Hamiltonian, the MM Hamiltonian or force field, the boundary and coupling between the QM and MM regions, the size of the QM region and the boundary condition for the MM region. In this study, we systematically investigate the influence of a variation of the molecular mechanics force field and the size of the QM region in QM/MM MD simulations on properties of the photoactive part of the blue light photoreceptor protein AppA. For comparison, we additionally performed classical MD simulations and studied the effect of a variation of the type of spatial boundary condition. The classical boundary conditions and the force field used in a QM/MM MD simulation are shown to have non-neglegible effects upon the structural and energetic properties of the protein which makes it advisable to minimize computational artifacts in QM/MM MD simulations by application of periodic boundary conditions and a thermodynamically calibrated force field. A comparison of the structural and energetic properties of MD simulations starting from two alternative, different X-ray structures for the blue light utilizing flavin protein in its dark state indicates a slight preference of the two force fields used for the so-called Anderson structure over the Jung structure.

A conformational analysis and vibrational spectroscopic investigation on 1,2- bis( o- carboxyphenoxy) ethane molecule

The minima on the potential energy surface of 1,2- bis( o- carboxyphenoxy) ethane ( CPE) molecule in its electronic ground state were searched by a molecular dynamics simulation performed with MM2 force field. For each of the found minimum-energy conformers, the corresponding equilibrium geometry, charge distribution, HOMO–LUMO energy gap, force field, vibrational normal modes and associated IR and Raman spectral data were determined by means of the density functional theory (DFT) based electronic structure calculations carried out by using B3LYP method and various Pople-style basis sets. The obtained theoretical data confirmed the significant effects of the intra- and inter-molecular hydrogen bonding interactions on the conformational structure, force field, and group vibrations of the molecule. The same data have also revealed that two of the determined stable conformers, both of which exhibit pseudo-crown structure, are considerably more favorable in energy to the others and accordingly provide the major contribution to the experimental spectra of CPE. In the light of the improved vibrational spectral data obtained within the “SQM FF” methodology and “Dual Scale Factors” approach for the monomer and dimer forms of these two conformers, a reliable assignment of the fundamental bands observed in the experimental room-temperature IR and Raman spectra of the molecule was given, and the sensitivities of its group vibrations to conformation, substitution and dimerization were discussed.

SYNTHESIS OF VARIOUS DITHIOCARBAMATE ESTERS THROUGH REGIOSELECTIVE THIOLATION OF 2-AMINOTHIAZOLE

Dithiocarbamates are of significant importance as pesticides, vulcanizing agents and pharmaceuticals. In search of efficient regioselective strategy, the various dithiocarbamate esters were synthesized by regioselective thiolation at endo-cyclic nitrogen of 2- aminothiazole. The treatment of Boc-protected 2-aminothiazole with carbon disulphide and alkyl halide in the presence of potassium hydroxide in DMF afforded Boc-dithiocarbamate esters (3a-h). The removal of protecting group with TFA provided corresponding dithiocarbamate esters (4a-h). The formation of regioselective product was confirmed by spectral and molecular mechanical MM2 force field and PM3 approximation MOPAC energy minimization data. The strategy provides an important methodology for regioselective thiolation.

Conformational study of the structure of free 12-thiacrown-4 and some of its cation metal complexes

Conformational search of 12-thiacrown-4, 12t4, was performed using the CONFLEX method and the MMFF94S force field whereby 156 conformations were predicted. Optimized geometries of the 156 predicted conformations were calculated at the HF, B3LYP, CAM-B3LYP, M06, M06L, M062x and M06HF levels using the 6-311G** basis set. The correlation energy was recovered at the MP2 level using the same 6-311G** basis set. Optimized geometries at the MP2/6-311G** level and G3MP2 energies were calculated for some of the low energy conformations. The D4 conformation was predicted to be the ground state conformation at all levels of theory considered in this work. Comparison between the dihedral angles of the predicted conformations indicated that for the stability of 12t4, a SCCS dihedral angle of 180° requirement is more important than a gauche CSCC dihedral angle requirement. Conformational search was performed also for the 12t4–Ag+, Bi3+, Cd2+, Cu+ and Sb3+ cation metal complexes using the CONFLEX method and the CAChe-augmented MM3 and MMFF94S force fields. Conformations with relative energies less than 10 kcal/mol at the MP2/6-31+G*//HF/6-31+G* level, with double zeta quality basis set on the metal cations, were considered for computations at the same levels as those used for free 12t4, using also the 6-311G** basis set. The cc-pVTZ-pp basis set was used for the metal cations. The predicted ground state conformations of the 12t4–Ag+, Bi3+, Cd2+, Cu+ and Sb3+ cation metal complexes are the C4, C4, C4, C2v and C4 conformations, respectively. This is in agreement with the experimental X-ray data for the 12t4–Ag+ and Cd2+ cation metal complexes, but experimentally by X-ray, the 12t4–Bi3+ and Cu+ cation metal complexes have Cs and C4 structures, respectively.

Equilibrium and Growth Morphology of Oligoacenes: Periodic Bond Chains (PBC) Analysis of Tetracene Crystal

The athermal equilibrium and growth shapes of tetracene crystal were calculated, in a vacuum, by using computer codes (CSEHP, GULP, and TINKER) along with three different sets of empirical potentials (UNI, UFF, and MM3 force fields). The surface profiles were obtained by applying the Hartman–Perdok method of the periodic bond chains (PBC). The specific surface and attachment energy values were calculated, both for ideal and relaxed surface profiles. From calculations, it follows that the equilibrium shape of tetracene shows a nearly barrel habit: one form, the {001}, rules over the other ones, and several faces roughly belong to the zone axis perpendicular to the dominant form. Concerning the growth morphology, a plate-like habit with the dominant pinacoid is obtained. An excellent agreement among the three force fields employed was observed. A relevant result emerging from this work is the observation that surface relaxation in general weakly affects surface and attachment energies of molecular crystals.

Shapes of Sulfur, Oxygen, and Nitrogen Mustards

Thorough conformational analyses have been performed on representative sulfur, oxygen, and nitrogen mustards. A total of 23, 18, and 38 unique conformers have been located for SM, OM, and NM, respectively, at the MP2/aug-cc-pVDZ level of theory. Despite the fact that these molecules differ only in the identity of the central heteroatom, comparison of their low energy conformations reveals that the shapes they adopt are distinctive to each molecule. Potential energy surfaces for CH(2)-X (X = S, O, and N-CH(3)) and CH(2)-CH(2) bond rotations are presented and, where possible, compared with dihedral angle distributions observed in crystal structure data. These results were used to benchmark and improve the performance of the MM3 and MMFF94 force fields.