Received 30 April, 2001; accepted in revised form 30 July, 2001

Abstract

hybrid QM/MM methods, junction between subsystems, APSLG wave function

Introduction

The hybrid QM/MM schemes for calculating large molecular systems acquire an increasing popularity []. In the literature there exist a large variety of hybrid QM/MM approaches differing by the QM and MM methods chosen and by specific forms of junction between classical and quantum subsystems. The common objects for applying the QM/MM schemes are large biological systems. Also other problems in the area of computational chemistry are related to the QM/MM setting: these are the problem of embedding in the cluster calculations on solids and surfaces; the problem of description of solute/solvent effects for reactions in condensed media; the pseudopotential descriptions of molecular and atomic electronic structure to mention several most important.

The purpose of the present work is to elucidate the precise form of the junction between the QM and MM regions in the hybrid QM/MM methods. It is more or less covered by the polarization effects when it takes place ''through space'' and there is no covalent bond between the subsystems. An opposite less transparent situation arises when chemical bonds exist between the QM and MM regions like in the case of reactions of large organic molecules.

It is now generally accepted that interregion boundary must not cut bonds because it leads to large fluctuations of electron density between subsystems causing numerous computational problems. Theoretical substantiation of such a choice of the interregion boundary has been given in Ref. TokmTchMis. The most popular way to perform the junction is based on using link atoms for saturation of free valence in the QM subsystem []. The link atom correction is, however, a poorly defined contribution to the energy. It leads, for example, to difference between potential and total energies, to collapse of link atoms with the neighbour QM atoms and to unphysical polarization of the boundary atoms. Also the question arises how to set the parameters for the link atoms.

Another group of methods uses localized orbitals [] coming from SCF procedures. These orbitals can follow the changes of molecular geometry on the basis of purely geometrical criteria. The main problem of these methods is the freezing of boundary one-electron states which does not allow to cover the effect of changes of the local environment upon them. Additional uncertainty is caused by poorly defined procedures of ''tail cutting'' of localized one-electron states. Also two other important features seem to be missed in the literature: (i) the renormalization of the MM part of the system due to variations of the electronic structure parameters (ESPs) on the boundary atom due to its participation in the QM subsystem, and (ii) the need that the adjustment of the QM residing hybrid orbitals (HOs) follows the variational principle for the energy when the geometry of the MM part changes. Here we present the formulae describing the linear response version of the theory resulting in the renormalization of the QM and MM parameters.

APSLG-MINDO/3 method and MM

To perform sequential derivation of hybrid QM/MM schemes some kind of derivation of the MM itself from specially designed (underlying) QM method is necessary. The method used is based on the trial wave function in the form of antisymmetrized product of strictly localized geminals (APSLG). Each geminal represents either a chemical bond or an electron lone pair. This method requires special choice of the Arai subspaces in which the geminals are constructed. These subspaces are formed by two (or one in the case of a lone pair) HOs t_m which in their turn are the result of an orthogonal transformation of the valence (s-, p-) AO basis for each heavy (non-hydrogen) atom:

tms+ = å i Î s,p  hmiais+.

(0)

The SO(4) matrices of the basis transformation h^A depend on six angles, which can be divided on two triples of angles: the first one ([(w)\vec] _b) defines the form of HOs (the relative weights of the s- and p-AOs in the HO) while the second one ([(w)\vec] _l) defines their orientation. Each HO contains scalar part s and vector part [(v)\vec]. The relation between the first order variations of the HOs and quasi- and pseudorotation angles [(w)\vec] _l and [(w)\vec] _b is given by:

d(1)s = -(d ® w   b  , ® v   ),d(1) ® v   = sd ® w   b  +d ® w   l  × ® v   .

(0)

Hybridization tetrahedron is defined by four vector parts of the HOs.

The singlet geminals are:

gm+ = umrma+rmb++vmlma+lmb++wm(rma+lmb++lma+rmb+),

(0)

where r_m and l_m are the HOs belonging to right and left atoms of the bond m respectively and u_m, v_m, and w_m are amplitudes of contributions to the m-th geminal (u_m²+v_m²+2w_m² = 1). The overall wave function is:

