Intra‐ and intermolecular interactions in crystals of polar molecules. A study by the mixed quantum mechanical/molecular mechanical SCMP‐NDDO method

Stabilization energies of crystals of polar molecules were calculated with the recently developed NDDO‐SCMP method that determines the wave function of a subunit embedded in the symmetrical environment constituted by the copies of the subunit. The total stabilization energies were decomposed into four components. The deformation energy is the difference between the energy of the molecule in the geometries adopted in the crystal on the one hand, and in vacuo, on the other hand. Further energy components are derived from the molecular geometry found in the crystal phase. The electrostatic component is the interaction energy of the molecule with the crystal field, corresponding to the charge distribution obtained in vacuo. The polarization component is the energy lowering resulted in the self‐consistent optimization of the wave function in the crystal field. The rest of the stabilization energy is attributed to the dispersion–repulsion component, and is calculated from an empirical potential function. The major novelty of this decomposition scheme is the introduction of the deformation energy. It requires the optimization of the structural parameters, including the molecular geometry, the intermolecular coordinates, and the cell parameters of the crystal. The optimization is performed using the recently implemented forces in the SCMP‐NDDO method, and this new feature is discussed in detail. The calculation of the deformation energy is particularly important to obtain stabilization energies for crystals in which the molecular geometry differs considerably from that corresponding to the energy minimum of the isolated molecule. As an example, crystals of diastereoisomeric salts are investigated. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 1679–1690, 2001     

Generalized hybrid orbital for the treatment of boundary atoms in combined quantum mechanical and molecular mechanical calculations using the semiempirical parameterized model 3 method

The application of combined quantum mechanical (QM) and molecular mechanical methods to large molecular systems requires an adequate treatment of the boundary between the two approaches. In this article, we extend the generalized hybrid orbital (GHO) method to the semiempirical parameterized model 3 (PM3) Hamiltonian combined with the CHARMM force field. The GHO method makes use of four hybrid orbitals, one of which is included in the QM region in self-consistent field optimization and three are treated as auxiliary orbitals that do not participate in the QM optimization, but they provide an effective electric field for interactions. An important feature of the GHO method is that the semiempirical parameters for the boundary atom are transferable, and these parameters have been developed for a carbon boundary atom consistent with the PM3 model. The combined GHO-PM3/CHARMM model has been tested on molecular geometry and proton affinity for a series of organic compounds.Acknowledgement We thank the National Institutes of Health for support of this research.Contribution to the Jacopo Tomasi Honorary Issue     

Electrostatic model for infrared intensities in a spectroscopically determined molecular mechanics force field

A new electrostatic model for the calculation of infrared intensities in molecular mechanics and molecular dynamics is presented. The model is based on atomic charges, atomic charge fluxes, and internal coordinate dipoles and their fluxes. The internal coordinate dipoles are used instead of atomic dipoles, thus simplifying the derivation of parameters. The model is designed to reproduce ab initio dipole derivatives, and the parameters can be obtained by (iterative) transformations from these, or by linear least squares fitting to them. A first application to linear alkanes has been made. For these molecules, the intensities can be predicted with an average accuracy of 30–40%. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 754–768, 1998     

Molecular dynamics and free energy perturbation study of hydride-ion transfer step in dihydrofolate reductase using combined quantum and molecular mechanical model

