Linear-scaling quantum mechanical methods for excited states

The poor scaling of many existing quantum mechanical methods with respect to the system size hinders their applications to large systems. In this tutorial review, we focus on latest research on linear-scaling or O(N) quantum mechanical methods for excited states. Based on the locality of quantum mechanical systems, O(N) quantum mechanical methods for excited states are comprised of two categories, the time-domain and frequency-domain methods. The former solves the dynamics of the electronic systems in real time while the latter involves direct evaluation of electronic response in the frequency-domain. The localized density matrix (LDM) method is the first and most mature linear-scaling quantum mechanical method for excited states. It has been implemented in time- and frequency-domains. The O(N) time-domain methods also include the approach that solves the time-dependent Kohn-Sham (TDKS) equation using the non-orthogonal localized molecular orbitals (NOLMOs). Besides the frequency-domain LDM method, other O(N) frequency-domain methods have been proposed and implemented at the first-principles level. Except one-dimensional or quasi-one-dimensional systems, the O(N) frequency-domain methods are often not applicable to resonant responses because of the convergence problem. For linear response, the most efficient O(N) first-principles method is found to be the LDM method with Chebyshev expansion for time integration. For off-resonant response (including nonlinear properties) at a specific frequency, the frequency-domain methods with iterative solvers are quite efficient and thus practical. For nonlinear response, both on-resonance and off-resonance, the time-domain methods can be used, however, as the time-domain first-principles methods are quite expensive, time-domain O(N) semi-empirical methods are often the practical choice. Compared to the O(N) frequency-domain methods, the O(N) time-domain methods for excited states are much more mature and numerically stable, and have been applied widely to investigate the dynamics of complex molecular systems.     

Improved pseudobonds for combined ab initio quantum mechanical/molecular mechanical methods

The pseudobond approach offers a smooth connection at the quantum mechanical/molecular mechanical interface which passes through covalent bonds. It replaces the boundary atom of the environment part with a seven-valence-electron atom to form a pseudobond with the boundary atom of the active part [Y. Zhang, T. S. Lee, and W. Yang, J. Chem. Phys. 110, 46 (1999)]. In its original formulation, the seven-valence-electron boundary atom has the basis set of fluorine and a parametrized effective core potential. Up to now, only the Cps(sp3)-C(sp3) pseudobond has been successfully developed; thus in the case of proteins, it can only be used to cut the protein side chains. Here we employ a different formulation to construct this seven-valence-electron boundary atom, which has its own basis set as well as the effective core potential. We have not only further improved Cps(sp3)-C(sp3) pseudobond, but also developed Cps(sp3)-C(sp2,carbonyl) and Cps(sp3)-N(sp3) pseudobonds for the cutting of protein backbones and nucleic acid bases. The basis set and effective core potential for the seven-valence-electron boundary atom are independent of the molecular mechanical force field. Although the parametrization is performed with density functional calculations using hybrid B3LYP exchange-correlation functional, it is found that the same set of parameters is also applicable to Hartree-Fock and MP2 methods, as well as DFT calculations with other exchange-correlation functionals. Tests on a series of molecules yield very good structural, electronic, and energetic results in comparison with the corresponding full ab initio quantum mechanical calculations.     

Prediction of the mechanisms of enzyme-catalysed reactions using hybrid quantum mechanical/molecular mechanical methods

The use of hybrid methods, involving both quantum mechanics and molecular mechanics, to model the mechanism of enzyme-catalysed reactions, is discussed. Two alternative approaches to treating the electrostatic interactions between the quantum mechanical and molecular mechanical regions are studied, involving either the inclusion of this term in the electronic Hamiltonian (QM/MM), or evaluating it purely classically (MO + MM). In the latter scheme, possible problems of using force fields that are standard for macromolecular modelling are identified. The use of QM/MM schemes to investigate the mechanism of the enzymes thymidine phosphorylase (ThdPase) and protein tyrosine phosphatase (PTP) is described. For both systems, transition states have been identified using a PM3 Hamiltonian. For ThdPase, concerted motion of the enzyme during the course of the reaction is suggested and, for PTP, a two-step dephosphorylation reaction is indicated, both with quite low barriers.     

Comparison of static first hyperpolarizabilities calculated with various quantum mechanical methods

The prediction of nonlinear electro-optic (EO) behavior of molecules with quantum methods is the first step in the development of organic-based electro-optic devices. Typical EO molecules may require calculations with several hundred electrons, which prevents all but the fastest methods (semiempirical and density functional theory (DFT)) from being used for EO estimation. To test the reliability of these methods, we compare dipole moments, polarizabilities, and first-order hyperpolarizabilities for a wide range of structures of experimental interest with Hartree-Fock (HF), intermediate neglect of differential overlap (INDO), and DFT methods. The relative merits of molecules are consistently predictable with every method.     