| F ñ = Õ m  gm+|0 ñ .

(0)

All the ESPs of the wave function (matrices h and geminal amplitudes) are then determined variationally.

Semiempirical implementation of this method using the MINDO/3 type of the Hamiltonian results in the following expression for the energy []:

E = 1 2 å A ¹ B  QAQBgAB+ å m  { å t Î {r,l}  [2UmtPmtt+(tmtm|tmtm)TmGmtt]+ +2gRmLm[Gmrl-2PmrrPmll]-4brmlmRmLmPmrl}+2 å k < m  å tt¢ Î {r,l}  dTkTm¢gtktm¢TkPkttPmt¢t¢,

(0)

where

gtktm¢Tk = 2(tktk|tm¢tm¢)Tk-(tktm¢|tm¢tk)Tk; QA = 2 å tm Î A  Pmtt-ZA; Pmtt¢ = á 0| gmtms+tms¢gm+| 0 ñ , Gmtt¢ = á 0| gmtmb+tma¢+tma¢tmbgm+| 0 ñ , Pmrr = um2+wm2, Pmll = vm2+wm2, Pmrl = Pmlr = (um+vm)wm, Gmrr = um2, Gmll = vm2, Gmrl = Gmlr = wm2

(0)

are the reduced Coulomb integrals, Mulliken charges, and matrix elements of one- and two-electron density matrices, respectively.

It was shown [] that the precision of results obtained by the APSLG-MINDO/3 method for alkanes, alyphatic amines and alcohols is somewhat better than that obtained by the standard SCF-MINDO/3 method. Moreover, the APSLG-MINDO/3 method possesses some special attractive features: (i) it operates with intuitive chemical concepts; (ii) it is a special case of O(N)-scaling methods and (iii) it has correct behaviour of the trial wave function for the whole range of interatomic separations.

The first feature and MM-like representation of the energy Eq. (5) allows one to use the APSLG-MINDO/3 method as a starting point to theoretically derive additive schemes including the MM. The derivation of deductive molecular mechanics can be done by fixing the ESPs corresponding to the geminals which perfect transferability was proven []. In the vicinity of the energy minimum the linear response relation between the variation of geometry parameters q and the reaction of the ESPs x to the latter holds:

x-x0 = -( ÑxÑxE) -1ÑxÑqE(q-q0).

(0)

This approach allowed to obtain the principal force fields for the sp³ carbon - bond stretching and bond bending (including the off-diagonal ones) - and the numerical values of their constants. They were in perfect agreement with values typically accepted in the current MM schemes. The precion of thus derived MM scheme is interpolated from those of the underlying QM description and linear response approximation to ESPs. The precision of the APSLG-MINDO/3 scheme is comparable with that of the standard (SCF-)MINDO/3 method. The precision of different approximations to the ESPs is also evaluated. For example, the explicit perturbative estimate [] of density matrix elements Eq. (6) even in the case of very polar water molecule leads to error only 0.014 kcal/mol in the heat of formation.

In the case of sp³ nitrogen the explicit form of the pyramidalyzation potential as a function of pyramidalization angle d is obtained:

2( Us-Up) tan2d+ æ ç è 3 2 C2+C3+2C5 ö ÷ ø tan2d-C3tan4d- -4Prl Ö3 é ë bssRmLm Ö   1-2tan2d   +Ö2bzsRmLmsindtand+Ö2bzsRmLmcosd ù û ,

(0)

where coefficients C_n are combinations of the Slater-Condon parameters:

C2 = 2F0(sp)+4G1(sp)-2F0(pp)-8F2(pp), C3 = F0(ss)-2F0(sp)-4G1(sp)+F0(pp)+4F2(pp), C5 = 2F0(sp)-G1(sp)-2F0(pp)+7F2(pp)

(0)

and U_s, U_p are matrix elements of electron attraction to the core.

Effect of the QM subsystem on the MM subsystem