We used molecular dynamics simulation and free energy perturbation (FEP) methods to investigate the hydride-ion transfer step in the mechanism for the nicotinamide adenine dinucleotide phosphate (NADPH)-dependent reduction of a novel substrate by the enzyme dihydrofolate reductase (DHFR). The system is represented by a coupled quantum mechanical and molecular mechanical (QM/MM) model based on the AM1 semiempirical molecular orbital method for the reacting substrate and NADPH cofactor fragments, the AMBER force field for DHFR, and the TIP3P model for solvent water. The FEP calculations were performed for a number of choices for the QM system. The substrate, 8-methylpterin, was treated quantum mechanically in all the calculations, while the larger cofactor molecule was partitioned into various QM and MM regions with the addition of “link” atoms (F, CH₃, and H). Calculations were also carried out with the entire NADPH molecule treated by QM. The free energies of reaction and the net charges on the NADPH fragments were used to determine the most appropriate QM/MM model. The hydride-ion transfer was also carried out over several FEP pathways, and the QM and QM/MM component free energies thus calculated were found to be state functions (i.e., independent of pathway). A ca. 10 kcal/mol increase in free energy for the hydride-ion transfer with an activation barrier of ca. 30 kcal/mol was calculated. The increase in free energy on the hydride-ion transfer arose largely from the QM/MM component. Analysis of the QM/MM energy components suggests that, although a number of charged residues may contribute to the free energy change through long-range electrostatic interactions, the only interaction that can account for the 10 kcal/mol increase in free energy is the hydrogen bond between the carboxylate side chain of Glu30 (avian DHFR) and the activated (protonated) substrate. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 977–988, 1998     

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     

Polarization energies,transport gap and charge transfer states of organic molecular crystals

A self-consistent calculation of electronic polarization in organic molecular crystals and thin films is presented in terms of charge redistribution in nonoverlapping molecules in a lattice. The polarization energies P₊ and P₋ of a molecular cation and anion are found for anthracene and perelynetetracarboxylic dianhydride (PTCDA), together with binding energies of ion pairs and transport gaps of PTCDA films on metallic substrates. The 500 meV variation of P₊+P₋ with film thickness agrees with experiment, as do calculated dielectric tensors. Comparisons are made to submolecular calculations in crystals.     

Control of molecular weight of the polymers produced during the crystalline-state photopolymerization of diethyl cis,cis-muconate as studied by gel permeation chromatography and scanning electron micrography

The molecular weight and its distribution of the resulting polymers were investigated during the crystalline-state photopolymerization of diethyl cis,cis-muconate (EMU). EMU crystals were prepared by several methods, recrystallization, milling, freeze drying, and precipitation, to obtain the crystals with various sizes of 10⁻⁶ to 10⁻² m. After crystalline-state photopolymerization via a crystal-to-crystal process, polymer crystals were isolated and characterized by optical microphotography and scanning electron micrography. Molecular weight and its distribution were determined by gel permeation chromatography with 1,1,1,3,3,3-hexafluoro-2-propanol and by intrinsic viscosity in trifluoroacetic acid. It was revealed that the size of the EMU crystals depended on the method used for the crystal preparation, and that the molecular weight of the polymer decreased as the crystal size became small. © 1998 John Wiley & Sons, Inc. J Polym Sci A: Polym Chem 36: 3147–3155, 1998     

Modelling the active site of glyceraldehyde-3 phosphate dehydrogenase with the LSCF formalism

In the framework of a theoretical approach to the relationship between structure and reactivity of the catalytic centers of enzymes, glyceraldehyde-3 phosphate dehydrogenase (GAPDH) has been chosen as a model enzyme. In GAPDH, the proximity of His₁₇₆ increases the reactivity of Cys₁₄₉ at neutral pH; however, its presence alone is not sufficient to explain the reactivity of the catalytic Cys. In order to determine which other interactions play an important role, a study of the geometric and electronic structure of the catalytic site has been made using a hybrid quantum mechanics/molecular mechanics local self-consistent field method. This allows the computation of the electronic properties of amino acid residues in subsystems influenced by other parts of the macromolecule. The quantum subsystem was centered on the Cys₁₄₉ residue of GAPDH. The structures of GAPDH taken from the crystallographic database did not include hydrogen atoms and these had to be added taking into account the fact that, in the active site, His₁₇₆ has three tautomeric forms: δ-His protonated, ε-His protonated and His⁺. The results presented here suggest that the most stable His…Cys system in GAPDH is a strongly hydrogen-bonded Cys₁₄₉ ⁻/His₁₇₆ ⁺ ion pair. Received: 24 March 1998 / Accepted: 3 September 1998 / Published online: 23 November 1998     

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     

New approach for representation of molecular surface