Performance assessment of semiempirical molecular orbital methods in describing halogen bonding: quantum mechanical and quantum mechanical/molecular mechanical-molecular dynamics study

The performance of semiempirical molecular-orbital methods--MNDO, MNDO-d, AM1, RM1, PM3 and PM6--in describing halogen bonding was evaluated, and the results were compared with molecular mechanical (MM) and quantum mechanical (QM) data. Three types of performance were assessed: (1) geometrical optimizations and binding energy calculations for 27 halogen-containing molecules complexed with various Lewis bases (Two of the tested methods, AM1 and RM1, gave results that agree with the QM data.); (2) charge distribution calculations for halobenzene molecules, determined by calculating the solvation free energies of the molecules relative to benzene in explicit and implicit generalized Born (GB) solvents (None of the methods gave results that agree with the experimental data.); and (3) appropriateness of the semiempirical methods in the hybrid quantum-mechanical/molecular-mechanical (QM/MM) scheme, investigated by studying the molecular inhibition of CK2 protein by eight halobenzimidazole and -benzotriazole derivatives using hybrid QM/MM molecular-dynamics (MD) simulations with the inhibitor described at the QM level by the AM1 method and the rest of the system described at the MM level. The pure MM approach with inclusion of an extra point of positive charge on the halogen atom approach gave better results than the hybrid QM/MM approach involving the AM1 method. Also, in comparison with the pure MM-GBSA (generalized Born surface area) binding energies and experimental data, the calculated QM/MM-GBSA binding energies of the inhibitors were improved by replacing the G(GB,QM/MM) solvation term with the corresponding G(GB,MM) term.     

A complete quantum mechanical study of chlorine photodissociation

A fully quantum mechanical dynamical calculation on the photodissociation of molecular chlorine is presented. The magnitudes and phases of all the relevant photofragment T-matrices have been calculated, making this study the computational equivalent of a "complete experiment," where all the possible parameters defining an experiment have been determined. The results are used to simulate cross-sections and angular momentum polarization information which may be compared with experimental data. The calculations rigorously confirm the currently accepted mechanism for the UV photodissociation of Cl(2), in which the majority of the products exit on the C(1)Π(1u) state, with non-adiabatic couplings to the A(3)Π(1u) and several other Ω = 1 states, and a small contribution from the B(3)Π state present at longer wavelengths.     

Comparison of linear-scaling semiempirical methods and combined quantum mechanical/molecular mechanical methods for enzymic reactions. II. An energy decomposition analysis

QM/MM methods have been developed as a computationally feasible solution to QM simulation of chemical processes, such as enzyme-catalyzed reactions, within a more approximate MM representation of the condensed-phase environment. However, there has been no independent method for checking the quality of this representation, especially for highly nonisotropic protein environments such as those surrounding enzyme active sites. Hence, the validity of QM/MM methods is largely untested. Here we use the possibility of performing all-QM calculations at the semiempirical PM3 level with a linear-scaling method (MOZYME) to assess the performance of a QM/MM method (PM3/AMBER94 force field). Using two model pathways for the hydride-ion transfer reaction of the enzyme dihydrofolate reductase studied previously (Titmuss et al., Chem Phys Lett 2000, 320, 169-176), we have analyzed the reaction energy contributions (QM, QM/MM, and MM) from the QM/MM results and compared them with analogous-region components calculated via an energy partitioning scheme implemented into MOZYME. This analysis further divided the MOZYME components into Coulomb, resonance and exchange energy terms. For the model in which the MM coordinates are kept fixed during the reaction, we find that the MOZYME and QM/MM total energy profiles agree very well, but that there are significant differences in the energy components. Most significantly there is a large change (approximately 16 kcal/mol) in the MOZYME MM component due to polarization of the MM region surrounding the active site, and which arises mostly from MM atoms close to (<10 A) the active-site QM region, which is not modelled explicitly by our QM/MM method. However, for the model where the MM coordinates are allowed to vary during the reaction, we find large differences in the MOZYME and QM/MM total energy profiles, with a discrepancy of 52 kcal/mol between the relative reaction (product-reactant) energies. This is largely due to a difference in the MM energies of 58 kcal/mol, of which we can attribute approximately 40 kcal/mol to geometry effects in the MM region and the remainder, as before, to MM region polarization. Contrary to the fixed-geometry model, there is no correlation of the MM energy changes with distance from the QM region, nor are they contributed by only a few residues. Overall, the results suggest that merely extending the size of the QM region in the QM/MM calculation is not a universal solution to the MOZYME- and QM/MM-method differences. They also suggest that attaching physical significance to MOZYME Coulomb, resonance and exchange components is problematic. Although we conclude that it would be possible to reparameterize the QM/MM force field to reproduce MOZYME energies, a better way to account for both the effects of the protein environment and known deficiencies in semiempirical methods would be to parameterize the force field based on data from DFT or ab initio QM linear-scaling calculations. Such a force field could be used efficiently in MD simulations to calculate free energies.     

