The dynamo library for molecular simulations using hybrid quantum mechanical and molecular mechanical potentials

The Dynamo module library has been developed for the simulation of molecular systems using hybrid quantum mechanical (QM) and molecular mechanical (MM) potentials. Dynamo is not a program package but is a library of Fortran 90 modules that can be employed by those interested in writing their own programs for performing molecular simulations. The library supports a range of different types of molecular calculation including geometry optimizations, reaction‐path determinations and molecular dynamics and Monte Carlo simulations. This article outlines the general structure and capabilities of the library and describes in detail Dynamo's semiempirical QM/MM hybrid potential. Results are presented to indicate three particular aspects of this implementation—the handling of long‐range nonbonding interactions, the nature of the boundary between the quantum mechanical and molecular mechanical atoms and how to perform path‐integral hybrid‐potential molecular dynamics simulations. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 1088–1100, 2000     

IR spectra of N-methylacetamide in water predicted by combined quantum mechanical/molecular mechanical molecular dynamics simulations

We applied the combined quantum mechanical (QM)/molecular mechanical (MM) molecular dynamics (MD) simulation method in assessing IR spectra of N-methylacetamide and its deuterated form in aqueous solutions. The model peptide is treated at the Austin Model 1 (AM1) level and the induced dipole effects by the solvent are incorporated in fluctuating solute dipole moments, which are calculated using partial charges from Mulliken population analyses without resorting to any available high-level ab initio dipole moment data. Fourier transform of the solute dipole autocorrelation function produces in silico IR spectra, in which the relative peak intensities and bandwidths of major amide bands are quantitatively compatible with experimental results only when both geometric and electronic polarizations of the peptide by the solvent are dealt with at the same quantum-mechanical level. We cast light on the importance of addressing dynamic charge fluctuations of the solute in calculating IR spectra by comparing classical and QM/MM MD simulation results. We propose the adjustable scaling factors for each amide mode to be directly compared with experimental data.     

Simulated scaling method for localized enhanced sampling and simultaneous "alchemical" free energy simulations: a general method for molecular mechanical, quantum mechanical, and quantum mechanical/molecular mechanical simulations

A potential scaling version of simulated tempering is presented to efficiently sample configuration space in a localized region. The present "simulated scaling" method is developed with a Wang-Landau type of updating scheme in order to quickly flatten the distributions in the scaling parameter lambdam space. This proposal is meaningful for a broad range of biophysical problems, in which localized sampling is required. Besides its superior capability and robustness in localized conformational sampling, this simulated scaling method can also naturally lead to efficient "alchemical" free energy predictions when dual-topology alchemical hybrid potential is applied; thereby simultaneously, both of the chemically and conformationally distinct portions of two end point chemical states can be efficiently sampled. As demonstrated in this work, the present method is also feasible for the quantum mechanical and quantum mechanical/molecular mechanical simulations.     

Classical and quantum mechanical/molecular mechanical molecular dynamics simulations of alanine dipeptide in water: comparisons with IR and vibrational circular dichroism spectra

We have implemented the combined quantum mechanical (QM)/molecular mechanical (MM) molecular dynamics (MD) simulations of alanine dipeptide in water along with the polarizable and nonpolarizable classical MD simulations with different models of water. For the QM/MM MD simulation, the alanine dipeptide is treated with the AM1 or PM3 approximations and the fluctuating solute dipole moment is calculated by the Mulliken population analysis. For the classical MD simulations, the solute is treated with the polarizable or nonpolarizable AMBER and polarizable CHARMM force fields and water is treated with the TIP3P, TIP4P, or TIP5P model. It is found that the relative populations of right-handed alpha-helix and extended beta and P(II) conformations in the simulation trajectory strongly depend on the simulation method. For the QM/MM MD simulations, the PM3/MM shows that the P(II) conformation is dominant, whereas the AM1/MM predicts that the dominant conformation is alpha(R). Polarizable CHARMM force field gives almost exclusively P(II) conformation and other force fields predict that both alpha-helical and extended (beta and P(II)) conformations are populated with varying extents. Solvation environment around the dipeptide is investigated by examining the radial distribution functions and numbers and lifetimes of hydrogen bonds. Comparing the simulated IR and vibrational circular dichroism spectra with experimental results, we concluded that the dipeptide adopts the P(II) conformation and PM3/MM, AMBER03 with TIP4P water, and AMBER polarizable force fields are acceptable for structure determination of the dipeptide considered in this paper.     

