Combined quantum-mechanics/molecular-mechanics (QM/MM) approaches have become the method of choice for modeling reactions in biomolecular systems. Quantum-mechanical (QM) methods are required for describing chemical reactions and other electronic processes, such as charge transfer or electronic excitation. However, QM methods are restricted to systems of up to a few hundred atoms. However, the size and conformational complexity of biopolymers calls for methods capable of treating up to several 100,000 atoms and allowing for simulations over time scales of tens of nanoseconds. This is achieved by highly efficient, force-field-based molecular mechanics (MM) methods. Thus to model large biomolecules the logical approach is to combine the two techniques and to use a QM method for the chemically active region (e.g., substrates and co-factors in an enzymatic reaction) and an MM treatment for the surroundings (e.g., protein and solvent). The resulting schemes are commonly referred to as combined or hybrid QM/MM methods. They enable the modeling of reactive biomolecular systems at a reasonable computational effort while providing the necessary accuracy.

QM/MM Methods for Biomolecular Systems

We propose a new long-range correction scheme that combines generalized-gradient-approximation (GGA) exchange functionals in density-functional theory (DFT) with the ab initio Hartree–Fock exchange integral by using the standard error function. To develop this scheme, we suggest a new technique that constructs an approximate first-order density matrix that corresponds to a GGA exchange functional. The calculated results of the long-range correction scheme are found to support a previous argument that the lack of the long-range interactions in conventional exchange functionals may be responsible for the underestimation of 4s−3d interconfigurational energies of the first-row transition metals and for the overestimation of the longitudinal polarizabilities of π-conjugated polyenes in DFT calculations.

A long-range correction scheme for generalized-gradient-approximation exchange functionals

Molcas: a program package for computational chemistry.

A simple modification of second-order Moller–Plesset perturbation theory (MP2) to improve the description of molecular ground state energies is proposed. The total MP2 correlation energy is partitioned into parallel- and antiparallel-spin components which are separately scaled. The two parameters (scaling factors), whose values can be justified by basic theoretical arguments, have been optimized on a benchmark set of 51 reaction energies composed of 74 first-row molecules. It is found, that the new method performs significantly better than standard MP2: the rms [mean absolute error (MAE)] deviation drops from 4.6 (3.3) to 2.3 (1.8) kcal/mol. The maximum error is reduced from 13.3 to 5.1 kcal/mol. Significant improvements are especially observed for cases which are usually known as MP2 pitfalls while cases already described well with MP2 remain almost unchanged. Even for 11 atomization energies not considered in the fit, uniform improvements [MAE: 8.1 kcal/mol (MP2) versus 3.2 kcal/mol (new)] are found. The results are furthermore compared with those from density functional theory (DFT/B3LYP) and quadratic configuration interaction [QCISD/QCISD(T)] calculations. Also for difficult systems including strong (nondynamical) correlation effects, the improved MP2 method clearly outperforms DFT/B3LYP and yields results of QCISD or sometimes QCISD(T) quality. Preliminary calculations of the equilibrium bond lengths and harmonic vibrational frequencies for ten diatomic molecules also show consistent enhancements. The uniformity with which the new method improves upon MP2, thereby rectifying many of its problems, indicates significant robustness and suggests it as a valuable quantum chemical method of general use.

Improved second-order Møller–Plesset perturbation theory by separate scaling of parallel- and antiparallel-spin pair correlation energies

Some of the new unique features of the MOLCAS quantum chemistry package version 7 are presented inthis report. In particular, the Cholesky decomposition method applied to some quantum chemical methods is described. This approach is used both in the context of a straight forward approximation of the two-electron integrals and in the generation of so-called auxiliary basis sets. The article describes how the method is implemented for most known wave functions models: self-consistent field, density functional theory, 2nd order perturbation theory, complete-active space self-consistent field multiconfigurational reference 2nd order perturbation theory, and coupled-cluster methods. The report further elaborates on the implementation of a restricted-active space self-consistent field reference function in conjunction with 2nd order perturbation theory. The average atomic natural orbital basis for relativistic calculations, covering the whole periodic table, are described and associated unique properties are demonstrated. Furthermore, the use of the arbitrary order Douglas-Kroll-Hess transformation for one-component relativistic calculations and its implementation are discussed. This section especially focuses on the implementation of the so-called picture-change-free atomic orbital property integrals. Moreover, the ElectroStatic Potential Fitted scheme, a version of a quantum mechanics/molecular mechanics hybrid method implemented in MOLCAS, is described and discussed. Finally, the report discusses the use of the MOLCAS package for advanced studies of photo chemical phenomena and the usefulness of the algorithms for constrained geometry optimization in MOLCAS in association with such studies. (c) 2009 Wiley Periodicals, Inc. J Corn put Chem 31: 224-247, 2010

