MoBioTools: A toolkit to setup quantum mechanics/molecular mechanics calculations

We present a toolkit that allows for the preparation of QM/MM input files from a conformational ensemble of molecular geometries. The package is currently compatible with trajectory and topology files in Amber, CHARMM, GROMACS and NAMD formats, and has the possibility to generate QM/MM input files for Gaussian (09 and 16), Orca (≥4.0), NWChem and (Open)Molcas. The toolkit can be used in command line, so that no programming experience is required, although it presents some features that can also be employed as a python application programming interface. We apply the toolkit in four situations in which different electronic-structure properties of organic molecules in the presence of a solvent or a complex biological environment are computed: the reduction potential of the nucleobases in acetonitrile, an energy decomposition analysis of tyrosine interacting with water, the absorption spectrum of an azobenzene derivative integrated into a voltage-gated ion channel, and the absorption and emission spectra of the luciferine/luciferase complex. These examples show that the toolkit can be employed in a manifold of situations for both the electronic ground state and electronically excited states. It also allows for the automatic correction of the active space in the case of CASSCF calculations on an ensemble of geometries, as it is shown for the azobenzene derivative photoswitch case.     

A list-processing approach to quantum mechanics

On the basis of several earlier papers undiuview current availability of logic programming facilities, the authors propose a list-processing approach to the modelling of algebraic quantum field theory methods in which the noncommutative algebra of quantum -mechanical operators is emulated by lists. The processing produces reordered sequences of elements of a ring with a unit commutator and generates dynamic structures which for some initial arrangements correspond to partially ordered graphs characterized by resources relations and combinatorial identities. The approach is illustrated by reviewing the simple case of a forced harmonic oscillator. The programming aspects are briefly described.     

Structure and electronic properties of hydrated mesityl oxide: a sequential quantum mechanics/molecular mechanics approach

The hydration of mesityl oxide (MOx) was investigated through a sequential quantum mechanics/molecular mechanics approach. Emphasis was placed on the analysis of the role played by water in the MOx syn–anti equilibrium and the electronic absorption spectrum. Results for the structure of the MOx–water solution, free energy of solvation and polarization effects are also reported. Our main conclusion was that in gas-phase and in low-polarity solvents, the MOx exists dominantly in syn-form and in aqueous solution in anti-form. This conclusion was supported by Gibbs free energy calculations in gas phase and in-water by quantum mechanical calculations with polarizable continuum model and thermodynamic perturbation theory in Monte Carlo simulations using a polarized MOx model. The consideration of the in-water polarization of the MOx is very important to correctly describe the solute–solvent electrostatic interaction. Our best estimate for the shift of the π–π^* transition energy of MOx, when it changes from gas-phase to water solvent, shows a red-shift of −2,520 ± 90 cm⁻¹, which is only 110 cm⁻¹ (0.014 eV) below the experimental extrapolation of −2,410 ± 90 cm⁻¹. This red-shift of around −2,500 cm⁻¹ can be divided in two distinct and opposite contributions. One contribution is related to the syn → anti conformational change leading to a blue-shift of ~1,700 cm⁻¹. Other contribution is the solvent effect on the electronic structure of the MOx leading to a red-shift of around −4,200 cm⁻¹. Additionally, this red-shift caused by the solvent effect on the electronic structure can by composed by approximately 60 % due to the electrostatic bulk effect, 10 % due to the explicit inclusion of the hydrogen-bonded water molecules and 30 % due to the explicit inclusion of the nearest water molecules.     

A sequential molecular mechanics/quantum mechanics study of the electronic spectra of amides

We report gas-phase electronic spectra of formamide, N-methyformamide, acetamide, and N-methylacetamide at 300 K calculated using a combination of classical molecular dynamics and time-dependent density functional theory (TDDFT). In comparison to excitation energies computed using the global minima structures, the valence npi* and pi(nb)pi* states show a significant red-shift of 0.1-0.35 eV, while smaller shifts are found for the n3s and pi(nb)3s Rydberg states. In this work, we have identified the physical origin of these shifts arising from variations of the molecular structure. We present simple relationships between key geometrical parameters and spectral shifts. Consequently, electronic spectra can be generated directly from ground-state structures, without additional quantum chemical calculations. The electronic spectrum of formamide in aqueous solution is computed using TDDFT using an explicit solvent model. This provides a quantitative determination of the condensed-phase spectrum. In general, this study shows that temperature effects can change the predicted excitation energies significantly and demonstrates how electronic spectra at elevated temperatures can be computed in a computationally efficient way.     