Infrared spectrum and absorption coefficient of the cofactor flavin in water

Flavins play a key role as redox cofactors of enzymes involved in important metabolic processes. Moreover, they undergo photochemical reactions as chromophores in sensors of blue light or magnetic field in many organisms. The reaction mechanisms of flavoproteins have been investigated by infrared spectroscopy and theoretical studies. However, basic information on flavins in the infrared spectral range has been missing, such as absorption spectra in water and absorption coefficients. Here, the cofactors flavin adenine dinucleotide (FAD) and flavin mononucleotide (FMN) were investigated in aqueous medium by Fourier transform infrared spectroscopy. Transmission and attenuated total reflection (ATR) configuration were employed in direct comparison. Absorption spectra in the range of 920–1800 cm⁻¹ were determined after accurate subtraction of the contributions from the water vibrations. The important carbonyl vibrations were resolved at 1661 and 1712 cm⁻¹. The absorption spectra may serve as a reference for theoretical and experimental studies on the effect of the microenvironment on the flavin cofactor. Furthermore, the molar absorption coefficient of FAD at 1547 cm⁻¹ was determined to 2200 L mol⁻¹ cm⁻¹ with an integral absorption coefficient of ∼50,000 L mol⁻¹ cm⁻². These values are prerequisite for the determination of reaction yields in flavoproteins from reaction-induced difference spectra.     

Semiempirical treatment of electrostatic potentials and partial charges in combined quantum mechanical and molecular mechanical approaches

A semiempirical treatment of electrostatic potentials and partial charges is presented. These are the basic components needed for the evaluation of electrostatic interaction energies in combined quantum mechanical and molecular mechanical approaches. The procedure to compute electrostatic potentials uses AM1 and MNDO wave functions and is based on one previously suggested by Ford and Wang. It retains the NDDO approximation and is thus both easy to implement and computationally efficient. Partial atomic charges are derived from a semiempirical charge equilibration model, which is based on the principle of electronegativity equalization. Large sets of ab initio restricted Hartee-Fock (RHF/6-31G*) reference data have been used to calibrate the semiempirical models. Applying the final parameters (C, H, N, O), the ab initio electrostatic potentials are reproduced with an average accuracy of 20% (AM1) and 25% (MNDO), respectively, and the ab initio potential derived charges normally to within 0.1 e. In most cases our parameterized models are more accurate than the much more expensive quasi ab initio techniques, which employ deorthogonalized semiempirical wave functions and have generally been preferred in previous applications. © 1996 John Wiley & Sons, Inc.     

Femtosecond dynamics of flavin cofactor in DNA photolyase: radical reduction, local solvation, and charge recombination

We report here our femtosecond studies of the photoreduction dynamics of the neutral radical flavin (FADH) cofactor in E. coli photolyase, a process converting the inactive form to the biologically active one, a fully reduced deprotonated flavin FADH(-). The observed temporal absorption evolution revealed two initial electron-transfer reactions, occurring in 11 and 42 ps with the neighboring aromatic residues of W382 and F366, respectively. The new transient absorption, observed at 550 nm previously in photolyase, was found from the excited-state neutral radical and is probably caused by strong interactions with the adenine moiety through the flavin U-shaped configuration and the highly polar/charged surrounding residues. The solvation dynamics from the locally ordered water molecules in the active site was observed to occur in approximately 2 ps. These ultrafast ordered-water motions are critical to stabilizing the photoreduction product FADH(-) instantaneously to prevent fast charge recombination. The back electron-transfer reaction was found to occur in approximately 3 ns. This slow process, consistent with ultrafast stabilization of the catalytic cofactor, favors photoreduction in photolyase.     

Adaptive partitioning in combined quantum mechanical and molecular mechanical calculations of potential energy functions for multiscale simulations