The problem of describing the boundary atoms in the QM/MM junction can be solved in the present framework by noticing that the boundary atom keeps some of its HOs in the MM region and thus the geminal ESPs for them are fixed, whereas other orbitals are lended to the QM region.

Carbon boundary atom

Let us consider a boundary sp³ atom with one HO pointing to the QM region and others related to the MM one. The one- and two-electron density matrix elements for the QM residing HO are evaluated by the QM procedure of choice and thus differ from the fixed ESPs values accepted in the MM region. Setting m = 1 for the HO belonging to the QM subsystem and assuming the atom under consideration to be the ''right-end'' atom for the QM bond the perturbation to the one-center energy for the boundary atom due to participation of one of its HOs in the QM subsystem can be written as:

E1¢ = 2U1rdP1rr+(r1r1 | r1r1)R1dG1rr+dP1rr å k ¹ 1  gtkr1R1.

(0)

The derivatives of the latter expression with respect to the angles [(w)\vec] _l, [(w)\vec] _b yield additional quasi- and pseudotorques:

® K   ¢ 1  = Ñ[(w)\vec] lE1¢ = 0,  ® N   ¢ 1  = Ñ[(w)\vec] bE1¢.

(0)

In the QM subsystem the bond orders vary as well. The atoms in the QM part have nonvanishing off-diagonal elements of the one-electron density matrix between arbitrary one-electron states ascribed to the QM subsystem. The corresponding contribution to the energy reads:

Eres¢ = -2 å A  å m  dPr1mbr1mR1A,

(0)

where dP_l1m are the elements of the one-electron density matrix between the r₁-th HO assigned to the QM subsystem and the m-th AO on an atom A in the QM subsystem and b_r1m^R₁A are the corresponding resonance integrals. The resonance integrals take the following form (in the corresponding R₁A-diatomic coordinate frame (DCF) with the orts [(e)\vec]_R1A^x , [(e)\vec]_R1A^u , and [(e)\vec]_R1A = [(e)\vec]_R1A^z ):

br1sR1A = bssR1As1R1+bzsR1Av1zR1, br1zR1A = bszR1As1R1+bzzR1Av1zR1, br1xR1A = bppR1Av1xR1, br1uR1A = bppR1Av1uR1.

(0)

It leads to explicit expressions for the resonance contribution to the pseudo- and quasitorques affecting the hybridization tetrahedron at the boundary atom:

® K   ¢ res  = Ñ[(w)\vec]lEres¢,  ® N   ¢ res  = Ñ[(w)\vec]bEres¢.

(0)

The total pseudo- and quasitorques due to quantum behaviour of electrons in the QM subsystems are:

® N   ¢   = ® N   ¢ 1  + ® N   ¢ res  , ® K   ¢   = ® K   ¢ res  .

(0)

Pseudo- and quasirotations of the hybridization tetrahedron are then obtained in the linear response approximation by action of inverse of the Ñ_[(w)\vec] ²E matrix upon these torques. In symmetric sp³ case this matrix is diagonal with eigenvalues corresponding to pure variations of the angles [(w)\vec] _b and [(w)\vec] _l:

Db = 8Prl é ê ê ê ë sR æ è bssRLsL-bszRL   _____ Ö1-sL2   ö ø +   _____ Ö1-sR2 3 æ è bzsRLsL-bzzRL   _____ Ö1-sL2   ö ø ù ú ú ú û , Dl = 16 3 Prl   _____ Ö1-sR2   æ è bzsRLsL-bzzRL   _____ Ö1-sL2   ö ø .

(0)

These corrections result both in a new form and the orientation of the hybridization tetrahedron. The latter are inconsistent with the positions of the atoms bonded to the boundary atom from the MM side of the system. Multiplying the quasi- and pseudorotation angle corrections by the ( Ñ_[(w)\vec] Ñ_{[(j)\vec] m}E) ^f matrix results in a torque acting upon the L_m-th atom of the m-th bond incident to the frontier atom R₁ = R_m. Also the additional one- and two-electron densities on the boundary atom give additional forces acting upon its MM neighbors. They are obtained by multiplication of the mixed second derivatives matrix ( Ñ_[(w)\vec] Ñ_RLmRmE) ^f on the variations of the quasi- and pseudorotation angles. We derived the explicit expressions for the matrices of second derivatives of the energy, which are however too cumbersome to be given here. They allow to have explicit form for reaction of the MM part upon the variation of electron densities in the QM subsystem.