A new algorithm is proposed for approximation to the molecular surface. It starts with a triangular mesh built on an ellipsoid embracing the whole molecular surface. The triangular mesh is obtained from an icosahedron subdivision sphere with highly uniform vertex distribution, and the embracing surface is deflated stepwise to the best adherence of its triangles onto the surface of the molecule. The deflating direction of each vertex of a triangle is defined by the vector normal at this point to the previous deflated embracing surface. Our results show that the speed of the triangulation embracing ellipsoid method and the quality of the surface obtained by the method are faster and better than the method that starts with a quadrilateral mesh built from meridian and parallel representations on an embracing sphere to get the molecular surface. Furthermore, the surface obtained by the method can be used directly to approximate the molecular surface by spherical harmonic expansions. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1805–1815, 1998     

Quantum mechanical computations on very large molecular systems: The local self-consistent field method

Quantum chemical computations on a subset of a large molecule can be performed, at the neglect of diatomic differential overlap (NDDO) level, without further approximation provided that the atomic orbitals of the frontier atoms are replaced by parametrized orthogonal hybrid orbitals. The electrostatic interaction with the rest of the molecule, treated classically by the usual molecular mechanical approximations, is included into the self-consistent field (SCF) equations. The first and second derivatives of energy are obtained analytically, allowing the search for energy minima and transition states as well as the resolution of Newton equations in molecular dynamics simulations. The local self-consistent field (LSCF) method based on these approximations is tested by studying the intramolecular proton transfer in a Gly-Arg-Glu-Gly model tetrapeptide, which reveals an excellent agreement between a computation performed on the whole molecule and the results obtained by the present method, especially if the quantum subsystem includes the side chains and the peptidic unit in between. The merits of the LSCF method are examplified by a study of proton transfer in the Asp⁶⁹—Arg⁷¹ salt bridge in dihydrofolate reductase. Simulations of large systems, involving local changes of electronic structure, are therefore possible at a good degree of approximation by introducing a quantum chemical part in molecular dynamics studies. This methodology is expected to be very useful for reactivity studies in biomolecules or at the surface of covalent solids. © 1994 by John Wiley & Sons, Inc.     

Influence of environment on proton-transfer mechanisms in model triads from theoretical calculations

We have carried out theoretical calculations to analyze molecular interactions and proton transfer mechanisms in the formate–imidazole–water system, which may be considered the simplest model of catalytic triads in serine proteases. Computations were carried out at the density functional theory level. The effect of a dielectric environment on energy surfaces is considered using a polarizable continuum model and the self-consistent reaction field approach. The role played by inertial and noninertial polarization of this environment is emphasized. Nonequilibrium solvation effects have been estimated. The results show that there are different reaction mechanisms, concerted or stepwise, that may be competitive, depending on the nature of the molecular environment. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1675–1688, 1998     

Modeling molecular weight distribution of comb-branched graft copolymers

Grafting one type polymer onto a different polymer type may yield a comb-branched copolymer. The branching density has a significant effect on its overall molecular weight distribution. A general model is derived to describe the bivariate distribution of molecular weight and branching density for such comb copolymers. The model is applicable for various grafting mechanisms provided the side chains are randomly grafted onto the backbone. The determining parameters are the molecular weight distributions of backbone and side chains, and the branching density. Analytical expressions are obtained for the cases of the side chains having uniform and Schulz–Zimm distributions. © 1998 John Wiley & Sons, Inc. J Polym Sci B: Polym Phys 36: 705–714, 1998     

Stochastic molecular optimization using generalized simulated annealing

We propose a stochastic optimization technique based on a generalized simulated annealing (GSA) method for mapping minima points of molecular conformational energy surfaces. The energy maps are obtained by coupling a classical molecular force field (THOR package) with a GSA procedure. Unlike the usual molecular dynamics (MD) method, the method proposed in this study is force independent; that is, we obtain the optimized conformation without calculating the force, and only potential energy is involved. Therefore, we do not need to know the conformational energy gradient to arrive at equilibrium conformations. Its utility in molecular mechanics is illustrated by applying it to examples of simple molecules (H₂O and H₂O₃) and to polypeptides. The results obtained for H₂O and H₂O₃ using Tsallis thermostatistics suggest that the GSA approach is faster than the other two conventional methods (Boltzmann and Cauchy machines). The results for polypeptides show that pentalanine does not form a stable α-helix structure, probably because the number of hydrogen bonds is insufficient to maintain the helical array. On the contrary, the icoalanine molecule forms an α-helix structure. We obtain this structure simulating all Φ, Ψ pairs using only a few steps, as compared with conventional methods. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 647–657, 1998     