Electronic properties of liquid ammonia: a sequential molecular dynamics/quantum mechanics approach

The electronic properties of liquid ammonia are investigated by a sequential molecular dynamics/quantum mechanics approach. Quantum mechanics calculations for the liquid phase are based on a reparametrized hybrid exchange-correlation functional that reproduces the electronic properties of ammonia clusters [(NH3)n; n=1-5]. For these small clusters, electron binding energies based on Green's function or electron propagator theory, coupled cluster with single, double, and perturbative triple excitations, and density functional theory (DFT) are compared. Reparametrized DFT results for the dipole moment, electron binding energies, and electronic density of states of liquid ammonia are reported. The calculated average dipole moment of liquid ammonia (2.05+/-0.09 D) corresponds to an increase of 27% compared to the gas phase value and it is 0.23 D above a prediction based on a polarizable model of liquid ammonia [Deng et al., J. Chem. Phys. 100, 7590 (1994)]. Our estimate for the ionization potential of liquid ammonia is 9.74+/-0.73 eV, which is approximately 1.0 eV below the gas phase value for the isolated molecule. The theoretical vertical electron affinity of liquid ammonia is predicted as 0.16+/-0.22 eV, in good agreement with the experimental result for the location of the bottom of the conduction band (-V 0=0.2 eV). Vertical ionization potentials and electron affinities correlate with the total dipole moment of ammonia aggregates.     

A quantum mechanics/molecular mechanics study of the protein-ligand interaction for inhibitors of HIV-1 integrase

Human immunodeficiency virus type-1 integrase (HIV-1 IN) is an essential enzyme for effective viral replication. Diketo acids such as L-731,988 and S-1360 are potent and selective inhibitors of HIV-1 IN. In this study, we used molecular dynamics simulations, within the hybrid quantum mechanics/molecular mechanics (QM/MM) approach, to determine the protein-ligand interaction energy between HIV-1 IN and L-731,988 and 10 of its derivatives and analogues. This hybrid methodology has the advantage that it includes quantum effects such as ligand polarisation upon binding, which can be very important when highly polarisable groups are embedded in anisotropic environments, as for example in metal-containing active sites. Furthermore, an energy decomposition analysis was performed to determine the contributions of individual residues to the enzyme-inhibitor interactions on averaged structures obtained from rather extensive conformational sampling. Analysis of the results reveals first that there is a correlation between protein-ligand interaction energy and experimental strand transfer into human chromosomes and secondly that the Asn-155, Lys-156 and Lys-159 residues and the Mg(2+) ion are crucial to anti-HIV IN activity. These results may explain the available experimental data.     

Protein/ligand binding free energies calculated with quantum mechanics/molecular mechanics

The calculation of binding affinities for flexible ligands has hitherto required the availability of reliable molecular mechanics parameters for the ligands, a restriction that can in principle be lifted by using a mixed quantum mechanics/molecular mechanics (QM/MM) representation in which the ligand is treated quantum mechanically. The feasibility of this approach is evaluated here, combining QM/MM with the Poisson-Boltzmann/surface area model of continuum solvation and testing the method on a set of 47 benzamidine derivatives binding to trypsin. The experimental range of the absolute binding energy (DeltaG = -3.9 to -7.6 kcal/mol) is reproduced well, with a root-mean-square (RMS) error of 1.2 kcal/mol. When QM/MM is applied without reoptimization to the very different ligands of FK506 binding protein the RMS error is only 0.7 kcal/mol. The results show that QM/MM is a promising new avenue for automated docking and scoring of flexible ligands. Suggestions are made for further improvements in accuracy.     

An efficient method for the calculation of quantum mechanics/molecular mechanics free energies

The combination of quantum mechanics (QM) with molecular mechanics (MM) offers a route to improved accuracy in the study of biological systems, and there is now significant research effort being spent to develop QM/MM methods that can be applied to the calculation of relative free energies. Currently, the computational expense of the QM part of the calculation means that there is no single method that achieves both efficiency and rigor; either the QM/MM free energy method is rigorous and computationally expensive, or the method introduces efficiency-led assumptions that can lead to errors in the result, or a lack of generality of application. In this paper we demonstrate a combined approach to form a single, efficient, and, in principle, exact QM/MM free energy method. We demonstrate the application of this method by using it to explore the difference in hydration of water and methane. We demonstrate that it is possible to calculate highly converged QM/MM relative free energies at the MP2/aug-cc-pVDZ/OPLS level within just two days of computation, using commodity processors, and show how the method allows consistent, high-quality sampling of complex solvent configurational change, both when perturbing hydrophilic water into hydrophobic methane, and also when moving from a MM Hamiltonian to a QM/MM Hamiltonian. The results demonstrate the validity and power of this methodology, and raise important questions regarding the compatibility of MM and QM/MM forcefields, and offer a potential route to improved compatibility.     