Effect of QM upon the nitrogen pyramidalization

It is possible that the boundary atom serving as the QM/MM junction is the sp³ nitrogen atom supplying its lone pair to the QM subsystem whereas three chemical bonds remain in the MM subsystem. In this case the same happens as above: the one- and two-electron density matrix elements change which produces the pseudo- and quasitorques acting on the nitrogen hybridization tetrahedron. The effect of these pseudo- and quasitorques is quite different since the second derivatives matrix Ñ_[(w)\vec] ²E must be now calculated for the sp³ nitrogen atom. We omit intermediate derivation because it is rather cumbersome and give here only the explicit expression for the deformation (pyramidalization) momentum:

dP1rr( 4Us-4Up+3C2+2C3+4C5-4C3tan2d) tand(1+tan2d),

(0)

which produces an additional deformation (flattening) of the nitrogen pyramid since for the lone pair the value of dP₁^rr can not be positive.

Effect of the MM subsystem on the QM subsystem

Any deformation in the MM system results in the variation of the pseudo- and quasirotation angles by Eq. (7). The shifts of the positions of the MM neighbours of the boundary atom result in quasi- and pseudotorques acting upon its hybridization tetrahedron. This produces variations of both one-center parameters correspoding to the QM residing HO and of the resonance parameters for the QM residing HO and all other orbitals in the QM region. The variation of the one-center matrix elements of the Hamiltonian corresponding to the QM HO is:

dU1r = 2(Us-Up)s1R1(d ® w   b  , ® v   R1 1  ), d( r1r1 | r1r1) R1 = ( 2C2s1R1+2C3(s1R1)3) (d ® w   b  , ® v   R1 1  ).

(0)

The renormalization of the QM resonanse intergrals to which the HO at hand is involved is:

dbr1sR1A = bssR1Ad(1)s1R1+bzsR1Ad(1)v1zR1, dbr1zR1A = bszR1Ad(1)s1R1+bzzR1Ad(1)v1zR1, dbr1xR1A = bppR1Ad(1)v1xR1,dbr1uR1A = bppR1Adv1uR1.

(0)

If the variations Eqs. (18), (19) are taken into account in the calculations on the QM part the effect of the geometry variations in the MM subsystem on the parameters of the electronic effective Hamiltonian for the QM part is taken into account without further reparameterization.

Acknowledgements

This work has been performed with financial support of RAS through the grant # 6-120 dispatched by its Young Researchers' Commission. Financial support for AMT on the part of the Haldor Topsoe A/S is acknowledged as well.

References

[]

CECAM-NSF Meeting on QC/MM methods, International Journal of Quantum Chemistry 60 (1996) No 6, Special Issue.

[]

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, Journal of Molecular Structure (Theochem) 506, 17-34(2000).

[]

D. Bakowies and W. Thiel, Hybrid models for combined quantum mechanical and molecular mechanical approaches, Journal of Physical Chemistry 100, 10580-10594 (1996).

[]

D.M. Philipp and R.A. Friesner, Mixed ab initio QM/MM modeling using frozen orbitals and tests with alanine dipeptide and tetrapeptide, Journal of Computational Chemistry 20, 1468-1494 (1999).

[]

A.M. Tokmachev and A.L. Tchougreeff, Semiempirical implementation of strictly localized geminals approximation for analysis of molecular electronic structure, Journal of Computational Chemistry 22, 752-764 (2001).

[]

A.L. Tchougréeff and A.M. Tokmachev, Submitted.

[]

A.M. Tokmachev and A.L. Tchougréeff, Generic molecular mechanics as based on local quantum description of molecular electronic structure, International Journal of Quantum Chemistry 88, 403-413 (2002).