A quantum mechanical approach to electrochemical behavior of spirochromics

Fully optimized semiempirical quantum‐chemical calculations of photochromic spiropyrans are presented. The vertical ionization potentials are calculated and their variation with substitutions are correlated to experimental oxidation potentials. The effects of the substitutions are studied and the partial charges on indoline and pyran components generated by HOMO are found to be responsible for the variations. The deactivating groups on the indoline ring system and deactivating groups on the pyran system increase the ionization potential and, consecutively, the oxidation potential. ©1999 John Wiley & Sons, Inc. Int J Quant Chem 75: 111–117, 1999     

A systematic treatment of three-dimensional quantum mechanical reaction coordinates

Summary The rigorous, collinear, canonical point transformation method with hyper-hyperbolic coordinates is extended to the infinite central mass problem in three dimensions. The initial transformation performed is (x_A, y _A, z _A, x _C, y _C, z _C) (, , , r, R, ), where (, , ) are the Euler angles; r and R are the AB and BC interatomic distances, respectively, and is the angle between r and R. A second transformation is then performed to (, , , , , ), where is the reaction coordinate mimicking the reaction path, and is the vibrational coordinate of the diatom. The transformed spaces are all one-to-one mappings from the original spaces, and thus do not have any three-to-one regions. The transformed momenta and Hamiltonians are derived, and are Hermitian in their respective transformed spaces.     

Comparison of quantum mechanical methods to compute the biologically relevant reactivities of cyclopenta-polycyclic aromatic hydrocarbons

In computational studies to understand the interaction of polycyclic aromatic hydrocarbons (PAHs) with biomolecular systems, the semiempirical method AM 1 has been used previously to determine the geometry of the PAH and its metabolites and relevant intermediates. A number of studies have shown that AM 1 provides geometries for parent PAHs that are acceptably close to experimentally determined structures. However, many of the properties that determine the manner by which PAHs interact with biological nucleophiles depend on the structure of metabolites and reactive intermediates where less experimental information is available. In a previous study, we used AM 1 to obtain the molecular geometries of reactive intermediates of cyclopenta-PAHs (cPAHs) and then used single-point Hartree-Fock calculations, with the gaussian 3-21g basis set, to obtain molecular energies and charge distributions, in order to predict the direction of epoxide ring opening. Recent advances in the availability of computational hardware and software have provided other, more rigorous, methods for approaching this problem. In this study, we used hartree-fock methods in the gaussian series of programs employing the 3-21_g and 6-31_g basis sets and the local density functional method Dmol to obtain molecular geometries, energies, and charge distributions of the epoxides and the two potential hydroxycarbocations that could result from protonated ring opening, for a series of cPAHs. We have also performed the same calculations with AMSOL /SM 2, a semiempirical method that adds the effect of the aqueous environment to the AM 1 Hamiltonian. The division of the cPAHs into classes is not altered by these more rigorous calculations. The inclusion of water in the Hamiltonian has a greater effect on the results than using the ab initio methods to obtain the structure. © 1994 John Wiley & Sons, Inc.     

Pairwise decomposition of residue interaction energies using semiempirical quantum mechanical methods in studies of protein-ligand interaction

Pairwise decomposition of the interaction energy between molecules is shown to be a powerful tool that can increase our understanding of macromolecular recognition processes. Herein we calculate the pairwise decomposition of the interaction energy between the protein human carbonic anhydrase II (HCAII) and the fluorine-substituted ligand N-(4-sulfamylbenzoyl)benzylamine (SBB) using semiempirical quantum mechanics based methods. We dissect the interaction between the ligand and the protein by dividing the ligand and the protein into subsystems to understand the structure-activity relationships as a result of fluorine substitution. In particular, the off-diagonal elements of the Fock matrix that is composed of the interaction between the ionic core and the valence electrons and the exchange energy between the subsystems or atoms of interest is examined in detail. Our analysis reveals that the fluorine-substituted benzylamine group of SBB does not directly affect the binding energy. Rather, we find that the strength of the interaction between Thr199 of HCAII and the sulfamylbenzoyl group of SBB affects the binding affinity between the protein and the ligand. These observations underline the importance of the sulfonamide group in binding affinity as shown by previous experiments (Maren, T. H.; Wiley: C. E. J. Med. Chem. 1968, 11, 228-232). Moreover, our calculations qualitatively agree with the structural aspects of these protein-ligand complexes as determined by X-ray crystallography.     