A self-consistent molecular field theory for aggregates of neutral molecules. I

A self-consistent theory is presented for aggregates of neutral molecules. According to the LCAO Hartree-Fock formalism a set of effective Hartree-Fock equations for molecules in the aggregate is derived. The molecular orbitals of each molecule are to be determined from the effective H-F equation for the molecule in which the interactions between the molecule and the surrounding ones are included as an intermolecular interaction field (molecular field). A self-consistent treatment leads to the molecular orbitals which are self-consistent with the molecular field. By this method, then-molecule problem becomesn times of one-molecule problem.     

Fully polarizable QM/MM calculations: An application to the nonbonded iodine–oxygen interaction in dimethyl‐2‐iodobenzoylphosphonate

The compound dimethyl‐2‐iodobenzoylphosphonate is unusual in that it forms well‐ordered crystals that clearly show short iodine‐oxygen interactions in which both the iodine and the oxygen are in their normal oxidation states. These interactions were studied using a new hybrid quantum mechanical–molecular mechanical approach that employs a polarizable molecular mechanics component. The electric field at the molecular mechanics atoms was calculated from a distributed multipole expansion of the wave function; this induced dipoles on the molecular mechanics atoms. The electrostatic potential in a spherical shell around the induced dipoles was reproduced through induced charges on the atomic center and those bonded to it using an analytical (rather than numerical) procedure. The new atomic charges (induced charges plus permanent charges) were then able to interact with the quantum mechanical entity and polarize the wave function. The procedure was iterated to convergence. The calculations show that the iodine atom becomes more positive in the crystal environment (modeled by a chain of three molecules of dimethyl‐2‐iodobenzoylphosphonate). Thus, while the cooperative effects of the crystal environment may not be the only feature stabilizing this unusual interaction, they do play a significant role in reducing the otherwise unfavourable iodine–oxygen monopole–monopole interaction. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 478–482, 2000     

Electrostatic interactions in peptides. Polarisation effects due to an α-helix

A quantum mechanical/molecular mechanical study of a dodecapeptide made of 11 alanine and one asparagine residues in a helical conformation is carried out by means of the local self-consistent field/molecular mechanical and integrated molecular orbital and molecular mechanics computational schemes. The electronic properties of the asparagine side chain are analysed to extract the influence of electrostatic and induction interaction. One finds that induction may play an important role in the energetic and structural features of the systems modelled with mixed methods. The importance of performing quantum computations which explicitly take account of the electrostatic interactions is pointed out.Contribution to the Jacopo Tomasi Honorary Issue     

Improved convergence in dual-topology free energy calculations through use of harmonic restraints

When calculating free energy differences between two molecular systems by means of molecular dynamics simulation, accessory potential functions can help eliminate uninteresting configurational entropy contributions, improve convergence, and facilitate reversibility. In this work, we demonstrate that the use of a harmonic potential function to restrain key portions of a molecular system in a free energy perturbation dual-topology molecular dynamics approach dramatically improves convergence and precision of the calculation. Limitations of this technique are illustrated, and its use in conjunction with a fixed bond-length constraint is developed. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1278–1283, 1998     

Time-reversible always stable predictor-corrector method for molecular dynamics of polarizable molecules

An improved method for classic molecular dynamics of polarizable molecules is proposed. The method uses a predictor, one evaluation of the electrostatic field per integration step, and relaxation (damping). The self-consistent solution is approximated with error of the second order (with respect to the timestep). The time reversibility (long-time energy conservation) error is of the (2n - 1)th order, where n is the predictor length. The method is easy to implement, efficient, accurate, and suitable for any model of polarizability.