In many applications of multilevel/multiscale methods, an active zone must be modeled by a high-level electronic structure method, while a larger environmental zone can be safely modeled by a lower-level electronic structure method, molecular mechanics, or an analytic potential energy function. In some cases though, the active zone must be redefined as a function of simulation time. Examples include a reactive moiety diffusing through a liquid or solid, a dislocation propagating through a material, or solvent molecules in a second coordination sphere (which is environmental) exchanging with solvent molecules in an active first coordination shell. In this article, we present a procedure for combining the levels smoothly and efficiently in such systems in which atoms or groups of atoms move between high-level and low-level zones. The method dynamically partitions the system into the high-level and low-level zones and, unlike previous algorithms, removes all discontinuities in the potential energy and force whenever atoms or groups of atoms cross boundaries and change zones. The new adaptive partitioning (AP) method is compared to Rode's "hot spot" method and Morokuma's "ONIOM-XS" method that were designed for multilevel molecular dynamics (MD) simulations. MD simulations in the microcanonical ensemble show that the AP method conserves both total energy and momentum, while the ONIOM-XS method fails to conserve total energy and the hot spot method fails to conserve both total energy and momentum. Two versions of the AP method are presented, one scaling as O(2N) and one with linear scaling in N, where N is the number of groups in a buffer zone separating the active high-level zone from the environmental low-level zone. The AP method is also extended to systems with multiple high-level zones to allow, for example, the study of ions and counterions in solution using the multilevel approach.     

Combined molecular mechanical and quantum mechanical potential study of a nucleophilic addition reaction in solution

The procedure of combined semiempirical quantum mechanical (AM1) and molecular mechanical potential⁷ was used to study the nucleophilic addition of hydroxide to formaldehyde in solution. The gas phase AM1 potential surface is approximately 26 kcal/mol more exothermic than the corresponding ab initio 6-31 + G* calculation results. The free energy profile for the reaction in solution was determined by means of molecular dynamic simulations. The resulting free energy of activation is approximately 5 kcal/mol. The difference of the free energy of solvation between the reactant and the product states is about 38 kcal/mol. As the reaction goes on, the number of hydrogen bonds formed by the hydroxide oxygen with the surrounding water molecules decreases, whereas the number of hydrogen bonds formed by the carbonyl oxygen increases. There is no significant change in the total number of hydrogen bonds between the solute and the solvent molecules, and the average number of these hydrogen bonds is between five and six during the entire reaction process. These results are consistent with previous studies using a model based on ad initio 6-31 + G* calculations in the gas phase. The reaction path in solution is different from the gas phase minimum energy reaction path. When the two reactants are at a large distance, the attack route of the hydroxide anion in solution is close to perpendicular to the formaldehyde plane, whereas in the gas phase the route is collinear with the carbonyl group. These results suggests that although AM1 does not yield accurate energies in the gas phase, valuable insights into the solvent effects can be obtained through computer simulations with this combined potential. This combined procedure could be applied to chemical reactions within macromolecules, in which a quantitative estimation of the effects of the environment would not be easily attainable by another technique. © 1994 by John Wiley & Sons, Inc.     

Dynamic structures of phosphodiesterase-5 active site by combined molecular dynamics simulations and hybrid quantum mechanical/molecular mechanical calculations

Various quantum mechanical/molecular mechanical (QM/MM) geometry optimizations starting from an x-ray crystal structure and from the snapshot structures of constrained molecular dynamics (MD) simulations have been performed to characterize two dynamically stable active site structures of phosphodiesterase-5 (PDE5) in solution. The only difference between the two PDE5 structures exists in the catalytic, second bridging ligand (BL2) which is HO- or H2O. It has been shown that, whereas BL2 (i.e. HO-) in the PDE5(BL2 = HO-) structure can really bridge the two positively charged metal ions (Zn2+ and Mg2+), BL2 (i.e. H2O) in the PDE5(BL2 = H2O) structure can only coordinate Mg2+. It has been demonstrated that the results of the QM/MM geometry optimizations are remarkably affected by the solvent water molecules, the dynamics of the protein environment, and the electronic embedding charges of the MM region in the QM part of the QMM/MM calculation. The PDE5(BL2 = H2O) geometries optimized by using the QM/MM method in different ways show strong couplings between these important factors. It is interesting to note that the PDE5(BL2 = HO-) and PDE5(BL2 = H2O) geometries determined by the QM/MM calculations neglecting these three factors are all consistent with the corresponding geometries determined by the QM/MM calculations that account for all of these three factors. These results suggest the overall effects of these three important factors on the optimized geometries can roughly cancel out. However, the QM/MM calculations that only account for some of these factors could lead to considerably different geometries. These results might be useful also in guiding future QM/MM geometry optimizations on other enzymes.     

