QM/MM methods, junction, physical principles, APSLG scheme, effective Hamiltonian
Modern computational quantum chemistry tends to cover with acceptable precision molecular systems of real interest. In the framework of ab initio methodology achieving a good quality results is usually concerned with extending basis sets of one-electron states and with explicit taking into account a great deal of electron correlation. It leads to extremely high computational costs for medium and large size systems due to unpleasant scalability of computational resources like N4¸N7 (where N is the dimension of a space spanned by one-electron basis functions). Even in the case of semiempirical methods the computational costs grow as N3 with growth of system size. This is a serious problem for application of quantum chemical procedures for calculating properties of many practically important molecular systems and especially of their chemical transformations Pople. The problem of high computational costs is especially actual when large systems of biological importance are considered or massive calculations of potential energy surfaces (PESs) are necessary, for example, in the framework of molecular dynamics simulations.
In the literature different means are proposed to significantly reduce computational costs without deteriorating the quality of obtained results. The first way is to smooth the dependence of computational costs on the size of the system. This type approaches usually exploit the localization of electronic degree of freedom, based on the ''principle of nearsightedness'' [], or the exponential decay of the one-electron density matrix elements []. The almost linear scalability is not impossible and, for example, authors of Ref. [] have shown that the scalability of the order N1.3 can be achieved. There is a reasonable number of successful attempts to achieve optimal N-scalability properties for the whole spectrum of quantum chemical schemes - ab initio, DFT, and semiempirical. The non-variational schemes are based on the ''divide-and-conquer'' methodology [,,], Fermi operator expansion [,], energy renormalization group Head-Gordon and recursion technique []. The variational approaches [,,,] are based on the substitution of energy minimization by the grand canonical potential minimization []. These techniques are thoroughly reviewed in Refs. [,]. At the same time these methods are mostly oriented on the tight-binding models of solids or solid clusters. The density matrix renormalization group method [] is quite important achievement since it competes in quality with most elaborated methods of conventional quantum chemistry []. We should stress that the use of local one-electron basis states can significantly reduce computational costs Werner. In this context direct determination of localized Hartree-Fock orbitals can be especially useful and workable []. Significant acceleration of computations can be achieved by using pseudodiagonalization procedures [] or by special choice of the trial electronic wave function [,,,] alternative to the standard SCF form. The application of the methods with good scalabilty properties to enzyme reactions, structure and properties of proteins and DNA is straightforward [,].
The molecular mechanics (MM) [] is another good candidate for using in calculations of properties of large molecular systems. Currently it is extensively used in the field of biochemical simulations. At the same time there are significant limitations for its use. Generally, they are consequences of classical nature of the MM force fields Tchougreeff1999. For example, the application of MM schemes to chemical reactions including the processes of bond formation or cleaveage, transition states, or to highly correlated compounds like transition metal complexes is impossible or, at least, questionable since the main prerequisite for the validity of MM procedures is its application to the ground electronic state of a molecule without low-lying excitations . Moreover, the MM schemes usually require rather complex parameterization procedure. The main advantages of the MM schemes are their low cost and high efficiency in the prediction of molecular geometry for organic compounds without significant electron correlation. The principal advantage of the QM procedures is a wide range of problems they potentially apply. The concert expoitation of advantages of the both (quantum and classical) approaches can be achieved by the incorporation of quantum mechanical (QM) description to the MM framework. This methodology, alternative to construction of the QM schemes with the O(N) scalability, leads to hybrid quantum mechanical/molecular mechanical (QM/MM) schemes allowing another way to bypass the bottleneck of Nn-scalability. The QM/MM schemes describe some relatively small part of the system where chemical transformations take place by an appropriate quantum chemical method while the rest (relatively inert environment) is covered by classical force fields (molecular mechanics). The practical usefulness and general validity of these approaches is based on the chemically motivated observation: chemical transformations usually affect only a small part of the whole system (reaction center) while the rôle of surrounding groups and molecules is reduced to modification of the PES due for example to some polarization or steric strains. This situation is characteristic for chemical reactions of biological interest (especially, for catalysis by enzymes).
Since the seminal work by Warshell and Levitt [] the hybrid QM/MM schemes of calculating large molecular systems acquired an increasing popularity. There is a big variety of different hybrid approaches described in the literature Singh,Field,Robb,Thiel,Hillier,Morokuma,Sauer,Thompson,Rivail,Gao,Sherwood. Evenmore, numerous cases of separating electronic variables like p-electron models [] or even taking into account explicitly only valence electrons in semiempirical methods [] can be considered as special cases of hybrid schemes since they also bear the family marks of the QM/MM approach, namely, (i) the separation of the system into parts, and (ii) treating these parts on quantum or classical levels, respectively. Very close to the QM/MM modeling is a problem of constructing effective fragment (group) potentials which can affect the quantum chemically described part of molecule modeling the influence of environment Ohta,Sanz,Colle,Peyerimhoff,Katsuki,Daudey,Jensen,Gordon,Abarenkov. In such a broader sense several other problems in the area of computational chemistry seem to be related to the QM/MM context: these are the problem of embedding in the cluster calculations on solids and their surfaces with special attention to adsorption and catalysis problems; the problem of description of solute/solvent effects for reactions in condensed media; the pseudopotential descriptions of of molecular and atomic electronic structure to list several most important. Usually all these areas are considered separately without an attempt of establishing relations among them. Also a great variety of different specific schemes referred to as ''protocols'' implemented in different packages are normaly considered only from the point of view of their practical feasibility and their fit for particular applied purpose, rather than in a context of their exact placement among other approximate methods and of evaluation of relative precision of that or another approximation.
The QM/MM modeling can be placed in the general context of combination of different level descriptions. We should mention that the hybrid QM/QM schemes are well known in the literature. Some schemes by construction are far enough from the QM/MM construction since the separation on parts is performed not on the geometrical principle but on the principle of necessity of correlated description to some delocalized orbitals. Well known examples of such an approach is the complete active spase self consistent field (CASSCF) scheme [] and GVBCAS method [] combining generalized valence bond and CASSCF approaches which can be considered as three-level scheme. The effective crystal field method developed for calculation of d-d spectra [] and its extension designed to treat catalytical systems [] can be considered as taking an intermediate position in this ierarchy since they use some strictly localized one-electron states (for example, d-orbitals of one transition metal ion) as a high-level subsystem treating it on the full configuration interaction level with delocalization corrections taken into account perturbatively. The QM/QM schemes similar by construction to QM/MM ones (geometric separation on parts) are also quite popular and should be specially mentioned since they have many common problems . For example, the IMOMO method Humbel,Svensson,Froese,Vreven integrating molecular orbital approximations of different level is very similar in the way of construction, main principles and physical contributions taken into account to the QM/MM IMOMM method []. Another example of such a construction is the QM/DIM (diatomics-in-molecules) scheme [] which takes the interactions between subsystems in the first order of perturbation theory. This scheme can be brought to a type of QM/MM one by relevant significant simplifications [,]: by ignoring, first, the differences in potential curves of diatomic fragments due to angular and spin momentum, second, the couplings between the states of the same symmetry, and multicenter electrostatic contributions. The scheme proposed by Náray-Szabó et al [,] is especially interesting here. It uses strictly localized bond orbitals (SLBOs) for the environment and the SCF procedure for the fragment of particular interest. Note that in the case of large systems the calculation of SLBOs can be very time-consuming process, so the transition to the MM scheme for the environment is an important and almost unavoidable step for calculations of biopolymers.
It should be noted that the hybrid quantum/classical schemes apply not only for determination of geometries, energies, and reaction mechanisms. The Monte Carlo [,] and molecular dynamics (MD) Berweger,Roux,Grochowski,Chalmet simulations are quite popular as frameworks for which various QM/MM procedures serve as ''subroutines''. Before employing hybrid schemes the large-scale MD simulations were performed only with low-level approximations for force fields. The use of hybrid schemes extends significantly the scope of their application, improve precision of the results that allows to improve the understanding of statistical properties and dynamical processes in liquids and biopolymers.
The area of the QM/MM modeling is not a steady one since there are many unsolved problems mostly caused by ad hoc character of construction of junction between subsystems. Despite a general claim of being stipulated by ''specific practical needs'' the ad hoc junction constructs present in the literature frequently lead to more or less serious technical (''practical'') problems in construction of hybrid schemes. This point is stressed in Refs. [,]. The most direct application of hybrid schemes is in the case of solute/solvent interactions since the border between quantum and classical parts of the whole system in this case is naturally defined by the division of the whole system into individual molecules or ions. In this case the environment (classical) variables can be chosen in many different ways, ranging from continual models with special properties [,,,,] to the discrete ones explicitly employing the information on the structure of solvent molecules [,,,] through some intermediate schemes adopting advantages of both extrema [,] are also possible. Nevertheless, even in such a transparent case the problems occur when the formation of complexes between solute metal ions and solvent molecules has to be considered. In this case some solvent molecules may have to be considered either as part of QM system when form a close contact with the metal ion or as a part of the MM system when they are absorbed by the solvent bulk. Approaches to constructing relevant force fields in this situation are discussed in Ref. []. Another important area of applications is solid state and surface chemistry Hillier,Sauer,Eichler,Greatbanks with catalysis and adsorption by oxide systems as the most popular objects. In this case the subsystems are very polar and usually charged, so a special attention should be paid to electrostatic interactions between subsystems. The most important problem in the field of embedding is the unphysical polarization of the QM subsystem (cluster) and the related necessity to adjust effective charges on the ions to obtain sensible results. It should be noted that the most difficult and not trivial case of separating quantum and classical subsystems is that of covalently bound systems since the electrons in chemical bonds connecting two regions are highly correlated. At the same time this case seems to be very actual (especially for biological applications) and requires thorough consideration. In this paper we consider the hybrid QM/MM schemes with covalent linkage between regions in more details, give account of the ''state-of-art'' in this field, and show how the junction between quantum and classical subsystems can be constructed in a consistent (not ad hoc) way by deriving it from an underlying QM description with an emphasis on the choice and refinement of one-electron states related to the intersubsystem boundary and the related parameterization. We consider the possibility of subsequent derivation (i.e. proof) of the model is an argument in favor it (in mathematical sense) while the absence of derivation (or its impossibility) as a contra argument.
The paper is organized as follows. In the next Section we formulate the physical principles of constructing hybrid QM/MM schemes and illustrate their application by separation of electronic variables in the framework of the effective Hamiltonian technique and by derivation of renormalization of QM and MM parameters in the framework of the deductive MM formalism. In the third Section the outline of existing approaches to QM/MM intersubsystem junction is given with numerical examples confirming theoretical conclusions and analysis from the general principles given in the second Section is performed. Finally, a brief conclusion is given.
The structure of the area of the QM/MM modeling is obscured by a great deal of the recipes proposed and by the lack of sequential derivations of them. The situation is partly caused by the structure of the MM itself - the MM scheme is not derived but is taken on the ground of some intuitive concepts of what classical terms must enter the energy energy with adjusted parameters for the latter. We extract here a minimal set of the most essential facts (physical principles) which are necessary prerequisites for constructing hybrid schemes:
These points are very general. Now we consider in more details the practical implementation of these principles.
As a first step for derivation of junction we consider a general formulation of separation of electronic variables of quantum subsystem from those describing electrons in the classical (MM) subsystem which provide a consistent form of intersubsystem junction. This separation is based on the Löwdin partition technique [] and the McWeeny group function formalism []. Generally, we construct an effective Hamiltonian for the QM subsystem in the presence of classical environment and the PES of combined system. This strategy was proposed and developed in Refs. [,,]. Practically it is based on the assumption of existence of a QM scheme underlying the MM one.
We consider two subsystems R- (reactive, quantum) and M- (inert,
classical). The electronic Hamiltonian for the whole system is a sum of
subsystem Hamiltonians and of their interaction which is taken to comprise
the terms of two types - the Coulomb Vc(q) and the resonance (electron
transfer) Vr(q) interactions:
| (0) |
| (0) |
| (0) |
The function Eq. (2) is very general and does not correspond to the
assumed wave function of the hybrid QM/MM method. First of all the numbers
of electrons in the subsystems must be fixed to apply computational schemes
to separate parts on the legal ground. Second, we assume that the M-subsystem is treated with use of the MM, i.e. the PES of the M-subsystem is evaluated without explicit invocaion of its wave function. The
parameters of the M-subsystem must be transferable, i.e. applicable
to combination with any R-subsystem and even in the absence of it. For
this purpose we should use the wave function of the ground state of the
effective Hamiltonian for the M-subsystem since it is in a certain sense
close the wave function calculated without any R-subsystem [].
Thus the required wave function is represented by the antisymmetrized
product of electronic wavefunction for the R-subsystem and that of the
ground state for the free M-subsystem:
| (0) |
After the first partition we obtain
| (0) |
| (0) |
| (0) |
| (0) |
| (0) |
| (0) |
The renormalization of the bare R-subsystem Hamiltonian leads to addition
of some interaction terms to the sum of energies of free subsystems. The
PESs for the whole molecule are obtained by formula:
| (0) |
| (0) |
| (0) |
| (0) |
General physical principles of constructing hybrid QM/MM approaches state that the subsequent derivation of junction between quantum and classical subsystems requires a QM wave function underlying the MM description of PESs. This QM method is necessary, for example, to perform averaging of the effective Hamiltonian in Eq. (7). At the same time more important that this method should produce in a consistent manner one-electron states necessary for explicit formation of boundary and its response on the changes in molecular geometry of fragments and/or electronic structure of the R-subsystem.
There are some general criteria for the QM method appropriate for underlying
the MM description: this method should express molecular electronic
structure and electronic energy in relevant local terms (i.e. to
distinguish bonded and nonbonded contributions) and reproduce molecular
properties with sufficient accuracy. The parameters characterizing the wave
function of this method (electronic structure parameters or ESPs) have to be
transferable in a broad sense of the term ''transferability'', i.e.
the form of any bond-related functions (e.g. the bond energy
dependence on interatomic separation) must be also transferable. As an
appropriate method satisfying these conditions we take the trial wave
function in the form of antisymmetrized product of strictly localized
geminals (APSLG) [] implemented with slightly modified
semiempirical MINDO/3 Hamiltonian [,]. The APSLG wave
function is constructed from two-electron building blocks - geminals:
| (0) |
| (0) |
The important feature of the APSLG-MINDO/3 energy is that it is a sum of
increments of different form for bonded and nonbonded atoms:
| (0) |
| (0) |
The results of Refs. [,] related to heats of formation and equilibrium geometries obtained by the APSLG-MINDO/3 method are somewhat more accurate than those of the standard (SCF-)MINDO/3 method. The APSLG form of the trial wave function also ensures its correct asymptotic behaviour under cleavage of chemical bonds which indirectly justifies some level of bonae fidelitatis of the wave function employed. The APSLG form of the wave function also allows to represent a renormalization of the bare R-subsystem Hamiltonian in terms of well-defined characteristics like atomic charges, bond polarizabilities and ionization potentials for the chemical bonds [,].
In the framework of QM/MM junction construction we should determine the structure of boundary one-electron states and their responses to variations of geometric parameters. The problem of constructing optimal boundary HOs is closely related to the more general problem of deduction of the MM description from some consistent QM description of electronic structure. This problem is actual since the form of the force fields and their particular sets of parameters are not justified that leads to a great variety of different MM schemes. The semiempirical APSLG-MINDO/3 method briefly described above is a perfect candidate for the MM derivation and it was successfully used for this purpose in Refs. [,]. Here we describe the main steps of the MM derivation and their consequences for the QM/MM modeling.
To derive the MM picture we must consider the question of transferability of
the ESPs. Above we defined two groups of parameters entering the expression
for the energy Eq. (17). Let us consider the first group of
parameters. In the case of lone pairs they are perfectly transferable. The
SCF approximation for non-polar bonds gives the geometry independent density
matrix elements:
| (0) |
| (0) |
| (0) |
The specific bond equilibrium distance can be obtained as the minimum
position for the bond-related energy. The core-core interaction in the
MINDO/3 scheme is not purely Coulomb one but is modified by a short-range
repulsion term:
| (0) |
| (0) |
| (0) |
General construction of the MM scheme and analysis of adjustment of
one-electron states on the QM/MM boundary can be performed using linear
response relations between the geometry parameters and the ESPs. These
relations can be obtained by expanding the energy up to the second order
with respect to both atomic coordinates q and the ESPs x with subsequent
minimization:
| (0) |
| (0) |
Mathematically the hybridization of the sp-basis at a given atom is
defined by a SO(4) rotation. Each HO is thus a normalized quaternion with
a scalar and a 3-vector parts (sm,[(v)\vec]m). The set of four vector parts
of HOs forms a hybridization tetrahedron with vectors [(v)\vec]m pointing to
its vertices. In principle, deductive MM can be constructed as a set of
classical type equations for these tetrahedra []. Here we
consider, however, more common case of the MM force fields. The dynamical SO(4) group can be parameterized by a pair of normalized quaternions
corresponding to fictitious rotations combining the hybridization (s-p
mixing, pseudorotation) and rotation of hybridization tetrahedron as a whole
(quasirotation). At the same time the first order variation of the HOs when
small quasi- and pseudorotations d[(w)\vec] l and d[(w)\vec] b are applied to the system of HOs can be derived in a very simple
form:
| (0) |
Let us demonstrate an application of formulae Eqs. (25) and (
26) on the example of methane molecule with sp3 hybridization. The
stereochemistry of carbon atom is characterized by four unit vectors [(e)\vec]RmLm corresponding to the directions of the CH bonds. In the diatomic
coordinate frame with the z axis directed along the [(e)\vec]RmLm
vector the resonance integrals are:
| (0) |
The equilibrium conditions for the form and orientation of hybridization
tetrahedron can be obtained by zeroing the pseudotorque [(N)\vec] and
quasitorque [(K)\vec] which are coefficients at d[(w)\vec] b and d[(w)\vec] l in the resonance energy expansion:
| (0) |
The fixation of parameters [(w)\vec] b (those of the shape of the
hybridization tetrahedron) can be considered as FO (fixed orbitals)
picture []. In this case only the resonance energy is the
orientation dependent contribution. The angular dependence of the energy
(bending) can be described by introducing small rotation vectors d[(j)\vec] m, which after applying them to vectors [(e)\vec]RmLm
lead to new (distorted) coordination tetrahedron. The final result of
substitution of second derivatives of the energy to Eq. (26) is
the explicit expression for the bending force field constant which can be
written as:
| (0) |
In general both orientation and form of hybridization tetrahedron must be
adjusted. This picture can be termed as TO (tuned orbitals)
GenMM. The reaction of hybridization tetrahedron on the changes of local
geometry can be considered in the linear response approximation Eq. (
25). In the case of valence angles deformation the algebraic structure
of the SO(4) manifold ensures quite special structure of geometry
deformations which can cause changes in the form of hybridization
tetrahedron. These deformations (hybridization compatible) correspond to
increasing of a valence angle by certain amount dq1m with
simultaneous decreasing of the opposite (spiro) angle by the same magnitude.
Substitution of matrices of second derivatives for the energy into Eq. (
25) gives the reaction of the form of hybridization tetrahedron on the
angular distortions of molecular geometry in the form:
| (0) |
| (0) |
| (0) |
| (0) |
The above results constitute a basis for construction of intersubsystem
junction in hybrid QM/MM schemes. There are two different effects arising
from separation of system on parts and consideration of one part as a MM one
- renormalization of the QM Hamiltonian parameters and imposing additional
forces and torqes on the MM subsystem. Let us exemplify the technique of
short-range intersubsystem junction construction by one special but quite
characteristic case of a boundary sp3 carbon atom with one HO
pointing to the QM region (bond m = 1) and others related to the MM one
(this case can be classified as the MM boundary atom). The transition to the
MM picture is performed by setting the FA and TO approximations. To
determine the effect of the QM part on the MM part we consider the
perturbation of molecular energy due to changes of electron densities in the
QM region. It affects both one-center energy of boundary atom and resonance
energy between HO r1 of the QM bond with number m = 1 and all other HOs
in the QM subsystem. The perturbation sets up quasi- and pseudotorques for
the boundary atom:
| (0) |
| (0) |
The additional pseudo- and quasitorques produce the pseudo- and
quasirotations of the hybridization tetrahedron on the boundary atom R1.
In the linear response approximation it corresponds to the treatment of the
corresponding pseudo- and quasitorques by the ( Ñ[(w)\vec]2E) -1 matrix which is simple (diagonal in tha basis of d[(w)\vec] b and d[(w)\vec] l) in the case of symmetric
hydride:
| (0) |
| (0) |
| (0) |
| (0) |
| (0) |
It should be noted that the MM subsystem also affects the parameters of the
QM subsystem since any geometry deformation in the MM subsystem induces
changes in pseudo- and quasirotation angles defining hybridization of the
boundary atom. It renormalizes the QM one-ceneter Hamiltonian parameters
according to:
| (0) |
| (0) |
To illustrate the formulae given above we performed the estimates of the
magnitude of the renormalization of the QM parameters and changes in the QM
forces and torques due to boundary taken as sp3 carbon atom. The changes in
the QM one-center Hamiltonian parameters due to elongation of one of the
boundary MM bonds are:
| (0) |
| (0) |
| (0) |
In the present Section we employ the theoretical framework developed above to rationale the state of art in the field of hybrid QM/MM modeling. It is instructive to invoke different classifications of QM/MM schemes developed. The simplest but not very much informative classification can be based on the types of quantum mechanical and molecular mechanical schemes used. There are no fundamental restrictions on the choice of the QM schemes and in fact very large diversity can be found in the literature. All types of QM procedures - ab initio Rivail,Truong,Eichler,Freindorf,Luque,Mo,Yoshida, DFT Tunon,Wei,Sierka,Stanton,Lyne,Long,French, and semiempirical Thiel,Cummins,Thery,Ferenczy - are quite popular and widely used as well as in usual quantum chemical studies. The usual form of wave function implied is the SCF one. At the same time methods based on the valence bond approximation are also well known [,,]. The choice of the MM scheme can be quite important since it affects the structure of bonding terms near (or in some schemes on) the border between the QM and MM parts and the electrostatic polarization of the QM part by the MM-treated environment environment depends on the particular charge scheme employed in the MM procedure. Practically, it is more convenient to work with the force field with electrostatics based on the atomic charges than with the force field using bond dipole interactions for this purpose since it can cause significant errors in the description of polar species. This notion have led the authors of Ref. [] employing the MM3 force field [] to replacing the the bond dipoles with potential-derived atomic point charges. We note that the sequential derivation of junction in the framework of the effective Hamiltonian approach reported in leads to representing electrostatic interactions in terms of atomic charges []. If the force field contains many explicit couplings between bonded interactions (like the MM3 force field) it may cause additional problems with construction of junction. Some words should be said about ionic force fields using formal ionic charges and employing electrostatic and short-range force fields. In the case of these force fields short-range interactions arising from the MM charges can not be separated from the long-range interactions []. These schemes often lead to incorrect electrostatic potentials so the most popular choice (especially in description of covalently bound QM/MM systems) is the valence force fields. Different implementations use different force fields: for example, the MM3 force field [] is used in Refs. [,], the CHARMM force field [] is used in Refs. [,], the AMBER force field [] is used in Refs. [,].
In the previous Section we obtained the formula for junction between quantum
and classical subsystems Eq. (13). The control for the types of
interactions which are taken into account is an important characteristic of
particular QM/MM scheme. The authors of Ref. [] have proposed a
classification of hybrid schemes based on the interaction between fragments.
According to it, the simplest type of model is mechanical embedding
(examples of this type of modeling are the IMOMM [] and IMOMO
[] schemes by Morokuma) when both QM and MM systems are not
polarized by each other and their interaction is represented by classical
force fields only. In this context the choice of parameters of intersystem
interaction can be crucially important, so, they are frequently optimized
[,]. More elaborated model is that including polarization
of the QM subsystem. This polarization can be covered by including the MM
charges into one-electron part of the Hamiltonian of the QM subsystem:
| (0) |
It should be noted that the classification given by authors of Ref. Bakowies is not complete since, for example, it does not include some self-consistent schemes like that of Ref. [] and does not consider the possibility of charge transfer between subsystems TokmTchMis,EnQMMM. The authors of Ref. [] have imposed a requirement of intersubsystem self-consistency on the construction of junction. It means that the charge transfer between subsystems should be taken into account. Practical implementation of this requirement was performed by using special iterative procedure of double (intrafragment and interfragment) self-consistent (DSC) calculations. It leads to explicit taking into account the electron transfer between subsystems (and also the polarization of the QM subsystem). This methodology, however, can not be justified since the self-consistent field procedures are applied to systems with strongly fluctuating numbers of electrons that also leads to poor definition of the fragments themselves (and of their quantum numbers). According to results of the previous Section the electron transfers should be considered as virtual ones and taken into account in the perturbative fashion. This point is confirmed numerically since the application of the DSC scheme to the iron picket-fence porphyrin has led to improbably large intersubsystem charge transfer of 3.6 electrons [].
In the case of intersubsystem junction represented by classical bonding terms an important question arises: which terms should be included and which should not? The most popular way is to include classical bonded force fields when at least one MM atom is involved in it [,]. At the same time it does not allow to avoid double counting of interactions computed quantum mechanically. To smooth this inconsistence the authors of Ref. Eurenius proposed to calculate only those classical bonding force fields where at least one central atom is from the MM subsystem or in the case of improper dihedral terms only those with both outer atoms from the MM subsystem.
The ways to construct the boundary between subsystems and general ideas concering the structure of the ''grey'' zone between quantum and classical fragments proposed in the literature are miscellaneous. For example, authors of Ref. [] pose the requirement of so called ''QM-MM continuity'' that leads to introducing intermediate fragments which are treated by both QM and MM methods. Since the early times of quantum chemistry and until now it was popular to make the problem of molecular structure investigation tractable by neglecting large substituents and saturating free valencies by hydrogen atoms. The QM/MM methodology tries to take the bulky substituents into account explicitly. The most straightforward way to treat covalently bound QM and MM parts is the link atom method stemming from saturating dangling bonds produced while cutting the system in parts across the bonds by additional atoms (usually, hydrogens). This methodology requires deleting the terms in the expression for the total energy of the system related to the link atoms. Practically it is very difficult since in the framework of this scheme serious artefacts appear quite expectedly. It should be noted that the polarization of the MM subsystem is often reproduced incorrectly. To cure this error some interaction terms for atoms near the boundary are typically artificially eliminated. Originally [], all the interactions (electrostatic and van-der-Waals) between link atoms and MM subsystem were not taken into account. Later, different recipes to omit some real physical interactions and by this compensate unphysical ones were proposed: the authors of Ref. [] neglect the Coulomb interactions between the QM subsystem and the MM group most close to it, the authors of Ref. [] force the charge on the boundary MM atom to be zero. It should be noted that such eliminations can not be justified and must be considered as special methods for masking the errors caused by the inconsequence in junction construction. In fact, the problem is often a consequence of the small distance between the link atom and the first MM atom. The more complex scheme [] does not neglect the interactions of link atom with the MM subsystem but sets charge on the boundary MM atom to be equal to zero. Another recipe used is shifting the values of the charges with subsequent compensation of this perturbation by introducing fictitious dipoles []. Typically, the possibility to manipulate by interactions of the artificial link atom are considered as an advantage of special flexibility of the link atom approach []. However, to our opinion the price of this flexibility is too high since it leads to a great uncertainty in the results obtained and thus marks down possible predictive capacity of the QM/MM approach.
The saturation of dangling bonds by hydrogen is not the only possible way proposed in the literature. The main reason of using other types of saturation groups is the intention to reproduce the polarity and other properties of the broken bond more correctly. The choices (suggestions) for saturation groups well known in the literature are pseudohalogens HyperChem,CumGr and ''dummy groups'' []. The important drawback of link atom method (and, especially of the ''dummy groups'' method) is introducing of additional nuclear degrees of freedom for which no reasonable equation of motion (or equivalently no energy minimum condition) can be derived. It requires special link atom corrections. To cope with problems incurred by introducing link atoms an adjusted connection atom approach [] was developed. The connection atom is a QM one but its interactions with the environment include some classical terms. So, the bonded interactions of the connection atom with MM atoms are covered by force field for a carbon atom. It is assumed that the C-C bond is cut and the connection atom mimics the methyl group. This method was calibrated against the NDDO calculations for methyl compounds. It should be noted that this scheme is valid only for equilibrium bond distances since the molecules used for calibration were taken in their equilibrium geometries.
The same problems can be addressed by using effective group potential techniques [,]. For example, in Ref. [] a special potential reproduces the interactions between the QM subsystem electrons and the missing valence electrons of one-electron boundary atom. In the pseudobond approach [,] the free valence of a QM atom is saturated with a special atom located exactly at the position of the neighbour MM atom. The basis set and number of electrons of this pseudoatom are set to be equal to those of fluorine atom. The C-C bond is mimicked by a specially adjusted effective core potential. Another approach based on the use of effective potentials is proposed in Ref. []. A series of one-electron quantum capping potentials replacing the link atom is formed by modification of conventional effective core potentials: the spherical shielding and Pauli terms replace the excess valence electrons. The capping potential is adjusted to reproduce all-electron geometries and charge distributions. The analysis of this scheme shows that the error induced by capping potential is significant (especially for angular distortions) but generally less than that in the simple link atom scheme. It is noteworthy that this scheme can be extended to construct many-electron quantum capping potentials.
The important problem of the link atom schemes is geometry optimization
since as it is shown in Ref. [] the difference between QM and
MM interaction energies for link atom enters into the total energy and
strongly depends on the link atom position since the forms of the QM and MM
force fields are very much different. Practically it leads to the collapse
of the fictitious link atom with the boundary atom. The characteristic
results is that the equilibrium position of the link atom is badly defined
and can not be reasonably rationalized. One of possible recipes for avoiding
the collapse of the link atom is to use for geometry optimization a
potential energy function not coinciding with the total energy of the
molecular system []:
| (0) |
It is seen that the implementation of the link atom method is generally not
a trivial task. The most simple for implementation model closely related to
the link atom method [] is the so-called ''subtractive scheme''
implemented for example in Refs. [,]. In the framework of
these schemes the dangling bonds of the QM fragment are sauturated by some
groups (hydrogens, methyl groups etc.) and form the so-called model
system. The whole system is calculated by a MM method of choice (or, more
generally, by whatever low-level method, may be a QM one) and the energy is
obtained by adding the QM (high-level) calculated energy of the model
subsystem and subtracting the MM calculated energy of the model subsystem:
| (0) |
| (0) |
The main features of the ''subtractive'' hybrid schemes can be illustrated by numerical examples. To understand better the problems caused by the ad hoc way of junction construction we consider characteristic examples where these schemes fail. In the most cases the sources of failures can be seen on the simple molecular objects which serve as tests for one or another hybrid scheme. As an example we consider an application of the IMOMM method to conformational properties of cis-butane performed in Ref. []. It was shown that marking two C end atoms of the molecule as QM ones leads to valence angles 129.9°, while pure QM calculation gives 117° and pure MM calculation predicts the value 116.1° for these angles. Therefore, we see the case when transition to QM/MM procedure spoils even the results of pure MM calculations. Problems also arise when multilayered schemes are used. For example, the energy of reaction of the oxidative addition of H2 to Pt(P(t-Bu)3)2 calculated by the B3LYP:HF:MM3 scheme is 7.9 kcal/mol smaller than that calculated by the B3LYP:HF:HF scheme [50]. The same difference in 7.9 kcal/mol is obtained for schemes B3LYP:B3LYP:MM3 and B3LYP:B3LYP:HF. Therefore, the choice of the description for the third layer (the most inert subsystem) turns out to be crucial for the description of the energy of this reaction which is estimated to be ~ 4 kcal/mol [50]. It speaks on the incorrect construction of junction in this scheme.
Principally different and more justified type of construction of junction
between QM and MM parts of the molecule is based on the use of local
one-electron functions. Only this type junctions cutting the nuclei and not
the bonds can be justified from the general point of view. In the local
self-consistent field (LSCF) approach [,,,]
the chemical bonds between QM and MM regions are represented by SLBOs. The
BOs can be obtained by whatever of a posteriori localization
procedures known in the literature. At the same time the localized orbitals
obtained in the localization procedures have some degree of delocalization,
i.e. they have non-zero contributions for the atoms not constituting a
given bond (or a lone pair) mimicked by this particular BO. These
contributions are called the ''tails'' of the localized orbitals and their
neglect lead to the SLBOs which are used in the LSCF scheme. An important
assumption made in the LSCF construction is that of the transferability of
the SLBOs in very wide ranges. The QM part of the system is described by a
set of delocalized MOs while the boundary is modeled by the frozen SLBOs.
The frozen character of the boundary SLBO may cause the sensitivity of the
LSCF results to the size of the QM region. The determination of the
electronic structure of the QM region is based on the diagonalization of the
modified Fock operator which includes Coulomb and exchange interaction with
boundary orbitals and interaction with the MM charges. In its original
implementation the LSCF method was developed by fixing the positions of the
atoms of the environment and thus it was not suitable for geometry
optimization procedures. This restriction was avoided in Ref. []
by special adjustment of the force field parameters. The authors of Ref.
[] have noticed that in the framework of original LSCF scheme the
ion-nuclei interactions are underestimated and the variation of the overlap
between boundary basis functions with respect to the boundary bond length is
not taken into account. To cure these problems they defined a 5-components
boundary bond potential
| (0) |
It is interesting to compare the possibilities and errors of different hybrid QM/MM schemes. The careful examination and comparison of link atom and LSCF techniques was performed in Ref. [] using the CHARMM force field [] and the AM1 method [] as a quantum chemical procedure. In the case of the link atom procedure two options were used: QQ - the link atom does not interact with the MM subsystem and HQ - link atom interacts with all MM atoms. The main conclusion of this consideration is that the LSCF and the link atom schemes are of similar quality. The error in the proton affinity determination induced by these schemes is several kcal/mol. It is noteworthy that all the schemes work rather badly in description of conformational properties of n-butane. The large charge on the MM atoms in the proximity of the QM subsystem (especially on the boundary atom) cause significant errors in the proton affinity estimates for all methods (especially, in the case of the LSCF approach where the error can be of tens of kcal/mol). This is not surprising since the stability and transferability of intrabond one- and two-electron density matrix elements Eq. (19) is broken here. It proves that the simple electrostatic model is not well appropriate for these schemes and that a detailed analysis of the junction form is necessary in general case. Moreover, the numerical analysis shows that the error induced by the HQ model is less than that induced by the QQ model. Since the HQ model explicitly includes unphysical interactions with artificial link atom it means that these interactions partly compensate errors in junction construction. Practically, the link atom interacts with the MM atoms even in the case of the QQ model but indirectly via the interactions of QM subsystem with the MM atoms. In the case of the QQ model non-compensated charge on the QM subsystem (without link atom) interacts with MM atomic charges. It causes significant polarization of the QM subsystem. It is confirmed by numerical estimates of charges on the QM subsystem [].
An attempt to escape the need of constructing large databases in the LSCF approach was made by Gao et al in the generalized hybrid orbital (GHO) method [,]. This approach is intended to interpolate the ESPs related to the shape and orientation of the HOs residing on the frontier atoms. The first step in the GHO method is dividing the hybrid orbitals into active and auxiliary ones - the first ones are added to the QM subsystem. In the original LSCF approach three orbitals of the boundary atom are included into the self-consistent procedure. In the GHO approach only one active orbital of the boundary atom participates in the set of orbitals subjected to the SCF procedure while the auxiliary orbitals mimic the effective core potential. Therefore, we can characterize the boundary atom in the GHO scheme as the MM one while the boundary atom in the LSCF scheme is the QM one. It allows to substitute the parameterization of the one-electron density on the HO for each molecule by choosing semiempirical parameters for boundary atom with assumption of their transferability. At the same time the possibility to choose reasonable and transferable semiempirical parameters for boundary atom must rely upon correct determination of the HOs. Conversely, the GHO scheme uses very crude assumptions [] about the form and the direction of HOs based on the pure geometry considerations which are equivalent to assumptions of (i) C3v symmetry of the local MM environment and (ii) coincidence of all the HOs directions with the directions of bonds. Practically, these conditions are satisfied only for methane molecule. Moreover, even for these assumptions the s-/p-ratio for active orbital is determined by incorrect formula working only for equivalent active and auxiliary orbitals or for purely p active orbital. The real structure of the HOs as a function of the geometry is given by the expressions of the previous Section and we note that neither of the assumptions accepted in Ref. [] is fulfilled. Neither asymmetry of the geometry of the boundary atom environment nor the chemical nature of neighbouring groups are taken into account in Ref. GHOOpt. In practical implementation of Refs. [,] all the HOs are determined by directions of the bonds between the boundary atom and its three MM neighbours. The coincidence of HOs and bond directions is a rather crude and uncontrollable assumption. In the case of molecules with large asymmetry or significant geometry constraints the difference between directions of bonds and HOs can exceed 20° with cyclopropane as the most characteristic example. Since the s-/p- ratio and direction of the active HO in the GHO method are totally the functions of only MM atoms positions the form of this HO may be far from being optimal to say nothing of its behaviour with geometry variations. In our opinion the non-variational form of the HOs has led the authors of Ref. [] to significant and hardly justifiable renormalization of the Hamiltonian parameters and the MM force field to reproduce correct bond lengths, direction of auxiliary orbitals along the corresponding bonds, and the Mulliken charges on the atoms, so, that the boundary could not be considered as a weak perturbation anymore. For example, in the case of carbon atom the resonance parameter bs should be changed by more than 10 eV; the MM C-C ideal bond length parameter r0 should be changed by 0.05 Å . It should be noted that the change of the delocalized description of electronic structure to the localized one does not lead to significant changes in the parameters of the semiempirical Hamiltonian as it was exemplified on the examples of MINDO/3 [], MNDO, AM1, and PM3 [] schemes.
The principles similar to those of the LSCF are used for junction construction in Ref. [] based on the fragment SCF method. Another model thoroughly elaborated to consider effect of motions of environment atoms on the ab initio level is that by Philipp and Friesner []. In the framework of this model the procedure of frozing the SLBO was refined and the degree of delocalization of SLBOs is extended from two boundary atoms to two boundary groups. At the same time an attempt to obtain numerical results of good quality has led authors to introducing some very artificial procedures in the construction of the junction like placement of very large charges on the bond. Moreover, this model requires extensive parameterization. To reproduce energetics of alanine dipeptide and tetrapeptide 27 parameters describing interface between subsystems were adjusted. Essential parameterization process is quite important problem of this model. It would be fair to say that reproduction of conformational energies for a polar system like polypeptide is a difficult test for any computational methodology. The authors of Ref. [] give the essential requirements for the bond at which molecule is cut: (i) this bond should not have significant multiple bonding character, and (ii) this bond should be far away from the region where significant electron re-distribution occurs. The first requirement seems to be quite reasonable for description of bond by single localized orbital, while the second one reflects a lack of some contributions to the energy which are important for small QM subsystems chosen. The estimate of these contributions is not possible due to ad hoc way of construction of intersubsystem junction.
The method based on the effective fragment potential (EFP) construction implemented in Ref. [] is close to the general LSCF methodology. In this case the boundary is modeled by a buffer region consisting of several localized molecular orbitals which are obtained by a QM calculation on all or a subset of the system. These orbitals are set frozen in the EFP calculation. The orbitals of the QM part are forced to be orthogonal to those of the buffer region and the environment is represented by an EFP. The important idea of this approach is to make the distance between the QM and EFP regions significant with the purpose to present the corresponding intersubsystem interactions as nonbonded ones. At the same time the frozen character of the buffer one-electron states can be important since the changes in polarization contributions from the buffer region are neglected during a geometry optimization. The numerical example given by authors of Ref. [] is quite characteristic. They calculated the proton affinities of lysine and H-bonded and non-H-bonded tripeptide Gly-Lys-Gly by the QM/buffer/EFP method. If the buffer region is chosen as g- and d-CH2 groups of lysine chain (i.e., quite far from the reaction center) the QM/buffer/EFP calculation gives the value of the proton affinity 2.2 kcal/mol higher than the reference QM one for all these molecules. It unequivocally testifies that these 2.2 kcal/mol constitute the error of junction construction in this case which seems to be quite large. Moreover, the difference between QM/buffer and QM/buffer/EFP results for proton affinities of lysine and non-H-bonded tripeptide Gly-Lys-Gly is only 0.2 kcal, i.e. the effect of environment described by the EFP is by an order of magnitude smaller than the error produced by the junction.
The hybrid QM/MM modeling is quite promising for study of large molecules especially in the rapidly growing field of computational biological chemistry. The case of covalent linkage between the QM and MM parts is especially important since it includes the enzymatic catalysis. At the same time this case is the most complex one since the boundary between subsystems is not well defined and construction of intersubsystem junction is not straightforward. In this case many ad hoc recipes are proposed in the literature. We have shown the state-of-the-art in this field in somewhat critical manner with special attention to the problems arising during the QM/MM modeling. At the same time the sequential derivations of QM/MM junctions are shown to be possible. We have given a series of physical principles which in our opinion should govern the construction of hybrid QM/MM schemes. As particular implementations of these principles we have demonstrated how the effective electronic Hamiltonian approach can be applied to the purpose of junction construction. It has required the use of quantum chemical description underlying the MM one. The APSLG-type trial wave function was taken for this description.
It should be noted that the important question of structure of boundary one-electron states (especially as a function of molecular geometry) remained unclear. We proposed a way of determination of these states based on the derivation of the MM description from the QM (APSLG) one. The derivation of the deductive MM allows to determine the effects of the MM subsystem on the QM one (renormalization of parameters) and of the QM subsystem on the MM one (torques and forces acting on the MM atoms). Explicit expressions are given for example of boundary sp3 carbon atom. Numerical estimates obtained for different QM/MM schemes illustrate the general points stressed in the text.
Financial support for AMT on the part of the Haldor Topsoe A/S is gratefully acknowledged. Authors are grateful to Professor J.-L. Rivail for the reprints and preprints of his works on the hybrid QM/MM methods, to Prof. I.V. Abarenkov for permission to use his work prior of publication and to Dr. X. Assfeld and Dr. A.A. Sokol for valuable discussions.
1
J. A. Pople Quantum chemical models (Nobel lecture), Angew. Chem. Int. Ed. 38, 1894-1902 (1999).
W. Kohn Density functional and density matrix method scaling linearly with the number of atoms, Phys. Rev. Lett. 76, 3168-3171 (1996).
S. Goedecker Decay properties of the finite-temperature density matrix in metals, Phys. Rev. B 58, 3501-3502 (1998).
K. N. Kudin and G. Scuseria A fast multipole method for periodic systems with arbitrary unit cell geometries, Chem. Phys. Lett. 283, 61-68 (1998).
W. Yang Direct calculation of electron density in density-functional theory, Phys. Rev. Lett. 66, 1438-1441 (1991).
S. L. Dixon, K. M. Merz Jr. Fast, accurate semiempirical molecular orbital calculations for macromolecules, J. Chem. Phys. 107, 879-893 (1997).
J. J. P. Stewart Application of localized molecular orbitals to the solution of semiempirical self-consistent field equations, Int. J. Quantum Chem. 58, 133-146 (1996).
S. Goedecker and L. Colombo Efficient linear scaling algorithm for tight-binding molecular dynamics, Phys. Rev. Lett. 73, 122-125 (1994).
S. Goedecker and M. Teter Tight-binding electronic-structure calculations and tight-binding molecular dynamics with localized orbitals, Phys. Rev. B 51, 9455-9464 (1995).
R. Baer and M. Head-Gordon Electronic structure of large systems: Coping with small gaps using the energy renormalization group method, J. Chem. Phys. 109, 10159-10168 (1998).
R. Haydock, V. Heine and M. J. Kelly Electronic structure based on the local atomic environment for tight-binding bands, J. Phys. C 8, 2591-2605 (1975).
W. Yang Absolute-energy-minimum principles for linear-scaling electronic-structure calculations, Phys. Rev. B 56, 9294-9297 (1997).
F. Mauri, G. Galli and R. Car Orbital formulation for electronic structure calculations with linear system size scaling, Phys. Rev. B 47, 9973-9976 (1993).
P. Ordejón, D. Drabold, M. Grumbach and R. M. Martin Unconstrained minimization approach for electronic computations that scales linearly with system size, Phys. Rev. B 48, 14646-14649 (1993).
M. Challacombe A simplified density matrix minimization for linear scaling self-consistent field theory, J. Chem. Phys. 110, 2332-2342 (1999).
X.-P. Li, R. W. Nunes and D. Vanderbilt Density matrix electronic structure method with linear system size scaling, Phys. Rev. B 47, 10891-10894 (1993).
S. Goedecker Linear scaling electronic structure methods, Rev. Mod. Phys. 71, 1085-1123 (1999).
S. Y. Wu and C. S. Jayanthi Order-N methodologies and their applications, Physics Reports 358, 1-74 (2002).
S. R. White Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69, 2863-2866 (1992).
S. R. White and R. L. Martin Ab initio quantum chemistry using the density matrix renormalization group, J. Chem. Phys. 110, 4127-4130 (1999).
M. Schütz, M. Hetzer and H.-I. Werner Low-order scaling local electron correlation methods. I. Linear scaling local MP2, J. Chem. Phys. 111, 5691-5705 (1999).
J. Rubio, A. Povill, J. P. Malrieu and P. Reinhardt Direct determination of localized Hartree-Fock orbitals as a step toward N scaling procedures, J. Chem. Phys. 107, 10044-10050 (1997).
M. C. Strain, G. E. Scuseria and M. J. Frisch Achieving linear scaling for the electronic quantum Coulomb problem, Science 271, 51-53 (1996).
A. Derecskei-Kovacs, D. E. Woon and D. S. Marynick Nonempirical wave functions for very large molecules. II. The PRDDO/M/FCP method, Int. J. Quantum Chem. 61: 67-76 (1997).
J. M. Cullen A rapid generalized valence-bond algorithm for semiempirical NDDO calculations, Int. J. Quantum Chem. 56, 97-113 (1995).
A. M. Tokmachev and A. L. Tchougréeff Semiempirical implementation of strictly localized geminals for analysis of molecular electronic structure, J. Comp. Chem. 22, 752-764 (2001).
A. M. Tokmachev and A. L. Tchougréeff Fast NDDO method for molecular structure calculations based on strictly localized geminals, J. Phys. Chem. A, accepted.
D. York, T. S. Lee and W. Yang Quantum mechanical study of aqueous polarization effects on biological macromolecules, J. Am. Chem. Soc. 118, 10940-10941 (1996).
J. P. Lewis, C. W. Carter Jr., J. Herman, W. Pan, T. S. Lee and W. Yang Active species for the ground state complex of cytidine deaminase: A linear scaling quantum mechanical investigation, J. Am. Chem. Soc. 120, 5407-5410 (1998).
U. Burkert and N. Allinger, Molecular Mechanics (ACS, Washington, DC, 1982).
M. B. Darkhovskii and A. L. Tchougréeff On unification of quantum chemistry and molecular mechanics. Potential energy surfaces of transition metal complexes as exemplified by spin transition in cis- [Fe(NCS)2(bipy)2], Khim. Fiz. 18, 73-79 (1999) [Chem. Phys. Reports 18, 149 (1999)].
A. Warshel and M. Levitt Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme, J. Mol. Biol. 103, 227-249 (1976).
U. C. Singh and P. A. Kollman A combined ab initio quantum mechanical and molecular mechanical method for carrying out simulations on complex molecular systems: Application to the CH3Cl+Cl- exchange reaction and gas-phase protonation of polyethers, J. Comp. Chem. 7, 718-730 (1986).
M. J. Field, P. A. Bash and M. Karplus A combined quantum mechanical and molecular mechanical potential for molecular dynamics simulation, J. Comp. Chem. 11, 700-733 (1990).
F. Bernardi, M. Olivucci and M. A. Robb Simulation of MC-SCF results on covalent organic multi-bond reactions: Molecular mechanics with valence bond (MM-VB), J. Am. Chem. Soc. 114, 1606-1616 (1992).
D. Bakowies and W. Thiel Semiempirical treatment of electrostatic potentials and partial charges in combined quantum mechanical and molecular mechanical approaches, J. Comp. Chem. 17, 87-108 (1996).
I. H. Hillier Chemical reactivity studied by hybrid QM/MM methods, J. Mol. Struct. (Theochem) 463, 45-52 (1999).
F. Maseras and K. Morokuma IMOMM: A new integrated ab initio + molecular mechanics geometry optimization scheme of equilibrium structures and transition states, J. Comp. Chem. 16, 1170-1179 (1995).
J. Sauer and M. Sierka Combining quantum mechanics and interatomic potential functions in ab initio studies of extended systems, J. Comp. Chem. 21, 1470-1493 (2000).
M. A. Thompson QM/MMpol: A consistent model for solute/solvent polarization. Application to the aqueous solvation and spectroscopy of formaldehyde, acetaldehyde, and acetone, J. Phys. Chem. 100, 14492-14507 (1996).
X. Assfeld and J. L. Rivail Quantum chemical computations on parts of large molecules: The ab initio local self consistent field method, Chem. Phys. Lett. 263, 100-106 (1996).
J. Gao ''Methods and Applications of Combined Quantum Mechanical and Molecular Mechanical Potentials.'' In Reviews in Computational Chemistry (VCH Publishers, Inc., New York, 1996).
P. Sherwood, ''Hybrid Quantum Mechanics/Molecular Mechanics Approaches.'' In Modern Methods and Algorithms of Quantum Chemistry, Proceedings, Second edition, J. Grotendorst, ed., NIC Series, 3 (John von Neumann Institute for Computing, Jülich, 2000).
T. E. Peacock, Electronic Properties of Aromatic and Heterocyclic Molecules (Academic Press, London, 1965).
K. F. Freed, ''Theoretical Basis for Semiempirical Theories.'' In Semiempirical Methods of Electronic Structure Calculations. Part A. Techniques. G.A. Segal, ed. (Plenum Press, New York, 1977).
K. Ohta, Y. Yoshioka, K. Morokuma and K. Kitaura The effective fragment potential method. An approximate ab initio MO method for large molecules, Chem. Phys. Lett. 101, 12-17 (1983).
J. A. Mejias and J. F. Sanz Unrestricted compact model potentials for ab initio embedded cluster calculations: Magnetic interactions in KNiF3, J. Chem. Phys. 102, 850-858 (1995).
R. Colle, A. Curioni and O. Salvetti A non-local representation of the effective potential due to a molecular fragment, Theoret. Chim. Acta 86, 451-465 (1993).
I. Frank, S. Grimme, M. von Armin and S. D. Peyerimhoff The solvent shift in the n® p* excitation of CH2OnH2O: An MRD-CI investigation using effective potentials for the representation of the water molecules, Chem. Phys. 199, 145-153 (1995).
S. Katsuki The spectral representation technique for an active-electron-only molecular orbital calculation applied to the case of nonspherical frozen charge, J. Chem. Phys. 98, 496-501 (1993).
I. V. Abarenkov and I. M. Antonova Chemical bond modeling with the energy driven orbital localization, Int. J. Quantum Chem., Submitted.
R. Poteau, I. Ortega, F. Alary, A. Ramirez Solis, J.-C. Barthelat and J. P. Daudey Effective group potentials. I. Method, J. Phys. Chem. A 105, 198-205 (2001).
V. Kairys and J. H. Jensen QM/MM boundaries across covalent bonds: A frozen localized molecular orbital-based approach for the effective fragment potential method, J. Phys. Chem. A 104, 6656-6665 (2000).
M. S. Gordon, M. A. Freitag, P. Bandyopadhyay, J. H. Jensen, V. Kairys and W. J. Stevens The effective fragment potential method: A QM-based MM approach to modeling environmental effects in chemistry, J. Phys. Chem. A 105, 293-307 (2001).
B. O. Roos The complete active space self-consistent field method and its applications in electronic structure calculations, Adv. Chem. Phys. 69, 399-446 (1987).
S. Clifford, M. J. Bearpark and M. A. Robb A hybrid MC-SCF method: generalised valence bond (GVB) with complete active space SCF (CASSCF), Chem. Phys. Lett. 256, 320-326 (1996).
A. V. Soudackov, A. L. Tchougréeff and I. A. Misurkin Ground-state multiplicities and d-d excitations of transition-metal complexes by effective Hamiltonian method, Int. J. Quantum Chem. 58, 161-173 (1996).
A. M. Tokmachev, A. L. Tchougréeff and I. A. Misurkin Effective Hamiltonian approach to catalytic activity of transition metal complexes, Int. J. Quantum Chem. 84, 99-109 (2001).
S. Humbel, S. Sieber and K. Morokuma The IMOMO method: Integration of different levels of molecular orbital approximations for geometry optimization of large systems: Test for n-butane conformation and SN2 reaction: RCl+Cl-, J. Chem. Phys. 105, 1959-1967 (1996).
M. Svensson, S. Humbel and K. Morokuma Energetics using the single point IMOMO (integrated molecular orbital + molecular orbital) calculations: Choices of computational levels and model system, J. Chem. Phys. 105, 3654-3661 (1996).
R. D. J. Froese, S. Humbel, M. Svensson and K. Morokuma IMOMO(G2MS): A new high-level G2-like method for large molecules and its applications to Diels-Adler reactions, J. Phys. Chem. A 101, 227-233 (1997).
Th. Vreven and K. Morokuma On the application of the IMOMO (integrated molecular orbital + molecular orbital) method, J. Comp. Chem. 21, 1419-1432 (2000).
A. V. Nemukhin, B. L. Grigorenko, E. Ya. Skasyrskaya, I. A. Topol and S. K. Burt A new hybrid approach for modeling reactions in molecular clusters: Application for the hydrogen bonded systems, J. Chem. Phys. 112, 513-521 (2000).
B. L. Grigorenko, A. V. Nemukhin and V. A. Apkarian Many-body potentials and dynamics based on diatomics-in-molecules: Vibrational frequency shifts in ArnHF (n=112,62) clusters, J. Chem. Phys. 104, 5510-5516 (1996).
G. Náray-Szabó and P. Surján Bond orbital framework for rapid calculation of environmental effects on molecular potential energy surfaces, Chem. Phys. Lett. 96, 499-501 (1983).
G. G. Ferenczy, J. L. Rivail, P. R. Surján and G. Náray-Szabó NDDO fragment SCF approximation for large electronic systems, J. Comp. Chem. 13, 830-837 (1992).
Th. N. Truong and E. V. Stefanovich Development of a perturbative approach for Monte Carlo simulations using a hybrid ab initio QM/MM method, Chem. Phys. Lett. 256, 348-352 (1996).
I. Tuñón, I. T. C. Martins-Costa, C. Millot, M. F. Ruiz-López and J. L. Rivail A coupled density functional-molecular mechanics Monte Carlo simulation method, J. Comp. Chem. 17, 19-29 (1996).
Ch. D. Berweger, W. F. van Gunsteren and F. Müller-Plathe Finite element interpolation for combined classical/quantum mechanical molecular dynamics simulations, J. Comp. Chem. 18, 1484-1495 (1997).
O. A. Sharafeddin, K. Hinsen, T. Carrington Jr. and B. Roux Mixing quantum-classical molecular dynamics methods applied to intramolecular proton transfer in acetylacetone, J. Comp. Chem. 18, 1760-1772 (1997).
P. Grochowski, B. Lesyng, P. Bala and J. A. McCammon Density-functional based parametrization of a valence bond method and its applications in quantum-classical molecular dynamics simulations of enzymatic reactions, Int. J. Quantum Chem. 60, 1143-1164 (1996).
S. Chalmet, D. Rinaldi and M. F. Ruiz-López A QM/MM/continuum model for computations in solution: Comparison with QM/MM molecular dynamics simulations, Int. J. Quantum Chem. 84, 559-564 (2001).
A. L. Tchougréeff Group functions, Löwdin partition, and hybrid QC/MM methods for large molecular systems, Phys. Chem. Chem. Phys. 1, 1051-1060 (1999).
A. M. Tokmachev, A. L. Tchougréeff and I. A. Misurkin Effective electronic Hamiltonian for quantum subsystem in hybrid QM/MM methods as derived from APSLG description of electronic structure of classical part of molecular system, J. Mol. Struct. (Theochem) 506, 17-34 (2000).
M. V. Basilevsky and G. E. Chudinov Dynamics of charge transfer chemical reactions in a polar medium within the scope of the Born-Kirkwood-Onsager model, Chem. Phys. 157, 327-344 (1991).
J. Tomasi and M. Perisco Molecular interactions in solution: An overview of methods based on continuous distributions of the solvent, Chem. Rev. 94, 2027-2094 (1994).
T. Mineva, N. Russo and E. Sicilia Solvation effects on reaction profiles by the polarizable continuum model coupled with the Gaussian density functional method, J. Comp. Chem. 19, 290-299 (1998).
V. Barone, M. Cossi and J. Tomasi Geometry optimization of molecular structures in solution by the polarizable continuum model, J. Comp. Chem. 19, 404-417 (1998).
D. Wei and D. R. Salahub A combined density functional and molecular mechanics simulation of a quantum water molecule in aqueous solution, Chem. Phys. Lett. 224, 291-296 (1994).
M. Orozco, J. M. López, C. Colomines, C. Alhambra, M. A. Busquets and F. J. Luque Theoretical representation of solvation in biochemical systems: From discrete solute-solvent interactions to bulk solvation, Int. J. Quantum Chem. 60, 1179-1187 (1996).
H. Takahashi, T. Hori, H. Hashimoto and T. Nitta A hybrid QM/MM method employing real space grids for QM water in the TIP4P water solvents, J. Comp. Chem. 22, 1252-1261 (2001).
J. Gao Energy components of aqueous solution: Insight from hybrid QM/MM simulations using a polarizable solvent model, J. Comp. Chem. 18, 1061-1071 (1997).
P. Th. van Duijnen and A. H. de Vries Direct reaction field force field: A consistent way to connect and combine quantum-chemical and classical descriptions of molecules, Int. J. Quantum Chem. 60, 1111-1132 (1996).
G. L. Heard and B. F. Yates Hybrid supermolecule-polarizable continuum approach to solvation: Application to the mechanism of the Stevens rearrangement, J. Comp. Chem. 17, 1444-1452 (1996).
A. L. Tchougréeff SO(4) group and deductive molecular mechanics, J. Mol. Struct. (Theochem), submitted.
U. Eichler, Ch. M. Kölmel and J. Sauer Combining ab initio techniques with analytical potential functions for structure predictions of large systems: Method and application to crystalline silica polymorphs, J. Comp. Chem. 18, 463-477 (1997).
S. P. Greatbanks, I. H. Hillier and P. Sherwood Comparison of embedded cluster models to study zeolite catalysis: Proton transfer reactions in acidic chabazite, J. Comp. Chem. 18, 562-568 (1997).
R. Constanciel, In Localization and Delocalization in Quantum Chemistry, O. Chalvet, R. Daudel, S. Diner, J.P. Malrieu, eds. (D. Reidel Publishers, Dordrecht, 1975).
P.-O. Löwdin Studies in perturbation theory. IV. Solution of eigenvalue problem by projection operator formalism, J. Math. Phys. 3, 969-982 (1962).
R. McWeeny, Methods of Molecular Quantum Mechanics, Second edition (Academic Press, London 1992).
A. M. Tokmachev and A. L. Tchougréeff Potential energy surfaces in hybrid quantum mechanical/molecular mechanical methods, Int. J. Quantum Chem. 84, 39-47 (2001).
P. R. Surján, In Theoretical Models of Chemical Bonding, Part 2, Maksi\'c, Z.B., ed., (Springer, Heidelberg, 1989).
A. M. Tokmachev and A. L. Tchougréeff Generic molecular mechanics as based on local quantum description of molecular electronic structure, Int. J. Quantum Chem. 88, 403-413 (2002).
A. L. Tchougréeff and A. M. Tokmachev Deductive molecular mechanics of sp3 carbon atom, Int. J. Quantum Chem., submitted.
A. L. Tchougréeff and A. M. Tokmachev, in preparation.
R. Hoffmann An extended Hückel theory. I. Hydrocarbons, J. Chem. Phys. 39, 1397-1412 (1963).
M. Freindorf and J. Gao Optimization of the Lennard-Jones parameters for a combined ab initio quantum mechanical and molecular mechanical potential using the 3-21G basis set, J. Comp. Chem. 17, 386-395 (1996).
F. J. Luque, M. Bachs, C. Alemán and M. Orozco Extension of MST/SCRF method to organic solvents: Ab initio and semiempirical parametrization for neutral solutes in CCl4, J. Comp. Chem. 17, 806-820 (1996).
Y. Mo and J. Gao Ab initio QM/MM simulations with a molecular orbital-valence bond (MOVB) method: Application to an SN2 reaction in water, J. Comp. Chem. 21, 1458-1469 (2000).
T. Yoshida, N. Koga and K. Morokuma A combined ab initio MO-MM study on isotacticity control in propylene polymerization with silylene-bridged group 4 metallocenes. C2 symmetrical and asymmetrical catalysts, Organometallics 15, 766-777 (1996).
M. Sierka and J. Sauer Structure and reactivity of silic and and zeolite catalysts by a combined quantum mechanics-shell-model potential approach based on DFT, Faraday Discussions 106, 41-62 (1997).
R. V. Stanton, D. S. Hartsough and K. M. Merz Calculation of solvation free energies using a density functional/molecular dynamics coupled potential, J. Phys. Chem. 97, 11868-11870 (1993).
P. D. Lyne, M. Hodoscek and M. Karplus A hybrid QM-MM potential employing Hartree-Fock or density functional methods in the quantum region, J. Phys. Chem. A 103, 3462-3471 (1999).
X. P. Long, J. B. Nicholas, M. F. Guest and R. L. Ornstein A combined density functional theory/molecular mechanics formalism and its application to small water clusters, J. Mol. Struct. 412, 121-133 (1997).
S. A. French, A. A. Sokol, S. T. Bromley, C. R. A. Catlow, S. C. Rogers, F. King and P. Sherwood From CO2 to methanol by hybrid QM/MM embedding, Angew. Chem. 113, 4569-4572 (2001).
P. L. Cummins and J. E. Gready Coupled semiempirical molecular orbital and molecular mechanics model (QM/MM) for organic molecules in aqueous solution, J. Comp. Chem. 18, 1496-1512 (1997).
V. Théry, D. Rinaldi, J.-L. Rivail, B. Maigret and G. G. Ferenczy Quantum mechanical computations on very large molecular systems: The local self-consistent field method, J. Comp. Chem. 15, 269-282 (1994).
G. G. Ferenczy and J. G. Ángyán Intra- and intermolecular interactions in crystals of polar molecules. A study by the mixed quantum mechanics/molecular mechanical SCMP-NDDO method, J. Comp. Chem. 22, 1679-1690 (2001).
J. Å qvist and A. Warshel Simulation of enzyme reactions using valence bond force fields and other hybrid quantum/classical approaches, Chem. Rev. 93, 2523-2544 (1993).
D. Bakowies and W. Thiel Hybrid models for combined quantum mechanical and molecular mechanical approaches, J. Phys. Chem. 100, 10580-10594 (1996).
N. L. Allinger, Y. H. Yuh and J.-H. Lii Molecular mechanics. The MM3 force field for hydrocarbons. 1, J. Am. Chem. Soc. 111, 8551-8565 (1989).
A.L. Tchougréeff Separation of electronic variables in quantum chemical problems on potential energy surfaces of large molecular systems, Khim. Fiz. 16, 62-77 (1997) [Chem. Phys. Reports 16, 1035 (1997)].
T. Matsubara, S. Sieber and K. Morokuma A test of the new ''integrated MO + MM'' (IMOMM) method for the conformational energy of ethane and n-butane, Int. J. Quantum Chem. 60, 1101-1109 (1996).
B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan and M. Karplus CHARMM - A program for macromolecular energy, minimization, and dynamics calculations, J. Comp. Chem. 4, 187-217 (1983).
A. J. Turner, V. Moliner and I. H. Williams Transition-state structural refinement with GRACE and CHARMM: Flexible QM/MM modeling for lactate dehydrogenase, Phys. Chem. Chem. Phys. 1, 1323-1331 (1999).
W. D. Cornell, P. Ciepak, C. I. Bayly, I. R. Gould, K. Merz, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell and P. Kollman A 2nd generation of force field for the simulation of proteins, nucleic acids, and organic molecules, J. Am. Chem. Soc. 117, 5179-5197 (1995).
B. Waszkowycz, I. H. Hillier, N. Gensmantel and D. W. Payling Combined quantum mechanical-molecular mechanical study of catalysis by the enzyme phospholipase-A2 - an investigation of the potential energy surface for amide hydrolysis, J. Chem. Soc., Perkin Trans. 2, 2025-2032 (1991).
D. Yin and A. D. MacKerell Jr. Combined ab initio/empirical approach for optimization of Lennard-Jones parameters, J. Comp. Chem. 19, 334-348 (1998).
V. V. Vasilyev, A. A. Bliznyuk and A. A. Voityuk A combined quantum chemical/molecular mechanical study of hydrogen-bonded systems, Int. J. Quantum Chem. 44, 897-930 (1992).
M. A. Thompson, E. D. Glendening and D. Feller The nature of K+/crown ether interactions: A hybrid quantum mechanical-molecular mechanical study, J. Phys. Chem. 98, 10465-10476 (1994).
M. Eichinger, P. Tavan, J. Hutter and M. Parrinello A hybrid method for solutes in complex solvents: Density functional theory combined with empirical force fields, J. Chem. Phys. 110, 10452-10467 (1999).
M. A. Thompson and G. K. Schenter Excited states of the bacteriochlorophyll b dimer of rhodopseudomonas viridis: A QM/MM study of the photosynthetic reaction center that includes MM polarization, J. Phys. Chem. 99, 6374-6386 (1995).
I. B. Bersuker, M. K. Leong, J. E. Boggs and R. S. Pearlman A method of combined quantum mechanical (QM)/molecular mechanics (MM) treatment of large polyatomic systems with charge transfer between the QM and MM fragments, Int. J. Quantum Chem. 63, 1051-1063 (1997).
K. P. Eurenius, D. C. Chatfield, B. R. Brooks and M. Hodoscek Enzyme mechanisms with hybrid quantum and molecular mechanical potentials. I. Theoretical considerations, Int. J. Quantum Chem. 60, 1189-1200 (1996).
M. J. Harrison, N. A. Burton and I. H. Hillier Catalytic mechanism of the enzyme papain: Predictions with a hybrid quantum mechanical/molecular mechanical potential, J. Am. Chem. Soc. 119, 12285-12291 (1997).
V. V. Vasilyev J. Mol. Struct. (Theochem) 304, 129 (1994).
A. H. de Vries, P. Sherwood, S. J. Colins, A. M. Rigby, M. Rigutto and G. J. Kramer Zeolite structure and reactivity by combined quantum-chemical-classical calculations, J. Phys. Chem. B 103, 6133-6141 (1999).
N. Reuter, A. Dejaegere, B. Maigret and M. Karplus Frontier bonds in QM/MM methods: A comparison of different methods, J. Phys. Chem. A 104, 1720-1735 (2000).
Hyperchem, Inc. Hyperchem Users Manual; Computational Chemistry; Hypercube, Inc.: Ontario, Canada, 1994.
P. Cummins and J. Gready Combined quantum and molecular mechanics (QM/MM) study of the ionization state of 8-methylpterin substrate bound to dihydrofolate reductase, J. Phys. Chem. B 104, 4503-4510 (2000).
S. Ranganathan and J. Gready Hybrid quantum and molecular mechanical (QM/MM) studies on the pyruvate to L-lactate interconversion in L-lactate dehydrogenase, J. Phys. Chem. B 101, 5614-5618 (1997).
I. Antes and W. Thiel Adjusted connection atoms for combined quantum mechanical and molecular mechanical methods, J. Phys. Chem. A 103, 9290-9295 (1999).
Y. Zhang, T.S. Lee and W. Yang A pseudobond approach to combining quantum mechanical and molecular mechanical methods, J. Chem. Phys. 110, 46-54 (1999).
Y. Zhang, H. Liu and W. Yang Free energy calculation on enzyme reactions with an efficient iterative procedure to determine minimum energy paths on a combined ab initio QM/MM potential energy surface, J. Chem. Phys. 112, 3483-3492 (2000).
J. A. DiLabio, M. M. Hurley and P. A. Christiansen Simple one-electron quantum capping potentials for use in hybrid QM/MM studies of biological molecules, J. Chem. Phys. 116, 9578-9584 (2002).
D. Bakowies, PhD Thesis, Universität Zürich, 1994.
G. J. Kramer and R. A. van Santen Theoretical determination of proton affinity differences in zeolites, J. Am. Chem. Soc. 115, 2887-2897 (1993).
M. Svensson, S. Humbel, R. D. J. Froese, T. Matsubara, S. Sieber and K. Morokuma ONIOM: A multilayered integrated MO + MM method for geometry optimizations and single point energy predictions. A test for Diels-Alder reactions and Pt(P(t-Bu)3)2 + H2 oxidative addition, J. Phys. Chem. 100, 19357-19363 (1996).
J. R. Shoemaker, L. W. Burggraf and M. S. Gordon SIMOMM: An integrated molecular orbital/molecular mechanics optimization scheme for surfaces, J. Phys. Chem. A 103, 3245-3251 (1999).
G. Monard, M. Loos, V. Théry, K. Baka and J. L. Rivail Hybrid classical quantum force field for modeling very large molecules, Int. J. Quantum Chem. 58, 153-159 (1996).
S. Antonczak, G. Monard, M. F. Ruiz-López and J. L. Rivail Modeling of peptide hydrolysis by thermolysin. A semiempirical and QM/MM study, J. Am. Chem. Soc. 120, 8825-8833 (1998).
N. Ferré, X. Assfeld and J.-L. Rivail Specific force field parameters determination for the hybrid ab initio QM/MM LSCF method, J. Comp. Chem. 23, 610-624 (2002).
M. J. S. Dewar, E. G. Zoebisch, E. Healy and J. J. P. Stewart AM1: A new general purpose quantum mechanical molecular model, J. Am. Chem. Soc. 107, 3902-3909 (1985).
J. Gao, P. Amara, C. Alhambra and M. J. Field A generalized hybrid orbital (GHO) method for the treatment of boundary atoms in combined QM/MM calculations, J. Phys. Chem. A 102, 4714-4721 (1998).
P. Amara, M. J. Field, C. Alhambra and J. Gao The generalized hybrid orbital method for combined quantum mechanical/molecular mechanical calculations: formulation and tests of the analytic derivatives, Theoret. Chem. Acc. 104, 336-343 (2000).
G. Náray-Szabó Chemical fragmentation in quantum mechanical models, Computers & Chemistry 24, 287-294 (2000).
D. M. Philipp and R. A. Friesner Mixed ab initio QM/MM modeling using frozen orbitals and tests with alanine dipeptide and tetrapeptide, J. Comp. Chem. 20, 1468-1494 (1999).