Integrating quantum and molecular mechanics

A computer algorithm is developed for integrating density functional quantum mechanics into a molecular mechanics program. The computationally infeasible aspects of the standard LCAO-MO approach (the iterative calculation of eigenvectors and the requirement of orthogonal expansions for the orbitals) are replaced with an efficient use of optimization via the trace theorem of linear algebra. The construction of a basis is also described for expanding the electron density that transforms with the molecular geometry. The combination of the trace method and the basis allow the solution for one configuration of atoms and electrons to be tracked over a wide range of internal conformations. The approach is readily adaptable to being used in the context of an imposed classical field that allows it to be used on part of a macromolecular complex. The initial implementation in the program AMMP is described. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1618–1633, 1999     

Structural forms of green fluorescent protein by quantum mechanics/molecular mechanics calculations

The molecular modeling of structural forms of the green fluorescent protein (GFP) with the Ser65Thr single-site mutation was performed by the quantum mechanics/molecular mechanics (QM/MM) method. Two model systems were constructed based on the crystallographic structure from the Protein Data Bank (PDB entry code 1EMA.) The model systems differ in the initial protonation state of the side chain of the amino acid residue Glu222 near the chromophore. The atomic coordinates of the protein macromolecule corresponding to the equilibrium geometric configurations were determined by total energy minimization using the QM/MM method within the density functional theory approximation PBE0/cc-pVDZ for the quantum subsystem that consists of the chromophore, a water molecule, and the side chains of Arg96, Glu222, and Ser205, and with the parameters of the AMBER force field for the molecular mechanics subsystem. In the analysis of the results, particular attention was given to the hydrogen bond redistribution in the chromophore-containing region of the protein caused by a change in the protonation state of the chromophore. The results obtained from the model containing the initially protonated side chain of Glu222 suggest a new interpretation of the photophysical processes in the green fluorescent protein.     

Development of a model for the rational design of molecular imprinted polymer: Computational approach for combined molecular dynamics/quantum mechanics calculations

A new rational approach for the preparation of molecularly imprinted polymer (MIP) based on the combination of molecular dynamics (MD) simulations and quantum mechanics (QM) calculations is described in this work. Before performing molecular modeling, a virtual library of functional monomers was created containing forty frequently used monomers. The MD simulations were first conducted to screen the top three monomers from virtual library in each porogen-acetonitrile, chloroform and carbon tetrachloride. QM simulations were then performed with an aim to select the optimum monomer and progen solvent in which the QM simulations were carried out; the monomers giving the highest binding energies were chosen as the candidate to prepare MIP in its corresponding solvent. The acetochlor, a widely used herbicide, was chosen as the target analyte. According to the theoretical calculation results, the MIP with acetochlor as template was prepared by emulsion polymerization method using N,N-methylene bisacrylamide (MBAAM) as functional monomer and divinylbenzene (DVB) as cross-linker in chloroform. The synthesized MIP was then tested by equilibrium-adsorption method, and the MIP demonstrated high removal efficiency to the acetochlor. Mulliken charge distribution and ¹H NMR spectroscopy of the synthesized MIP provided insight on the nature of recognition during the imprinting process probing the governing interactions for selective binding site formation at a molecular level. We think the computer simulation method first proposed in this paper is a novel and reliable method for the design and synthesis of MIP.     