Generalized quantum mechanical two-centre problems

If χ _i (χ _k ) is an exact generalized diatomic orbital (solution of Eq. (1) of text), a sequence of functions χ _i ^(N) converging to χ _i may be constructed so that matrix elements of frequently occurring operators between χ _i ^(N) and χ _k ^(N) may be computed without any numerical integration. Exact expectation values are given for kinetic and potential energy, dipole moment, θ ²=x ²+y ², and quadrupole moment 3z ²?r ², for various ratios of nuclear charges Z ₁,Z ₂ and for several distances R. Special subjects discussed in terms of computed expectation values are:

1. R-dependence of the contributions to total energy of HeH²⁺ in state 2pσ and of LiH³⁺ in state 3dσ
2. RZ-and λ-dependence of dipole and quadrupole moment functions in state 1sσ
3. Some properties of those generalized diatomic orbitals which approach, for R going to 0, Slater-type atomic functions.

     

Generalized quantum mechanical two-centre problems

Correlation diagrams for the lowest electronic states of the systems (LiH)³⁺, (HeH)⁺⁺, (LiHe)⁴⁺, H ₂ ⁺ , (He, –1), and (H, –1)^– (finite dipole with one electron) have been computed exactly and are discussed with special regard to the non-crossing rule and to the asymptotic behaviour of generalized diatomic orbitals.

---

Zusammenfassung Korrelationsdiagramme für die tiefsten elektronischen Zustände der Systeme (LiH)³⁺, (HeH)⁺⁺, (LiHe)⁴⁺, H ₂ ⁺ , (He, –1) und (H, –1)^– (endlicher Dipol mit einem Elektron) wurden exakt berechnet und werden im Hinblick auf die Nichtüberschneidungsregel und auf das asymptotische Verhalten verallgemeinerter zweiatomiger Bahnfunktionen diskutiert.

---

Résumé On a calculé exactement les diagrammes de corrélation pour les plus bas états électroniques des systèmes (LiH)³⁺, (HeH)⁺⁺, (LiHe)⁴⁺, H ₂ ⁺ , (He, –1) et (H, –1) (dipole fini à 1 électron). La discussion porte plus spécialement sur la loi de non-intersection et sur le comportement asymptotique des orbitales diatomiques généralisées.

---

---

Dedicated to Prof. Dr. Hermann Dänzer on occasion of his 65 th birthday on October 21 st.

---

The authors express their thanks to the Deutsche Forschungsgemeinschaft for financing computer time on the IBM 7094 and the TR 440 of the Deutsches Rechenzentrum Darmstadt and the CD 3300 of the University of Mainz.     

Generalized quantum mechanical two-centre problems

For the one-electron Schrödinger equation among the solutions of which the Slater-Zener-type functions can be found, it is shown, that it can be generalized to the two-centre case only in one way, if one demands separability in prolate spheroidal coordinates, and if in addition to the Coulomb term of the potential energy there shall be an additional function of the product r ₁ · r ₂ only. The generalized problem with a potential energy of the form V(r) = ? Z₁/r₁ ? Z₂/r₂ ? Q(R)/r₁r₂ is studied for the case of two equal centres Z ₁=Z₂=Z≧0 with regard to the existence and number of bound states. The results are extended as far as possible also to the case with unequal centres. For some examples with equal centres wave functions and correlation diagrams have been computed exactly for the lowest electronic states.     

Nonadiabatic transitions in chemical reactions: A quantum mechanical study

This work presents an exact quantum mechanical treatment of a reactive three-atom collinear model system incorporating nonadiabatic couplings. It was assumed that nonadiabatic transitions are induced by the vibrational motion only. The main findings are: (i) The reaction process can create conditions in which weak nonadiabatic couplings terms ( for which the Massey parameter was round 10) may cause large probabilities (～0.5) for transitions from one electronic surface to the other. In other words, the reaction process is able in certain cases to create a near resonance situation which makes the non-adiabatic transition almost independent of the magnitude of the coupling term. For this to happen the two surfaces need not be proximate, nor need they “almost” cross along a certain line (ii) In cases where the main nonadiabatic transitions take place outside the interaction region one may, at least qualitatively, decouple the reaction process from the nonadiabatic one. Thus, under the conditions specified one may first treat the reactive system on the ground state surface without including the excited interacting surface and then treat the nonadiabatic process independently.