Quantum mechanical/molecular mechanical simulations of the Tl(III) ion in water

Classical molecular dynamics (MD) and combined quantum mechanical/molecular mechanical (QM/MM) MD simulations have been performed to investigate the structural and dynamical properties of the Tl(III) ion in water. A six-coordinate hydration structure with a maximum probability of the Tl-O distance at 2.21 A was observed, which is in good agreement with X-ray data. The librational and vibrational spectra of water molecules in the first hydration shell are blue-shifted compared with those of pure liquid water, and the Tl-O stretching force constant was evaluated as 148 Nm(-1). Both structural and dynamical properties show a distortion of the first solvation shell structure. The second shell ligands' mean residence time was determined as 12.8 ps. The Tl(III) ion can be classified as "structure forming" ion; the calculated hydration energy of -986 +/- 9 kcal mol agrees well with the experimental value of -986 kcal mol.     

Semiempirical quantum chemical treatment of the standard reduction potentials of quinone and plastoquinone in water

The standard electrode potential for the quinone (Q)-hydroquinone (QH₂) couple in aqueous acidic media has been explicitly calculated. Molecular geometries of Q and QH₂ have been optimized. Protonation of Q, i.e., the formation of QH⁺ and QH, have been considered. Molecular geometries of these species have been thoroughly optimized. The energy of complexation of these molecules with water have been calculated by optimizing the structures of the hydrated complexes Q · 6H₂O, QH₂ · 6H₂O, QH⁺ · 6H₂O. and QH · 6H₂O. The ion–solvent interaction energy of QH⁺ · 6H₂O, in turn, has been calculated by considering the complex QH⁺ · 6H₂O…? 2H₂O, where the two extra water molecules approach the charge center of the complex QH⁺ · 6H₂O vertically from top and bottom of the quinonoid ring. The standard reduction potential calculated by the CNDO method, 0.8548 V, is somewhat larger than the experimental potential, 0.6998 V, at 25°C. But the INDO value, 0.7085 V, is in excellent agreement with the observed potential. The electrode potential for the plastoquinone (PQ)-plastohydroquinone (PQH₂) couple present in the aqueous pool in chloroplast has been calculated by the INDO method. The basic geometries of PQ, PQH⁺, and PQH_2/sb have been synthesized by adopting the optimized geometries of Q, QH⁺, and QH₂ and considering methyl substituents as well as an isoprenoid side chain containing up to 3 isoprene units with possible geometrical isomerism. The hydrated species PQ · 6H₂O, PQH⁺ · 6H₂O, and PQH₂ · 6H₂O are unstable compared to the isolated species PQ, PQH⁺, and PQH₂, respectively. In fact, we have found that the hydration of PQH⁺ and PQH₂ is much less extensive, and stability arises only when the hydroxyl groups in these two molecules are hydrogen-bonded to water molecules. But PQH⁺ is also stabilized through the association of two more water molecules in the vertical direction. For this reason, we have calculated the reduction potential of the PQ/PQH₂ system from the energies of the isolated molecules PQH₂ and the hydrated species PQH⁺ · 2H₂O. The computed standard reduction potential is 0.2785 V and it yields a potential of 0.07V at pH 7 at 25°C, which is in good agreement with the reduction potential 0.11 V observed for plastoquinone in the aqueous pool in chloroplast. © 1994 John Wiley & Sons, Inc.     

Structure and ultrafast dynamics of liquid water: a quantum mechanics/molecular mechanics molecular dynamics simulations study