Quantum mechanics and mixed quantum mechanics/molecular mechanics simulations of model nerve agents with acetylcholinesterase

 The accurate modeling of biological processes presents major computational difficulties owing to the inherent complexity of the macromolecular systems of interest. Simulations of biochemical reactivity tend to require highly computationally intensive quantum mechanical methods, but localized chemical effects tend to depend significantly on properties of the extended biological environment – a regime far more readily examined with lower-level classical empirical models. Mixed quantum/classical techniques are gaining in popularity as a means of bridging these competing requirements. Here we present results comparing two quantum mechanics/molecular mechanics implementations (the SIMOMM technique of Gordon et al. as implemented in GAMESS, and the ONIOM technique of Morokuma et al. found in Gaussian 98) as performed on the enzyme acetylcholinesterase and model nerve agents. This work represents part of the initial phase of a DoD HPCMP Challenge project in which we are attempting to reliably characterize the biochemical processes responsible for nerve agent activity and inhibition, thereby allowing predictions on compounds unrelated to those already studied. Received: 10 October 2001 / Accepted: 13 November 2002 / Published online: 1 April 2003 Contribution to the Proceedings of the Symposium on Combined QM/MM Methods at the 222nd National Meeting of the American Chemical Society, 2001 Correspondence to: M. M. Hurley e-mail: hurley@arl.army.mil     

Pattern recognition in molecular quantum mechanics

It is shown that the classical concept of an open system does not encompass quantal systems but has to be replaced by the non-Boolean notion of an entangled system. Molecular, chemical, or biological phenomena can be considered to be reduced to a fundamental theory like quantum mechanics only if the fundamental and the phenomenological theories are formally and interpretatively connected, and if the classifications used in the empirical sciences are shown to follow from a single set of fundamental dynamical laws. These conditions enforce a non-statistical and ontic interpretation of quantum mechanics, hence a non-Boolean calculus of propositions. In this interpretation the notion of a world state is well-defined, its Schmidt-decomposition defines a background-dependent model state for molecular systems and creates the phenomena we can observe. To any molecular system there is associated in an objective way a nonnegative number which we call the integrity. The integrity measures the inherent fuzziness of the system concept in a holistic theory, and is used to define recognizable molecular patterns.     

Interpretation of nuclear resonant vibrational spectra of rubredoxin using a combined quantum mechanics and molecular mechanics approach

Nuclear resonant vibrational spectra of the reduced and oxidized form of a mutant of rubredoxin from Pyrococcus abyssii were measured and are compared with simulated spectra that were calculated by a combined quantum mechanics (QM) and molecular mechanics (MM) method. Density functional theory was used for the QM level. Calculations were performed for different models of rubredoxin. Realistic spectra were simulated with reduced models that include at least the iron center, the four cysteins coordinating it, and the residues connected to the cysteins together with a QM layer that comprises the first two coordination shells of the iron center. Larger QM layers did not lead to significant changes of the simulated spectra.     

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.     

Cytochrome P450CAM enzymatic catalysis cycle: a quantum mechanics/molecular mechanics study

The catalytic pathway of cytochrome P450cam is studied by means of a hybrid quantum mechanics/molecular mechanics method. Our results reveal an active role of the enzyme in the different catalytic steps. The protein initially controls the energy gap between the high- and low-spin states in the substrate binding process, allowing thermodynamic reduction by putidaredoxin reductase and molecular oxygen addition. A second electron reduction activates the delivery of protons to the active site through a selective interaction of Thr252 and the distal oxygen causing the O--O cleavage. Finally, the protein environment catalyzes the substrate hydrogen atom abstraction step with a remarkably low free energy barrier ( approximately 8 kcal/mol). Our results are consistent with the effect of mutations on the enzymatic efficacy and provide a satisfactory explanation for the experimental failure to trap the proposed catalytically competent species, a ferryl Fe(IV) heme.     

The non-crossing rule in molecular quantum mechanics

Previous mathematical proofs of the non-crossing rule concerning the potential energy curves of a diatomic molecule are examined and criticised. A new treatment is presented which is simple and mathematically rigorous.     

Mixed basis functions in molecular quantum mechanics

A series of calculations have been carried out for the linear system of four equidistant hydrogen atoms for an internuclear separation of 1.7 a.u. The configuration interaction technique was used, the orbital basis consisting of a mixed set of 1s Slater and floating spherical Gaussian functions. The results obtained are encouraging although the effects of unoptimized non-linear parameters are noticeable.

---

Zusammenfassung Es wurde eine Reihe von Rechnungen für ein lineares System von vier äquidistanten Wasserstoffatomen (Kernabstand jeweils 1.7 a.u.) durchgeführt. Dabei benutzt man die Konfigurationswechselwirkungs-Technik, die Orbitalbasis bestand aus einem gemischten Satz von 1s Slater- und floating sphärischen Gauß-Orbitalen. Die Resultate sind ermutigend, obwohl die Effekte der nicht optimierten nichtlinearen Parameter merklich sind.