Software news and update MOLCAS 7 : The Next Generation

A redundant internal coordinate system for optimizing molecular geometries is constructed from all bonds, all valence angles between bonded atoms, and all dihedral angles between bonded atoms. Redundancies are removed by using the generalized inverse of the G matrix; constraints can be added by using an appropriate projector. For minimizations, redundant internal coordinates provide substantial improvements in optimization efficiency over Cartesian and nonredundant internal coordinates, especially for flexible and polycyclic systems. Transition structure searches are also improved when redundant coordinates are used and when the initial steps are guided by the quadratic synchronous transit approach. © 1996 by John Wiley & Sons, Inc.

Using redundant internal coordinates to optimize equilibrium geometries and transition states

A procedure that automatically identifies internal rotation modes and rotating groups during the normal mode vibrational analysis is outlined, and an improved approximation to the corrections for the thermodynamic functions is proposed. The identification and the characterization of the internal rotation modes require no user intervention and make extensive use of the information imbedded in the redundant internal coordinates. Rigid-rotor internal rotation modes are obtained by fixing stretching, bending, and out-of-plane bending motions and solving the vibrational problem for the constrained system. Normal vibrational modes corresponding to internal rotations are identified by comparing them with the constrained modes. The atomic composition of the rotating groups is determined automatically and the kinetic energy matrix for internal rotation is given by either the constrained Wilson-G matrix or the Kilpatrick and Pitzer protocol. The potential periodicity, the rotating tops’ symmetry numbers, and the well-multiplicity are obtained using simple rules. These parameters can be altered by user input. An improved analytical approximation to the partition function for a one-dimensional hindered internal rotation has been developed that reproduces the accurate values tabulated by Pitzer and Gwinn to ±0.4% with a maximum error of 2.1%. This approximation is shown to behave better than previously available approximations over a wider range of regimes. The one-dimensional rotor treatment is generalized to give useful approximations to the multidimensional rotor thermodynamic functions that can be a good start for more thorough studies.

Identification and treatment of internal rotation in normal mode vibrational analysis

We have used Almlof and Haser’s Laplace transform idea to eliminate the energy denominator in second-order perturbation theory (MP2) and obtain an energy expression in the atomic orbital basis. We show that the asymptotic computational cost of this method scales quadratically with molecular size. We then define atomic orbital domains such that selective pairwise interactions can be neglected using well-defined thresholding criteria based on the power law decay properties of the long-range contributions. For large molecules, our scheme yields linear scaling computational cost as a function of molecular size. The errors can be controlled in a precise manner and our method reproduces canonical MP2 energies. We present benchmark calculations of polyglycine chains and water clusters containing up to 3040 basis functions.

Linear scaling second-order Moller–Plesset theory in the atomic orbital basis for large molecular systems

We present a reformulation of the coupled cluster equations in the atomic orbital (AO) basis that leads to a linear scaling algorithm for large molecules. Neglecting excitation amplitudes in a screening process designed to achieve a target energy accuracy, we obtain an AO coupled cluster method which is competitive in terms of number of amplitudes with the traditional molecular orbital (MO) solution, even for small molecules. For large molecules, the decay properties of integrals and excitation amplitudes becomes evident and our AO method yields a linear scaling algorithm with respect to molecular size. We present benchmark calculations to demonstrate that our AO reformulation of the many-body electron correlation problem defeats the “exponential scaling wall” that has characterized high-level MO quantum chemistry calculations for many years.

Linear scaling coupled cluster and perturbation theories in the atomic orbital basis

Mapping out a reaction mechanism involves optimizing the reactants and products, finding the transition state and following the reaction path connecting them. Transition states can be difficult to locate and reaction paths can be expensive to follow. We describe an efficient algorithm for determining the transition state, minima and reaction path in a single procedure. Starting with an approximate path represented by N points, the path is iteratively relaxed until one of the N points reached the transition state, the end points optimize to minima and the remaining points converged to a second order approximation of the steepest descent path. The method appears to be more reliable than conventional transition state optimization algorithms, and requires only energies and gradients, but not second derivative calculations. The procedure is illustrated by application to a number of model reactions. In most cases, the reaction mechanism can be described well using 5 to 7 points to represent the transition state...

Philippe Y. Ayala

Papers

Using redundant internal coordinates to optimize equilibrium geometries and transition states

Identification and treatment of internal rotation in normal mode vibrational analysis

Linear scaling second-order Moller–Plesset theory in the atomic orbital basis for large molecular systems

Linear scaling coupled cluster and perturbation theories in the atomic orbital basis

A combined method for determining reaction paths, minima, and transition state geometries