A quantum mechanics/molecular mechanics molecular dynamics simulation was performed for liquid water to investigate structural and dynamical properties of this peculiar liquid. The most important region containing a central reference molecule and all nearest surrounding molecules (first coordination shell) was treated by Hartree-Fock (HF), post-Hartree-Fock [second-order Moller-Plesset perturbation theory (MP2)], and hybrid density functional B3LYP [Becke's three parameter functional (B3) with the correlation functional of Lee, Yang, and Parr (LYP)] methods. In addition, another HF-level simulation (2HF) included the full second coordination shell. Site to site interactions between oxygen-oxygen, oxygen-hydrogen, and hydrogen-hydrogen atoms of all ab initio methods were compared to experimental data. The absence of a second peak and the appearance of a shoulder instead in the gO-O graph obtained from the 2HF simulation is notable, as this feature has been observed so far only for pressurized or heated water. Dynamical data show that the 2HF procedure compensates some of the deficiency of the HF one-shell simulation, reducing the difference between correlated (MP2) and HF results. B3LYP apparently leads to too rigid structures and thus to an artificial slow down of the dynamics.     

Conserving energy during molecular dynamics simulations of water,proteins, and proteins in water

Molecular dynamics simulations have been carried out for a series of systems of increasing complexity including: pure water, a model polypeptide (α-helical decaglycine) in vacuo, a protein (Pancreatic Trypsin Inhibitor, PTI) in vacuo, and a fully solvated protein (PTI in water). The equations of motion were integrated using Andersen's velocity version of the Verlet algorithm with internal contraints (the RATTLE algorithm). The accuracy with which the equations of motion are integrated has been analyzed for several different simulation conditions. The effects of various nonbonded interaction truncation schemes on the conservation of energy have been examined, including the use of atomic cutoffs, and (neutral group) residue cutoffs. The use of a smoothing function to eliminate the discontinuities in the potential at the cutoff leads to a significant improvement in the accuracy of the integration for each of the systems studied. The accuracy with which the equations of motion are integrated using the RATTLE algorithm for pure water and for the solvated protein are found to be comparable when the nonbonded interactions are tapered with a smoothing function at the cutoff.     

Structure and dynamics of methanol in water: a quantum mechanical charge field molecular dynamics study

An ab initio quantum mechanical charge field molecular dynamics simulation was carried out for one methanol molecule in water to analyze the structure and dynamics of hydrophobic and hydrophilic groups. It is found that water molecules around the methyl group form a cage-like structure whereas the hydroxyl group acts as both hydrogen bond donor and acceptor, thus forming several hydrogen bonds with water molecules. The dynamic analyses correlate well with the structural data, evaluated by means of radial distribution functions, angular distribution functions, and coordination number distributions. The overall ligand mean residence time, τ identifies the methanol molecule as structure maker. The relative dynamics data of hydrogen bonds between hydroxyl of methanol and water molecules prove the existence of both strong and weak hydrogen bonds. The results obtained from the simulation are in excellent agreement with the experimental results for dilute solution of CH(3)OH in water. The overall hydration shell of methanol consists in average of 18 water molecules out of which three are hydrogen bonded.     

Force-field modeling through quantum mechanical calculations: molecular dynamics simulations of a nematogenic molecule in its condensed phases

Interaction energy of the 4-n-pentyloxy-4'-cyanobiphenyl (5OCB) dimer is computed at MP2 level, for many geometrical arrangements using the Fragmentation Reconstruction Method (FRM). DFT calculations are performed for a number of geometries of the monomer. The resulting database is used to parameterize an atomistic intra- and inter-molecular force-field suitable for classical bulk simulations. Several structural and dynamical properties in 5OCB isotropic and liquid crystalline phases are computed from molecular dynamics simulation mainly in the NPT ensemble. Lengthy runs (more than 70 ns) and large sample sizes (up to 806 molecules) were used to determine the nematic to isotropic transition temperature up to a precision of few K. Good agreement was found in most of the investigated properties, thus validating the accuracy of the proposed model potential, only derived by quantum mechanical calculations.     

Toward a new approach for determination of solute's charge distribution to analyze interatomic electrostatic interactions in quantum mechanical/molecular mechanical simulations

We present an alternative approach to determine "density-dependent property"-derived charges for molecules in the condensed phase. In the case of a solution, it is essential to take into consideration the electron polarization of molecules in the active site of this system. The solute and solvent molecules in this site have to be described by a quantum mechanical technique and the others are allowed to be treated by a molecular mechanical method (QM/MM scheme). For calculations based on this scheme, using the forces and interaction energy as density-dependent property our charges from interaction energy and forces (CHIEF) approach can provide the atom-centered charges on the solute atoms. These charges reproduce well the electrostatic potentials around the solvent molecules and present properly the picture of the electron density of the QM subsystem in the solution system. Thus, the CHIEF charges can be considered as the atomic charges under the conditions of the QM/MM simulation, and then enable one to analyze electrostatic interactions between atoms in the QM and MM regions. This approach would give a view of the QM nuclei and electrons different from the conventional methods.