A.L. Tchougréeff
Karpov Institute of Physical Chemistry
Vorontsovo pole 10, Moscow, 105064, RUSSIA
and
Center for Computational Chemistry at the M.V. Keldysh
Institute for Applied Mathematics of RAS

Deductive Molecular Mechanics as Applied to Develop
QM/MM Picture of Dative and Coordination Bonds

Abstract

1  Introduction

Chemical bonds are portrayed in textbooks as being nonpolar covalent, polar covalent, ionic, dative, donor-acceptor, coordination, and so on without giving clear reasons for ascribing whatever specific bond to that or another class []. A correct classification may, however, be important not only from the pedagogical point of view: from molecular mechanics (MM) experience we know that constructing a mechanistic description is not equally easily possible for different classes of compounds containing bonds of different types. For purely ''organic'' molecules with well defined bonds numerous empirical parameterizations have been developed quite successfully [,,]. Corresponding efforts when applied to metal containing compounds and hydrogen bonds until now did not give completely satisfactory result [,,]. Of course, a good number of works appeared which parameterize potential energy surfaces (PES) of some well defined classes of metal containing compounds, but some questions important also from the practical point of view remain unclear.

The problem addressed here concerns uniform description of the metal complexes with variable number of the ligands. Such a description becomes in demand in the molecular dynamics simulations (see e.g. MillotDehez) of metal ion behavior in solutions of complexing solvents containing also chelating ligands (e.g. crown-ethers, cyclic polyammines etc.). In such a system one may expect formation of numerous complexes with different number of solvent ligands and/or different mode (chelating number) of complexation which must be treated on equal footing in order to estimate their energies on a uniform scale. Clearly, the usual MM harmonic approximation for the metal-donor bond stretch energy accepted say in Ref. [] cannot account for such effects. Meanwhile, a simplistic replacement of a harmonic potential by another one with better asymptotic behavior (like the Morse potential) does not solve the problem since too many other factors seem to be important. These are the effective charge variation and nontrivial dependence of the interligand and metal-ligand interaction energy on metal ligand separations, to mention only a few.

The charge variability appears due to two types of interactions almost equally important in the case of metal ions binding by donor ligands: one is due to polarization of the ligands by the point metal ion and by the charges residing in other ligands as well; another is due to electron transfers from the donor atoms to the ions' empty shells (Lewis acid-base interactions). The importance of the former mechanism has been recently stressed in Ref. []. The same concept has been recently used while developing the COSMOS MM force field [] employed later for analysis of behavior of the Zn²⁺ complexes with nitrogen containing ligands Sternbergetal. The authors [] mentioned that the charge redistribution due to electron transfers is not important. This may be true for that class of objects the cited authors actually consider: the crown ether complexes of the Cs⁺ or Mg²⁺ ions, where the results of quantum chemical analysis reveal a negligibly weak transfer of electronic density from the oxygen donor atoms to the metal ions (though the calculated extent of this transfer is known to be ''method dependent''). Such a picture is not generally valid for all metals since some of them are much stronger Lewis acids than heavy alkali cations. For example, our older calculations on the charge distribution in the transition metal complexes revealed a general trend that the formally divalent transition metal cations bear an effective charge of about one unit charge, whereas for the trivalent cations the effective charge is less than two unit charges []. Even in less pronounced case of the Mg²⁺ ions coordinated through oxygens in xylose isomerase the effective charge obtained on Mg within the PM3 semiempirical calculation [] is close to unity. Similar picture has been reported in Ref. [] for Zn²⁺ complexes with imidazole. Remarkable role of the charge redistribution in close vicinity of the Ln³⁺ cations which does not reduce only to the polarization of the surrounding ligands has been reported by Malta et al. []. Thus the overall picture appears to be too confusing to hope that it can be disentangled by a combination of may be locally successful ad hoc recipes.

Previously we developed an approach which can be used to put the process of developing mechanistic descriptions of PES (i.e. of developing MM force fields) of different classes of compounds on a rational basis. It is the deductive molecular mechanics [,,] (DMM) which allows to develop a form of the MM force fields analyzing the form of the electronic wave function chosen in a form relevant to the physical picture of the considered class of molecules. In the present paper we apply the previously developed DMM approach to analytical derivation of the QM based form for the force fields involving the nontransition metal atoms. The DMM methodology which is a general framework for this description occupies a border position since it is designed to a bridge the gap between the QM and MM descriptions.

2  DMM Framework

Deductive molecular mechanics (DMM) is an assembly of approximations DMM,SO4DMM,SP3CDMM performed in the general QM/MM context which when applied lead to a mechanistic description of PES for molecular systems. The mentioned QM method [,,] involves the usage of the trial wave function in the from of the antisymmetrized product of strictly local geminals (APSLG, see [] and references therein) with semiempirical Hamiltonians. The DMM methodology consists of constructing direct estimates for the electronic structure parameters (ESP) which are defined in the APSLG framework in terms of the parameters of the effective Hamiltonian. The ESPs thus estimated are then inserted in the APSLG based expression for the energy of a molecule at hand. This gives the required mechanistic description of the PES.

Below we briefly describe the QM-APSLG method for molecular electronic structure underlying the MM description of molecular PES and describe the physical prerequisites and procedures previously used [] to construct the mechanistic description of the molecular PES in organic molecules. Next we analyze the physical picture corresponding to the ''dative'' or ''ionic'' bonds between metal ions (modeled by their acceptor orbitals) and molecules containing nitrogen and oxygen donor atoms and construct an APSLG-based QM/MM model for description for the corresponding contributions to the PES.

2.1  APSLG based QM method

In the semiempirical APSLG approximation the electronic wave function of the molecule has the form (see Ref. [,,]) of the antisymmetrized product of the two-electron functions (geminals):

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

(1)

The form of the geminals used in the second quantization notation is:

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

(2)

for (two-center) chemical bonds and:

gm+ = rma+rmb+

(3)

for lone pairs. The amplitudes of the configurations (u_m, v_m, and w_m) in Eq. (2) are determined with use of the variational principle. The normalization condition is imposed on the amplitudes:

um2+vm2+2wm2 = 1.

(4)

The operators t_ms⁺ (where t is either r or l) create electrons with spin projection s on the hybrid orbitals (HOs) | t_m ñ ascribed by the construction procedure for the function Eq. (1) to the m-th geminal.

The amplitudes (u_m,v_m,w_m) are determined by sequential diagonalizations of 3×3 matrices of effective Hamiltonians formed separately for each two-center bond (see below). Using the Hamiltonian of the MINDO/3 form MINDO/3 with the APSLG trial wave function Eqs. (1), (2) in Ref. [] resulted in a semiempirical APSLG-MINDO/3 QM method. The total energy then has a form closest to the standard MM energy expression []:

Etotal = å A  EA+ å A < B  EAB,  where EA = å m Î A  { å t Î {r,l} ÇA  [2UmtPmtt+(tmtm | tmtm)TkGmtt]+ + å k < m  å tt¢ Î { r,l} ÇA  2gtktm¢TkPkttPmt¢t¢}. ERmLmbond = 2gRmLm[Gmrl-2PmrrPmll]-4brmlmRmLmPmrl EABnonbond = 1 2 ( ZAZBDAB+QAQBgAB)

(5)

in that sense that the overall energy expression is a sum of local (atom and bond) contributions. In the above expression R_mL_m refer to the right- and left-end atoms of the m-th bond. Each contribution is written through the intrabond matrix elements of one- and two-electron densities:

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,

(6)

and molecular integrals transformed to the basis of strictly local HOs residing at each ''heavy'' (nonhydrogen) atom:

Umt = sm2(Us-Up)+Up ( tmtm | tmtm) = C1+C2sm2+C3sm4, gtktm¢ = C4+C5[sm2+sk2]+C3sm2sk2 brmlmRmLm = bssRmLmsmRmsmLm+bszRmLmsmRm( ® v   Lm m  , ® e   RmLm  )+ + bzsRmLm( ® v   Rm m  , ® e   RmLm  )smLm+bppRmLm( ® v   Rm m  , ® v   Lm m  )+ + ( bzzRmLm-bppRmLm) ( ® v   Rm m  , ® e   RmLm  )( ® v   Lm m  , ® e   RmLm  )

(7)

where C_n's are the combinations of the Slater-Condon parameters of the intraatomic Coulomb electron-electron interactions introduced in Ref. DMM, U_s and U_p are the parameters of the electron attraction to the core on s- and p-orbitals of the heavy atom; b_mn^R_mL_m are the resonance integrals in the diatomic coordinate frame where the z-axis coincides with the interatomic vector [(e)\vec]_RmLm PopleBeveridge (hereinafter we use an arrow to indicate a three-dimensional vector). For explanation of other elements of the above expression vide infra.

The density matrix elements Eq. (6) comprise the subset of the QM ESP's related to the geminal amplitudes. The diagonal geminal-related ESP's define also the effective atomic charges:

QA = 2 å tm Î A  Pmtt-ZA,

(8)

where Z_A's stand for the core charges.

The molecular integrals in Eq. (7) for the energy expression Eq. (5) are written in the basis of one-electron HOs given by the formula (Ref. []):

| tm ñ = å i Î AO  hmit|i ñ ;t = r,l;i = s,x,y,z

(9)

The h matrices (of their own for each ''heavy'' - nonhydrogen - atom) are the SO(4) matrices due to obvious orthonormality conditions. They result in another set of ESPs related to the APSLG picture of molecular electronic structure which is the set of sextuples of atom-specific Jacobi angles three of which ([(w)\vec] _b = (w_sx,w_sy,w_sz)) - for each heavy atom - define the overall shape of the system of the HOs residing at it and other three ([(w)\vec] _l = (w_yz,w_xz,w_xy)) define the orientation of the whole system of the HOs at a respective heavy atom. Following [,] the strictly local HOs | t_m ñ allow also for the quaternion representation. Indeed, a quaternion [] may be characterized as an entity comprising a scalar and a 3-vector parts: q = (s,[(v)\vec]). The coefficient s of the s-orbital does not change under the spatial rotation of the molecule, whereas the coefficients at the p-functions transform as if they were the components of the 3-dimensional vector [(v)\vec] leaving incidentally intact the total weight of the p-component of the HO which equals to its squared norm. Thus each of the HOs located at a heavy atom can be presented as a normalized quaternion:

qm = (sm, ® v   m  ), sm2+ ê ê ® v   m  ê ê 2   = 1.

(10)

Four quaternions residing at a given atom are orthogonal in the usual sense:

smsm¢+( ® v   m  , ® v   m¢  ) = 0;m ¹ m¢

For an arbitrary HO in the quaternion representation the first order correction to its components occurring when the Jacobi angles get small variations d[(w)\vec] _l and d[(w)\vec] _b acquire a particularly simple form [,,]:

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

(11)

The Jacobi angles at each atom together with the elements of the density matrices Eq. (6) comprise the set of the ESPs defined in the semiempirical APSLG framework. Both the amplitudes u_m,v_m,w_m for each bond (or an equivalent set of the geminal related density ESPs Eq.( 6)) and the Jacobi angles ([(w)\vec] _b,[(w)\vec] _l) for each heavy atom are determined variationally in Refs. [,,] by minimizing the energy Eq. (5).

The Jacobi angles themselves are not, however, visual enough for representing the systems of HOs at the heavy atoms of the molecule. For this reason in Refs. [,,] it was proposed to consider ''hybridization tetrahedra'' formed at each heavy atom by four corresponding vector parts [(v)\vec]_m,m = 1¸4 of HO quaternions residing there. The orientation and the shape of these tetrahedra are uniquely related with the corresponding properties of the represented system of the HOs. The interhybrid angles q_mm^¢ (key invariants of the shape of the HOs' system) are given by:

cosqmm¢ = ( ® v   m  , ® v   m¢  ) ê ê ® v   m  ê ê ê ê ® v   m¢  ê ê = - smsm¢   _____ Ö1-sm2 Ö 1-sm¢2

(12)

where the latter equality appears due to orthonormalization condition Eq. (10) imposed on the quaternions representing HOs. We also notice here that due to normalization condition:

4 å m = 1  sm2 = 1

(13)

for the weights of the s-AO in all HOs the shape of the hybridization tetrahedron is uniquely defined by setting any three of the four coefficients s_m. This together with Eq. (11) immediately results in an explicit linear dependency of the vector parts [(v)\vec]_m,m = 1¸4:

4 å m = 1  sm ® v   m  = 0

(14)

The DMM is developed in [] terms of interactions between the hybridization tetrahedra which are considered as true building blocks of the mechanistic description for the molecular PES instead of commonly accepted ''balls and springs''. The tetrahedra interact according to formulae Eqs. (5)-(7) which depend both on molecular geometry and on shapes and mutual orientation of hybridization tetrahedra. The precise form of the DMM force fields depends on the approximations employed for the density matrix ESPs, which are accepted on a basis of analysis of the physical situation. For example, the DMM of ''organic'' molecules is derived in Ref. [] for the bond geminals close to the SCF closed shell two-electron function residing in the Arai subspace spanned by one-electron functions | t_m ñ . In this case the density matrix ESPs can be naturally represented in the form:

Gmtt¢ = 1 4 +dGmtt¢, Pmtt¢ = 1 2 +dPmtt¢.

(15)

A series of approximations is then possible which either set the elements of the density matrices equal to their invariant values (dP_m^tt¢ = 0,dG_m^tt¢ = 0 for covalent bonds; dP_m^tt¢ = [1/2],dG_m^tt¢ = [3/4] for lone pairs) or treat the small corrections dP_m^tt¢,dG_m^tt¢ perturbatively. With these assumptions the optimization of the total energy becomes a mechanistic procedure in a course of which the relative positions, orientations, and shapes of the hybridization tetrahedra are adjusted to reach the energy minimum.

3  APSLG approximation and dative bond

The derivation of ''organic'' DMM was performed in Ref. [,] assuming that the single bond is close to the SCF two-electron function. This resulted in the approximation Eq. (15) for the density matrix elements. Validity of this picture can be illustrated by the date given in Table . However, the infinite bond-length asymptotic wave function of two electrons forming a single bond between two atoms is the Heitler-London singlet with two electrons with equal probability residing one by one on either end of this bond. This corresponds to the limit to which the covalent nonpolar bond flows when the interatomic distance increases (the homolytic cleavage of a s-bond). Dative bonds by contrast flow to the ionic limit Eq. (3) in the case of infinite bond elongation. This prevents from using the approximation Eq. (15) for the geminal ESPs since they correspond to the wave function with different asymptotic behavior. Our immediate purpose is to obtain estimates for the geminal ESPs in the ionic limit.

3.1  APSLG analysis of ionic limit

3.1.1  Effective bond Hamiltonian in the ionic limit

The geminal wave functions Eq. (2) in the APSLG approximation are by definition obtained by diagonalizing the effective Hamiltonian for the m-th bond:

æ ç ç ç è Hmeff ö ÷ ÷ ÷ ø æ ç ç ç è um zm vm ö ÷ ÷ ÷ ø = em æ ç ç ç è um zm vm ö ÷ ÷ ÷ ø ,

(16)

where z_m = Ö2w_m and

Hmeff = Hm0+Hm¢ Hm0 = æ ç ç ç è a 0 0 0 b 0 0 0 c ö ÷ ÷ ÷ ø ;Hm¢ = d æ ç ç ç è 0 1 0 1 0 1 0 1 0 ö ÷ ÷ ÷ ø

(17)

with

a = emR b = 1 2 (emR+emL)-Dgm c = emL d = -Ö2brmlmRmLm

where

emR = 2Umr+(rmrm|rmrm)Rm-4gRmLmPmll+ +2 å B ¹ Rm  gRmBQB+2 å \Sb tm1 Î Rm m1 ¹ m\endSb grmtm1RmPm1tt; emL = 2Uml+(lmlm|lmlm)Lm-4gRmLmPmrr+ +2 å B ¹ Lm  gLmBQB+2 å \Sb tm1 Î Lm m1 ¹ m\endSb glmtm1LmPm1tt Dgm = 1 2 å t = r,l  ( tmtm | tmtm)Tm-gRmLm > 0

(18)

The asymptotic ground state of H_m^eff is controlled by relative position of the H_m⁰ nonzero diagonal elements on the energy scale. If b is the lowest diagonal matrix element of H_m⁰ the asymptotic ground state is the Heitler-London one (''organic'' case). An alternative situation takes place if the lowest diagonal matrix element of H_m⁰ equals (for the sake of definiteness) to a. This is possible if:

0 < 1 2 (emL-emR)-Dgm

(19)

In the above inequality all components depend on interatomic distance. Thus the character of the bond may (as it is well known, see Ref. ShaikShurki for the recent analysis of this situation in similar terms) change in the course of its formation or cleavage. The critical distance r_RmLm^c defined by the relation:

Dgm(rRmLmc) = 1 2 (emL(rRmLmc)-emR(rRmLmc))

(20)

separates the regions with different physical regime: at the distances shorter than the critical one independently turning the resonance interaction off leads to the Heitler-London wave function; at the longer distances turning the resonance off leads to the ionic wave function. The latter is obviously the lone pair wave function Eq. (3). This consideration allows to precise the definition of dative bonds. These are those which flow to the single ionic configuration in the zero resonance limit without direct relation to the interatomic separation. Namely this function must be used as a zero approximation for constructing the estimates for the density ESPs for the dative bond geminals.

3.1.2  Density ESPs in the ionic limit

The previous application [,,] of the DMM approach was based on the fact that for the geminals having the Heitler-London wave function as their asymptotic limit the related ESPs can be taken as transferable quantities. For example, it was shown that the spin bond orders P_m^rl in Eq. (6) can be set equal to 1/2 which is equivalent to the well known picture of the localized single bonds in organic molecules. For a considerable range of the bond-lengths around the equilibrium ones the above value of the spin bond-order remains invariant up to second order with respect to a small parameter and the interatomic distance dependence of the bond energy is dominated by that of the resonance (one-electron hopping) integrals between the left-end and right-end HOs ascribed to the bond under consideration. The deviations of the geminal-related ESPs from their invariant (transferable) values could be treated as a relatively small corrections. In the case of the ''dative'' bond the ionic configuration corresponding to accommodation of both bond electrons on (for the sake of definiteness) the right-end atom is the asymptotic wave function for the separated metal ion and the ligand molecule. This differs significantly from the ''organic'' situation so that the results of Refs. [,,] cannot be employed.

As we mentioned previously the density ESPs must be evaluated through the effective bond Hamiltonian Eq. (17). If this is done directly (without intermediate calculation of quantum amplitudes) corresponding quantities are the measurables in the quantum mechanics terms. In terms of the geminal amplitudes the ESPs are given by Eq. (6). To get the required direct estimates we turn to the projection operator technique. According to Refs. [,] whatever projection operator P can be obtained from the not too distant unperturbed operator P₀ by the following formula:

P = (P0+V)(1+V+V)-1(P0+V+), dim\limfuncImP = dim\limfuncImP0

(21)

where the matrix V and its hermitean conjugate V⁺ satisfy the conditions:

P0V = 0;VP0 = V;(1-P0)VP0 = V; P0V+ = V+;P0V+(1-P0) = V+;V+P0 = 0

(22)

which is the formal manifestation for the ''off-diagonal'' block structure of the matrices V and V⁺. The matrix V is obviously an dim\limfuncImP₀×dim\limfuncIm(1-P₀) matrix. The projection operator P₀ reappears if V = V⁺ = 0 thus ensuring the correct asymptotic behavior of the proposed system of ESPs. The order of the square matrix V⁺V equals to dim\limfuncImP₀ so that at modest dimensionalities of the subspace in question it may be easy to get an inverse matrix of Eq. (21). Alternatively, making use of the recursion relation:

(1+V+V)-1 = 1- V+V 1+V+V

(23)

we arrive to the following form of the projection operator Eq. ( 21):

P = P0+(V+V+)+VV+- V+V 1+V+V - VV+V 1+V+V - V+VV+ 1+V+V - VV+VV+ 1+V+V

(24)

Iterating this move one can easily obtain the expansion of P in a power series in V and V⁺.

Now we address the projection operator upon the ground state of a geminal. It is known in terms of the geminal amplitudes:

P = æ ç ç ç è u2 uz uv uz z2 zv uv zv v2 ö ÷ ÷ ÷ ø

(25)

The zero approximation for the projection operator P₀ to the ionic limit ground state corresponds to u = 1;z = v = 0. To ensure the correct (ionic) limit of a general one dimensional projection operator in a three dimensional space we apply the above prescription with the notion that dim\limfuncImP₀ = 1 and dim\limfuncIm(1-P₀) = 2. In this case the matrix block V consists of one row and two columns so that there are two independent parameters. Clearly the one-times-one matrix can be easily inverted and this results in the following form of the ground state projection operator:

P = 1 1+x2+y2 æ ç ç ç è 1 x y x x2 xy y xy y2 ö ÷ ÷ ÷ ø ,

(26)

where x and y are arbitrary real parameters. The projection property P² = P is checked immediately as well as the fact that \limfuncSpP = 1. P₀ obviously reappears when x = y = 0.

The Schrödinger equation in terms of the projection operator reads PupyshevStepanov:

HP = PH

(27)

which can be recast to a system of nonlinear equations for x and y (details are given in the Appendix):

x = d (a-b) ( 1+y-x2) y = d (a-c) x( 1-y) xy = d (b-c) ( x2-y-y2) ü ï ï ï ï ï ï ï ý ï ï ï ï ï ï ï þ

(28)

The last equation for the product xy must be inserted into one for y and the system becomes one for x and y. Solving this system will be equivalent to solving the original 3×3 eigenvalue problem for the effective bond Hamiltonian. In a perturbative manner we get for the first order approximation:

x = d (a-b) ;y = 0

(29)

In the weak bonding limit the ESP parameter y is one of higher order in the small parameter d than x. Identifying the matrix elements of the projection operators Eqs. (25), ( 26) establishes the relation between two parameterizations of the three-dimensional projection operators and produces explicit forms for the ESPs Eq. (6) in terms of x:

Prr » 1- x2 2 ;Pll » x2 2 ;Prl » x Ö2 ; Grr » 1-x2;Grl » x2 2 ; dPrr » 1 2 - x2 2 ;dGrr » 3 4 -x2.

(30)

3.1.3  Bond energy in the ionic limit

In the previous Subsection we obtained the picture of the geminal-related ESPs for the ''dative'' or ''donor-acceptor'' bond. This picture is different from the picture of ''polar covalent'' bond which may be characterized as an asymmetric one but nevertheless retaining the Heitler-London limit when b® 0. With use of the above results concerning the geminal related ESPs in the ionic limit we notice that the key estimate of the ''organic'' DMM becomes invalid: the geminal related ESPs even approximately are not the transferable quantities. By contrast, the values of all density matrix elements become at least linearly dependent on the resonance integral between the donor HO and the acceptor orbital. Together with the exponential distance dependence of the resonance integral this produces strong anharmonic potential characteristic for large interatomic separations.

The energy itself of the geminal in the ionic limit (ionic bond) can be easily estimated:

EDAbond » - 8bDA2 (emA-emD)-2Dgm æ ç è 1+ gDA (emA-emD)-2Dgm ö ÷ ø

(31)

This has a natural form of the perturbative estimate of the resonance energy known in the theory.

3.1.4  Bonding contribution to the QM/MM description of dative bond

As it is seen from the perturbative estimates the density ESPs are not transferable and are rather sophisticated functions of other ESPs which define the shape and orientation of the hybridization tetrahedron on the donor atom. At this stage it is possible alternatively to trying to develop a mechanistic description to retain the variable x (and, may be, y) to reproduce details of the electronic structure. This is equivalent to using a QM description specifically for the dative bond. Incidentally this picture is an adequate form suitable for constructing a QM/MM junction. In mechanistic terms their meaning can be interpreted as a strongly anharmonic spring. This form must be used instead either of the energy of a covalent bond with constant (transferable) spin bond-orders characteristic for the usual covalent bonds or of the perturbative bond energy estimate. Though they are still rather sophisticated functions of the internuclear separation and of the shape and orientation of the hybridization tetrahedron residing on the donor atom they can be easily calculated since require only the elementary functions for their evaluation [] and have correct asymptotic behavior at infinite internuclear separation thus being potentially useful in the molecular dynamics modeling of complexes with variable number of coordination bonds between the metal atom and ligand donor atoms. This picture of the dative bonding is nevertheless a QM one though very much reduced.

We also notice that according to the perturbative estimates for the solutions of the system Eq. (28) the equilibrium value of the y variable is by one order of magnitude in b_DA higher than x. Since only the combinations xy and y² enter in the expression of the projection operator and thus in that for the energy we may set y = 0 without introducing too large error in energy and by this to further simplify the projection operator Eq. (26). With use of this approximation the bond energy becomes:

EDAbond » - 2Ö2x 1+x2 (bDA-gDAx)

(32)

where x is an equal in rights variable of a QM/MM description. It must be determined from the general energy optimization procedure.

3.1.5  Qualitative picture of the DMM force fields at donor atoms.

In the previous Section we presented a derivation of the DMM force fields related to the donor atom with lone pair interacting with metal ions. This derivation resulted in rather complicated formulae where all the terms depend on details of hybridization of the donor atom through the resonance integrals b_DA. One may foresee two situations: (i) the shape of the hybridization tetrahedron is to major extent defined by the shape of the molecule itself; (ii) the shape of the hybridization tetrahedron is rather flexible and its variation in the course of formation of the dative bond is significant.

Indeed, when the dative bonds are formed by the sp³-nitrogen atoms the shape of their hybridization tetrahedra do not change. This can be understood on the basis of the linear dependency condition Eq. ( 14). Three covalent bonds fix the norm and the directions of three vector parts [(v)\vec]_m of the four. Then the fourth is fixed by the cited linear condition. In this case the dependence of the dative bond energy on geometry parameters can be obtained, say, from analysis of that of the perturbative estimate Eq. (31). In the latter expression the molecular geometry dependence is dominated by the square of the resonance integral. Assuming as in the previous Section that the acceptor is represented only by its empty s-orbital we get as in the case of a hydrogen atom a significant simplification for the resonance integral of the dative bond:

bDA = bssDAs1+bzsDA( ® v   1  , ® e   1  )

(33)

The relative orientation of the acceptor atom and the donor atom hybridization tetrahedron enters through the scalar product of the vector part of the lone pair involved in the formation of the dative bond ([(v)\vec]₁ ) and the dative bond vector ([(e)\vec]₁) which is a unit vector directed along the line connecting the donor and acceptor atoms. The maximum of the above expression is obviously reached when [(v)\vec]₁||[(e)\vec]₁ and whatever escape from this line producing a ''lone-pair misdirection'' contributes to the corresponding MM force field for the A-D-X angles (where X stands for whatever atom covalently bound to the donor one). It must be noticed, however, that the corresponding contribution is by no means the leading one in terms of energetics since one may expect much greater importance of purely electrostatic forces in this situation. Indeed, the dative bonding contribution to the ''misdirect potential'' is scaled down by a small value of the ESP x Eqs. (26), ( 30) - spin bond order of the dative bond. At the same time the general context of the dative bond formation allows to expect that a local dipole moment [(m)\vec] _D resides on the donor molecule and the ion-dipole energy term:

QAe( ® m   D  , ® e   1  ) rDA2

where Q_A is the effective charge of the acceptor atom (close to its formal ionic charge) and r_DA is the donor-acceptor interatomic separation is larger than the dative bonding contribution. For the symmetry reasons (we refer here to a model C_3v symmetric ''ammonia'' molecule) one may expect that [(v)\vec]₁||[(m)\vec] _D and thus the resulting force field which tends to align [(m)\vec] _D and [(e)\vec]₁ and also acts to prevent the lone-pair misdirection, though does not have too much relation to the lone pair direction at all.

3.2  QM/MM model for metal complexes

In the previous Sections we applied the DMM methodology to derive formulae for the energy (force fields) for donor atoms interacting with complexing agents. A simplified representation for the acceptor with use of a single s-orbital was used. In this Section we consider the metal-ligand interactions from a slightly different point of view. The metal ion in a complex acquires some density not from one but from many lone pairs of the ligating donor atoms. Constructing of a mechanistic or at least an economic QM description for such a case would possibly help to rationale terms of interligand interaction force fields which are sometimes included in the standard MM picture to assure proper description of the metal complexes.

Metal ions (both with and without open d-shell) stay aside from the general MM picture based on the concept of localized transferable two-center bonds. The physical reasons are the specific properties of metal-ligand bonds such as lack of saturability, directionality (see[]), and transferability which are sufficient components of the standard MM picture. To be more precise the metal-ligand bonds lack directionality at metal center, though the directionality at donor atoms exists and the corresponding effects are termed as ''misdirect'' of lone pairs. These well known properties represent to our opinion the reasons why despite numerous attempts present in the literature (for review see e.g. Refs. Comba,Hay) the PES's of metal complexes are not easily covered by the MM-like schemes. In this context an attempt to deduce a QM/MM description for (nontransition) metal atoms using the general methodology (selecting the relevant form of the trial electronic wave function first, and extracting convenient ESPs to be used in the model energy expression) looks very much desirable. A necessity of namely QM/MM approach in this context is of course the flexibility of the corresponding ESPs which have to be recalculated for different molecular geometries in an economical way.

3.2.1  Physical concepts for selecting the form of wave function

The physical concepts of the metal-ligand bonding mentioned above are largely negative statements. The metal-ligand bonds are nontransferable, nonlocalized, nondirectional at the metal site and directional on the ligand side. Thus the trial wave function for the metal complexes is generally not that of the APSLG form Eq. (1). On the other hand it is clear (and we used this assumption in previous Sections) that for the ligands themselves (''organic'' part) the APSLG form of the trial wave function is a relevant approximation. We performed comparative study of electronic structures of simple amines and ethers on one hand and their cyclic polycounterparts on another hand by the semiempirical APSLG-MNDO method. The calculation results are given in Table . One can see that the relevant parameters of electronic structure (the bond orders and electron densities and the weights of the s-functions in the lone pairs) are fairly transferable from the low-molecular amines and ethers to their cyclic polyanalogs. Our calculations on cyclic chelating ligands have been performed at more or less arbitrary conformation of the molecule at hand. One can see that the dispersion of the values of all ESPs related to donor atoms entering the cyclic chelating ligands is always smaller than the dispersion of the same values in a series of ethers or amines ranging from water or ammonia to the corresponding alkyl di- or trisubstitutes, respectively. Thus the APSLG form (together with its semiempirical implementation) seems to be a relevant approximation for treating free chelating agents like crown ethers or cyclic polyamines.

This brings us to the situation we are already familiar with: to that one which requires different methods of description to be applied to different parts of a molecule. We faced such a necessity while describing the properties of the d-shell in transition metal complexes [] and while developing a general theory of the QM/MM junctions Separation,PCCP. In both cases the solution have been reached by employing the McWeeny group function approximation [] in specific physical conditions. The general method consists in (i) selecting the relevant Arai subspaces [], (ii) ascribing an adequate number of electrons to each of them, and (iii) by selecting an appropriate approximation for the electronic wave function in each of these subspaces. The physically substantiated picture of the metal-ligand bonding can be formalized by assuming the following form of the trial electronic wave function:

Y = FMLPÙFAPSLG

(34)

where Ù stands for the antisymmetrized product of the electronic functions; F_APSLG is the product of the bonding geminals Eq. ( 2) in the ''organic'' part of the complex; and F_MLP is a group function invoked to describe the metal ion with its closest environment. It is constructed in the Arai space spanned by the lone pair HOs defined by the Jacobi angles on the donor atoms and by the valence AOs of the metal ion. The number of electrons in F_MLP equals to that in all the lone pairs involved since the ion is assumed to provide only its vacant orbitals. According to Ref. [] the bonding geminals included in F_APSLG affect the effective Hamiltonian for the MLP Arai subspace only through the one-electron densities P_m^tt¢ residing in the HOs ascribed to the bonds. The same applies to the effect of the F_MLP function upon the effective Hamiltonians for the bonding geminals: only the one-electron densities in the lone pair HOs and in the metal AOs enter the expressions for the bond effective Hamiltonians. In that respect the form of the effective Hamiltonians for the bond geminals Eqs. (17), (18) does not change when the product of the lone pair functions Eq. (3) is replaced by F_MLP.

3.2.2  Electronic wave function near metal ion and the ESPs

Now let us address a possible approximation to be used for the F_MLP function. Two groups of concepts can be attracted to do the necessary choice. In contrast with the directionality and saturability characteristic for ''organic'' covalent bonds the bonds formed by metal ions do not posses these properties. Thus there is no need to invoke the HO formation on the metal ion. Also we notice that in the infinite separation limit the F_MLP wave function must flow to the antisymmetrized product of the lone pair geminals Eq. (3). The latter is in fact a single determinant function with all lone pair HOs doubly filled. Finally, we notice that important qualitative explanations concerning the structure of metal complexes have been given in Ref. [] with use of the self consistent field (SCF or single determinant) form of the electronic wave function in the Arai subspace spanned by the lone pair orbitals and metal s- and p-orbitals. With all these arguments we arrive to a conclusion that the single determinant (SCF) wave function is an appropriate form of the F_MLP function.

To get an economical description of the ESPs relevant to the F_MLP function we notice that the standard SCF wave function implies the MO expansion coefficients over the specified one-electron basis set to be the ESPs. This representation may be obtained by diagonalizing the matrix of the effective Fock operator (see below) in the specified Arai space. A more economical selection of the ESPs to be used in a QM/MM and eventually in a DMM-like description is possible. Indeed, the dimension of the Arai subspace we are interested in is N_M+N_LP, where N_M is the number of valence orbitals on the metal ion and N_LP is the number of lone pairs on the attached donor atoms. The number of the ESPs in the MO representation is (N_M+N_LP)², which is the number of the MO expansion coefficients. These coefficients are subject to (N_M+N_LP)(N_M+N_LP+1)/2 orthonormalization conditions, but, first, these conditions are very difficult to explicitly use (may be by introducing (N_M+N_LP)(N_M+N_LP-1)/2 Jacobi angles), and, second, even thus reduced number of parameters is superfluous. The reason is that the single determinant wave function is determined up to the subspace of the filled orbitals. Whatever unitary transformations applied separately to the filled and the empty one-electron states does not change the wave function. Since the numbers of parameters necessary to describe these irrelevant transformations equal, respectively, to N_M(N_M-1)/2 and N_LP(N_LP-1)/2 the minimal number of parameters required to describe the single determinant function with N_LP doubly filled and N_M empty orbitals is only N_MN_LP. This reduction is possible with use of the method given by Eqs. (21) - (24). Taking the product of the lone pair geminals (the ionic asymptotic limit) as a zero approximation for the F_MLP function we set the operator P₀ projecting to the subspace lone pair HOs (dim\limfuncImP₀ = N_LP) as a starting point for constructing the parameterization of the subspace of the filled orbitals according to Eq. (21). The projection operator P Eq. (21) is in its turn the one-electron density entering the effective bond Hamiltonians Eq.(18) for the bonding geminals and the semiempirical APSLG energy expression Eq.(5). Matrices V are obviously an N_LP×N_M matrices. They contain the relevant ESPs for the F_MLP function and ensure the correct asymptotic behavior: if a ligand goes to infinity the corresponding rows in the matrix V are set zero. One may check that in this case the corresponding off-diagonal matrix elements in the projection operator P remain zero and the diagonal ones remain unity as they were in the limiting projection operator P₀.

3.2.3  Effective Hamiltonian and DMM energy

Now let us turn to the contribution of the closest surrounding of the metal ion in the complex which consists of the metal centered AOs and the donor atoms' lone pair HOs to the energy. Notice first that the energy of the ''organic'' part of the complex is described by Eq. (5) with only variance that in the intraatomic intergeminal Coulomb terms (proportional to g_tmtk^¢) for the donor atoms whose lone pairs interact with the metal corresponding densities must be taken as diagonal matrix elements of the projection operator P instead of unity values characteristic for lone pairs in pure ''organic'' environment. The same values must be used in Eq. (8) for calculating the effective charges residing on the donor atoms.

The energy corresponding to the single determinant wave function with the occupied subspace \limfuncImP is given by []:

EMLP = 2( \limfuncSpheffP+\limfuncSpPS[P])

(35)

where h^eff the one-electron part of the effective Hamiltonian which contains also the mean Coulomb field induced by electrons from other electron groups (in our case - with those in the bond geminals of the ''organic'' part of the complex) in the MLP Arai subspace and S[P] is the average (mean-field) Coulomb electron-electron interaction of electrons described by the F_MLP function. The one-electron part of the effective group Hamiltonian can be presented as

heff = h0eff+h¢,

(36)

where h^¢ describes the one-electron transfers (resonance) between the lone pair HOs and the metal valence AOs. The effective operator h₀^eff can be defined as one commuting with the unperturbed projection operator P₀. For the metal AOs the matrix elements of h₀^eff are particularly simple:

( h0eff) mm = UmmM+ å B ¹ M,F  gMBQB+ å F  gMF æ è 2 å r Î APSLGÇF  Prr-ZF ö ø

(37)

m = s,p; and the first summation is extended to all atoms in the complex except the donor atoms involved in the resonance interaction with the metal atom (termed here as frontier atoms - F) and the metal atom itself. The second summation extends to those HOs on the frontier atoms which are responsible for the bond formation in the ''organic'' part of the complex. The diagonal matrix elements of h₀^eff for the m-th lone pair HOs residing on the F-th frontier atom have a form similar to that of the effective two-electron bond Hamiltonian Eq. (18) with obvious variance that the matrix elements of h₀^eff are ones of the one-electron operator:

( h0eff) lmlmF = Uml+ å B ¹ F  gMBQB+2 å tm1 Î APSLGÇF  glmtm1FPm1tt+ + å F¢  gFF¢ æ è 2 å r Î APSLGÇF¢  Prr-ZF¢ ö ø -gMFZM

(38)

By this the one-electron part of the effective Hamiltonian in the MLP Arai subspace is completely defined.

Now we address the average Coulomb interaction entering the energy expression. It has the standard SCF form for the metal orbitals:

Sss[P] = gssPss+2gsp å p  Ppp+2 å F  gMF å r Î MLPÇF  Prr; Spp[P] = gppPpp+2gspPss+2gpp¢ å p¢ ¹ p  Pp¢p¢+2 å F  gMF å r Î MLPÇF  Prr,

(39)

where the Coulomb parameters g_ss, g_sp, g_pp and g_pp^¢ relevant for the SCF description of the sp-shell of a metal atom have been introduced in Ref. []. The Coulomb interactions between electrons in the lone pars reduce to a standard contribution of the form:

Slmlm[P] = (lmlm|lmlm)FPmll+2 å tm1 Î MLPÇF  glmtm1FPm1tt+ +2gMF æ è Pss+ å p  Ppp ö ø +2 å F¢  gFF¢ å r Î MLPÇF¢  Prr

(40)

The latter expression flows to the same value as in the APSLG approximation for the free ligand since for the lone pairs the following (SCF-type) relation holds:

Gmll = (Pmll)2.

(41)

The SCF approximation implies also an appearance of the non-diagonal matrix elements of the Coulomb mean field. They have the form:

Smlm[P] = gMFPmlm

(42)

By this the energy operator for the electronic group describing the properties of the metal ion and its closest surrounding is defined.

Finally we address possible approximations to the energy of the metal ion with its closest surrounding on the basis of Eqs. (21), ( 35). Inserting these series in Eq. (35) and cutting the expansion at a particular overall order in V and V⁺ results in approximation to the energy as power series with respect to the set of the ESPs characteristic for the electron group of metal ion and its closest surrounding. The terms of odd overall order obviously correspond to the spin bond orders which appear between the lone pairs on the donor atoms and the metal ion. The necessity to take into account the variability of the metal-donor bond orders has been demonstrated recently in Ref. Sternbergetal in the context of an MM study. The terms of even overall order correspond to electron density transferred from the lone pairs to metal ion. While considering the dative bond from the point of view of a single donor atom we were restricting ourselves with the harmonic approximation relative to the corresponding parameter x appear only in the fourth order with respect to V. In the case of a metal complex it may be insufficient since the first terms responsible for the interactions between electrons transferred from the donor atoms to the metal. However, the Coulomb interactions between the donor atoms themselves vary already due to second order terms. The degree of this variation may be about a couple of electron volts by order of magnitude and should cause significant deviations from the points-on-a-sphere model and its modifications operating with constant interligand force fields.

4  Hybridization and dative bonding of the oxygen donor atom

In the previous Section the ESPs defined in the APSLG framework are determined by a variational procedure for the total energy Eq. (5). In the free ligand the shape and orientation of the system of HOs (hybridization tetrahedron) is determined by the geometry and interactions characteristic for the latter. It is only weakly perturbed by the resonance interaction with the empty (acceptor) orbitals. For the sake of simplicity modeling the acceptor was restricted by a single empty orbital of s-symmetry. This can be interpreted as a process of ''protonation'' or ''metallation'' of a donor atom. Forming additional (a superfluous from the point of view of usual valence) bond affects the shape of the hybridization tetrahedron. Considering the sp³-nitrogen atom [] reveals a correspondingly weak response of the nitrogen hybridization tetrahedron to the formation of an additional bond due to significant rigidity of the sp³-nitrogen hybridization tetrahedron. The formulae describing the reaction of the nitrogen hybridization tetrahedron are given in Ref. []. The case of doubly bonded oxygen donor atom remarkably differs from this simple picture. The covalent bonds formed by the oxygen atom fix only two of four HOs centered on it. Thus the linear dependency Eq.(14) does not suffice to fix the shape of the corresponding hybridization tetrahedron. To get around this uncertainty we consider here in detail the model ''water'' molecule i.e. of the oxygen atom with two hydrogen-like substituents which serve to represent usual covalent bonds. Then we consider what happens if an additional (dative) bonds are formed with an atom bearing single acceptor orbital of the s-character.

4.0.4  Model ''water'' molecule

The properties of the oxygen atom with two covalent bonds in the APSLG picture are determined by the interplay two energy contributions: (i) the one-center energy of the atom and (ii) the resonance energy of two covalent bonds it forms. The hybridization/density dependent part of the one-center energy Eq. (5) reads:

E¢ = 2 å m  ( sm2(Us-Up)+Up) dPm+ + å m  (C1+C2sm2+C3sm4)dGm+ + å k ¹ m  ( C4+C5[sm2+sk2]+C3sm2sk2) æ ç è dPkdPm+ 1 2 dPk+ 1 2 dPm ö ÷ ø .

(43)

where we assume that the corrections dP_m and dG_m are counted from their invariant values characteristic for covalently bonding geminals Eq.(15). In this case for the bonding geminals with m = 3,4 dP_m = dP and dG_m = dG correspond to their values in oxygen compounds. Any way, their variation is negligibly small as one can see from Table . We take them equal and do not vary in the following treatment for the sake of simplicity. For the lone pairs m = 1,2 and dP_m = [1/2] and dG_m = [3/4] (see above). The symmetry condition s₃ = s₄ = s also holds. One can check that the energy depends on the overall weight of the s-orbital in two lone pairs ( s₁²+s₂²) rather than on the specific distribution of the s-character between them. Adding the resonance energy of two covalent bonds we obtain the energy of a doubly bonded oxygen atom:

E = C3( 1 2 -dP)2[ 1-2s2] 2+ +[2[ Us-Up] ( 1 2 -dP)+2C3( 1 4 -dP2)+ +2C5( 3 4 -dP-dP2)+C2( 3 4 -dG)][1-2s2] + -2C3( dP+dP2-dG) s4+ -4 æ è 2bssOHs+bzsOH[( ® v   3  , ® e   3  )+( ® v   4  , ® e   4  )] ö ø POHrl

(44)

This expression contains all the molecular mechanics of a doubly bonded oxygen (at least as it comes from the semiempirical APSLG approximation). Let us consider characteristic features of this expression. As in the case of ammonia the minimum of the one-center energy corresponds to s = 0 which refers to no hybridization. The reason is that the expression Eq. ( 44) for the one-center energy is dominated by the U_s-U_p difference which amounts to -12 eV whereas the contributions from intraatomic Coulomb terms (both of second and fourth order in s_m's) entering with different signs (as in the case of nitrogen) covers the overall variation of no more than 2 eV. Taking into account that the range of variation of the variable s is restricted by the inequality 1-2s² ³ 0 due to normalization conditions allows to conclude that the one-center energy is a function with a single minimum at s = 0.

Inserting Eq. (11) to the expression for the resonance energy we find that the minimum of the resonance part alone is reached when:

bssOH(d ® w   b  , ® v   3  + ® v   4  )+bzsOHs(d ® w   b  , ® e   3  + ® e   4  ) = 0

The linear configuration with:

® v   3  = - ® v   4  ; ® e   3  = - ® e   4

is obviously a solution. It is the minimum for the resonance energy. Thus as in the case of an ''ammonia'' molecule modeling the sp³ hybridized nitrogen atom the actual (bent) form of our model ''water'' molecule is a result of an interplay between the one- and two-center contributions to the energy. The remarkable distinction from the nitrogen case is that the energy depends only on sum of the s-weights residing either in the bonding HOs or in the lone pairs, so that the weight 1-2s² of the s-orbital falling to two lone pairs on oxygen is arbitrarily distributed between them. Thus the shape of the hybridization tetrahedron on the oxygen atom remains undefined so far.

The form of the bonding HOs may be, nevertheless, specified on the basis of the orthonormality relations for the HOs which are consequences of the group SO(4) structure of the hybridization manifold. According to Eq. ( 12) the interhybrid angle for the bonding HOs is given by:

cosq34 = - s2 1-s2

For the symmetry reasons the p_y-AO of the oxygen atom does not contribute to the lone pair HOs and thus its weight is equally distributed between the bonding HOs:

v3y = -v4y = 1 Ö2

Since the overall p-weight residing in each bonding HO is 1-s² the x-components of the bonding HOs become:

v3x = v4x = -   æ  ú Ö 1 2 -s2

This comprises the description of the system of the HOs at the doubly bonded oxygen atom.

4.0.5  Donor-acceptor interactions of the model ''water'' molecule.

Now let us consider what happens when an additional ''dative'' bond is formed with an atom containing for the sake of simplicity only one empty s-orbital used for bonding. The linear response approximation previously constantly employed in the DMM framework cannot be employed in the present case since the matrix of the energy second derivatives with respect to variations d[(w)\vec] _b and d[(w)\vec] _l is degenerate (and thus cannot be inverted) due to above mentioned invariance of the energy with respect to the deformations of the oxygen hybridization tetrahedron between approximate sp² and sp³ hybridization of the lone pairs. For that reason we try to extract the information on the shape of the oxygen hybridization tetrahedron with an extra dative bond from the structure of the SO(4) hybridization manifold. This is an analog for the degenerate perturbation theory for the energy considered as a function of the set of hybridization parameters [(w)\vec] _b and [(w)\vec] _l.

To start with it we assume that the HO with m = 1 will be used for the dative bond. For it we use the estimates Eq.(30) for the density ESPs. This may be termed as a ''harmonic'' approximation in terms of the x ESP. The energy of the dative bond is then given by

EOAbond » -2Ö2bOAx-gOAx2

(45)

With these definitions we get a correction to the one-center energy of the ''water'' oxygen atom due to partial electron density transfer to the acceptor orbital:

DE¢ = = -(s12(Us-Up)+Up+(C1+C2s12+C3s14)- -( C4+C5[1-2s2]+C3s12[1-2s2-s12]) )x2+ +2( C4+C5[s12+s2]+C3s12s2) æ ç è dP+ 1 2 ö ÷ ø x2.

(46)

This contribution lifts the degeneracy of the one-center energy with respect to the distribution of the s-weight between the two lone pairs. This contribution describes the energy required to extract an amount of electron density proportional to x² from the oxygen lone pair with the weight of the s-function equal to s₁². This energy is the larger the larger is the s-weight of the HO involved. This is in agreement with the estimates of the ionization potential of the water molecule performed both in the semiempirical APSLG approximation Ref. [] and with independent ones.

The formulae Eqs. (44) - (46) present together a specific QM/MM force fields for the dative bond formed by the doubly bonded oxygen atom. It is so since the cited equations represent the energy components in terms of the parameters of the semiempirical QM Hamiltonian and of the ESPs characterizing the covalent and dative bond on one hand and the shape and orientation of the hybridization tetrahedron residing on the oxygen atom on the other hand. On a recipe level the sum of Eqs. (45), (46) must be added to the total energy and the latter must be optimized also with respect to x and s₁ as well as to all other ESPs at each value of the geometry parameters. A simplified treatment with somehow fixed values of dP and s is also possible. A remarkable feature of this approach is that it remains valid both at very large and shorter separations between the donor and acceptor atoms which allows to cover uniformly the regions normally treated by different methods: by QM at a shorter distances and by standard MM at longer ones.

Now we are equipped to study the shape of the hybridization tetrahedron on the oxygen donor atom. The structure of the SO(4) hybridization manifold does not pose enough restriction on its flexibility. We have to remind in this context that in the case of quadruply bonded carbon or triply bonded nitrogen atoms namely the structure of the hybridization manifold fixed the tetrahedral form of the model ''methane'' or ''ammonia'' molecules through the linear dependence relation Eq. (14). In the case of the doubly bonded oxygen the two defined vector parts ([(v)\vec]₃ and [(v)\vec]₄) do not suffice to determine the other two. We assume that the perturbation incurred by the dative bond formation does not decrease the overall s-weight (2s²) of the covalent bonding HOs since it would results in a too large energy increase of the ''water'' molecule. The energy of the latter is, however, independent on the actual value of s₁ which controls the distribution of the s-weight between the lone pairs. The latter can access only a restricted range of values:

0 £ s12 £ 1-2s2;s22 = 1-2s2-s12

The limiting values of s₁² correspond either to pure p-character of the dative bond HO (s₁² = 0) pointing aligned to the normal of the ''water'' molecular plane or to almost sp² HO directed along the C₂ axis lended for the bond formation. The hybridization manifold structure uniquely determines the direction of the dative bonding HO for each specific value of s₁. According to Eq. (12) we have:

cosq13 = cosq14 = - s1s   _____ Ö1-s12     ____ Ö1-s2

With the value of s fixed this defines the projection of the dative bonding HO (m = 1) on the C₂ axis (x); the rest of the p-weight comes from the z-component of the HO:

v1x = - Ö2s1s   _____ Ö1-2s2 ;v1y = 0;v1z =   æ  ú Ö 1- s12 1-2s2   .

The latter result means that the bonding HO stays in the mirror plane of the ''water'' molecule. For the acceptor ion the going out of this plane results in the corresponding restoring force. Thus obtained values of the components of the vector part of the dative bonding HO must be inserted in the expression Eq. (33) for the resonance integral. Analysis of the one-center term Eq. (46) describing the energy response of the oxygen hybridization tetrahedron to the formation of the dative bond leads to the conclusion that for the dative bond the energy minimum is reached if the pure p-HO is involved in its formation. This result can be easily understood since the resistance to the bond formation (the coefficient at the x²) is larger for the larger s-contribution to the bonding HO. We arrive to an interesting situation different from that for the nitrogen dative bond. In the case of oxygen we expect that the bond-related force field opposes the electrostatic (ion-dipole) forces which tend to place the acceptor atom on the C₂ axis which coincides with the direction of the donor effective dipole moment. The authors of Ref. HayAllinger report a possibility of such competition on the basis of analysis of the structures of crown-ether complexes of alkali and alkali earth ions. It turns out that the singly charged alkali ions much easier acquire a nonplanar coordination geometry at least with one of the ether oxygen atoms in the crown ether complexes than the doubly charged alkali earth ions. For the latter the planar trigonal geometry of the ether donor oxygen is preferred. This is precisely the trend one would expect on the basis of relative strength of the ion-dipole interactions of single and double charged ions.

In order to numerically test the above evaluations we performed a series of APSLG-MNDO calculations on a simple model of a complex of Li⁺ with H₂O where the cation has been represented by a single s-orbital with the standard MNDO parameterization for lithium. The calculation have been performed for the Li-O separation of 2.14 Å characteristic for Li⁺ complexes with ethers []. Results are presented in Table . In agreement with the above estimates we found first of all that the approximations of the ionic limit are valid for the ion-molecule interaction. In all cases the equilibrium value of x does not exceed 0.3 which can be shown to be a safe estimate for the validity of the ionic limit expansions. This ESP reaches its maximum for the "p"-coordination of the lithium ion to the water molecule. Nevertheless, the energy of this configuration is maximal. By contrast the minimal energy is reached for the planar configuration at the oxygen atom where the Li-O bond order is also minimal. This demonstrates the dominance of the electrostatic forces in shaping the molecular geometry of the model complexes. On the other hand the shape of the system of HO residing on the oxygen atom pretty much deviates from one expected on purely geometrical grounds. Almost true coincidence is observed for the "p"-coordination of the lithium ion. In this case the bonding HO is almost perpendicular to the water plane though the s-weight in it is not negligible: ca. 8%. At the same time in the configuration corresponding to the energy minimum the angle between the Li-O bond and the corresponding HO is also noticeable.

To conclude this Section we notice that the described discrepance between the molecular geometry at the oxygen donor atom and the shape of its hybridization tetrahedron is characteristic only for the ionic limit of the dative bond. If the additional bond reaches the covalent regime characteristic for example for real water protonation (formation of the H₃O⁺ cation) there is no such uncertainty and as in the isoelectronic case of ammonia the pyramidal C_3v geometry is with no doubt the equilibrium one and the misalignment between the bonds and HO does not exceed a couple of degrees [].

5  Conclusion

In the present paper we analyzed the behavior of the APSLG approximation at the frontier of its applicability area. Being originally designed for treating the systems with well defined localized two-center bonds it is employed here for analysis of dative (donor-acceptor) bonds and coordination compounds of metal ions. Two major results are acquired on this way. First, we developed a DMM description for the dative bonds formed by ether oxygen atom. It turned out that the DMM of the dative bonds differs from that of the usual covalent bonds in that respect that the ESPs corresponding to the bond orders are strongly distance dependent. Such a situation is described in the MM context (see e.g. Ref. []) by referring to the Pauling's bond-length-bond-order logarithmic relation []. Here we propose a sequential description for such a situation based on the standard QM technique. Also this treatment allowed to analyse the known flexibility of the coordination mode of the ether oxygen to acceptors and to rationale the observed dependence of the coordination trends on the charge of the metal cation. Second, we proposed to describe the metal ion and its closest surrounding with use of the McWeeny's parameterization for the occupied states projection operator in the space of one-electron functions spanned by the lone-pair HOs and the metal ion vacant valence orbitals. Within that sort of description one may analyze the sources and expected relative importance of different empirical force fields which appear in the The description in terms of the projection operator in the space of one-electronic states, of course, occupies a border position between the MM and QM ways of describing molecular structure.

One also has to notice that even in the ''organic'' MM realm the explicit reference to details of electronic structure (particularly when it concerns the redistribution of effective charges) is considered to be acceptable due to Gasteiger (see [] and reference therein). So the border between two types of description becomes more and more smeared with passage of time. On the other hand addressing the electronic structure in that or another way becomes even better substantiated when for the object at hand one may expect even much more significant charge redistribution characteristic for metal complexes. In any case the proposed treatment of the metal ion with its closest surrounding can be treated as an example of application of the general QM/MM methodology [,] to this specific problem.

6  Appendix

The Schrödinger equation in terms of projection operator Eq. ( 27) gives a start for the perturbation expansion. Indeed, for the projection operator Eq.(26)

P = P0+P¢

where

P¢ = 1 1+x2+y2 æ ç ç ç è -x2-y2 x y x x2 xy y xy y2 ö ÷ ÷ ÷ ø

we have

(H0+H¢)( P0+P¢) = (P0+P¢)(H0+H¢) H0P0 = P0H0

which can be rewritten as:

[H0,P¢] = [P0,H¢]+[P¢,H¢]

with square brackets standing for commutators. This can be explicitly written as:

æ ç ç ç è 0 (a-b)x (a-c)y -(a-b)x 0 (b-c)xy -(a-c)y -(b-c)xy 0 ö ÷ ÷ ÷ ø = d æ ç ç ç è 0 1 0 -1 0 0 0 0 0 ö ÷ ÷ ÷ ø + +d æ ç ç ç è 0 -x2+y -xy+x x2-y 0 -y-y2+x2 -x+xy y+y2-x2 0 ö ÷ ÷ ÷ ø

which gives the system of equations Eq. (28) for x and y by equating the left- and right-hand matrix elements.

7  Acknowledgments

The work is performed with partial financial support through the grant No 6-120 dispatched by the Young Researchers Commission of RAS. The author is thankful to Mr. A.M. Tokmachev and Mr. M.B. Darhovskii for their help in performing this work.

References

[]

M. Jones, Organic Chemistry, Northon Press, New York, 2000.

[]

U. Burkert, N.L. Allinger, Molecular Mechanics. ACS, Washington, D.C. 1982.

[]

T.A. Halgren. J Comp Chem 17 (1996) 490 - 519; ibid. 17 (1996) 520 - 552; ibid. 17 (1996) 553 - 586; T.A. Halgren, R.B. Nachbar. ibid. 17 (1996) 587 - 615; ibid. 17 (1996) 516 - 641.

[]

N.L. Allinger, K. Chen, J.-H. Lii. J Comp Chem 17 (1996) 642 - 668; N. Nevins, K. Chen, N.L. Allinger. ibid. 17 (1996) 669 - 694; N. Nevins, J.-H. Lii, N.L. Allinger. ibid. 17 (1996) 695 - 729; N. Nevins, N.L. Allinger. ibid. 17 (1996) 730 - 746; N.L. Allinger, K. Chen, J.A. Katzenellenbogen, S.R. Wilson, G.M. Anstead. ibid. 17 (1996) 747 - 755.

[]

P. Comba, T.W. Hambley. Molecular modeling of inorganic compounds. VCH, New York, 1995.

[]

V. Burton, R. Deeth, C. Kemp, P. Gilbert. J Am Chem Soc 117 (1995) 8407; R. Deeth, V. Paget. J Chem Soc Dalton Trans (1997) 537.

[]

I. V. Pletnev, V. L. Melnikov. Koord Khim 23 (1997) 205 [in Russian]; Russ Chem Bull 7 (1997) 1278.

[]

C. Millot, Dehez. J Comp Meth Sci Eng 2 (2002) 451 - 456.

[]

B.P. Hay, L. Yang, J.-H. Lii, N.L. Allinger. J Mol Struct THEOCHEM 428 (1998) 203 - 219.

[]

J. Golebiowski, V. Lamare, M.T.C. Martins-Costa, C. Millot, M.F. Ruiz-López. Chem Phys 272 (2001) 47 - 59.

[]

M. Möllhoff and U. Sternberg J Mol Model 7 (2001) 90-102.

[]

U. Sternberg, F.-T. Koch, M.Bräuer, M. Kunert, and E. Anders. J Mol Model 7 (2001) 54-64.

[]

A.V. Soudackov, A.L. Tchougréeff, and I.A. Misurkin. Russ J Phys Chem (1996); A.V. Sinitsky, M.B. Darkhovskii, A.L. Tchougréeff, and I.A. Misurkin. Int J Quant Chem 88 (2002) 370-379.

[]

M. Fuxreiter, Ö. Farkas, and G. Náray-Szabó. Protein Engineering (1995) 8, 925 - 933.

[]

R. Cini, D.G. Musaev, L.G. Marzilli, K. Morokuma. J Mol Struct THEOCHEM 392 (1997) 55 - 64.

[]

O.L. Malta, H.J. Batista, L.D. Carlos. Chem Phys 282 (2002) 21 - 30.

[]

A.L. Tchougréeff, A.M. Tokmachev. J Comp Meth Sci Eng 2 (2002) 309 - 314.

[]

A.L. Tchougréeff. J Mol Struct THEOCHEM Accepted.

[]

A.M. Tokmachev, A.L. Tchougréeff. Int J Quant Chem (2003) Accepted.

[]

M.B. Darhovskii, A.L. Tchougréeff. Russ J Phys Chem 75 (2000) 112; M.B. Darhovskii, M.G. Razumov, I.V. Pletnev, A.L. Tchougréeff. Int J Quant Chem 88 (2002) 588 - 605.

[]

A.M. Tokmachev, A.L. Tchougréeff. Russ J Phys Chem 73 (1999) 259.

[]

A.M. Tokmachev, A.L. Tchougréeff. J Comp Chem 22 (2001) 752-764.

[]

A.M. Tokmachev, A.L. Tchougréeff. J Phys Chem A 107 (2003) 358 -365.

[]

P. Surján. Top Curr Chem 203 (1999) 64 - 88.

[]

R.C. Bingham, M.J.S. Dewar, D.H. Lo. J Am Chem Soc 97 (1975) 1302, 1307.

[]

J.A. Pople, D.L. Beveridge. Approximate Molecular Orbital Theory. McGraw-Hill, NY, 1970.

[]

S.L. Altmann, Rotations, quaternions and double groups. Oxford Scientific Publications, Clarendon Press, Oxford University Press: NY 1986; A.V. Berezin, Yu.A. Kurochkin, and E.A. Tolkachev, Quaternions in relativistic physics. Nauka i Tehnika, Minsk, 1989 [in Russian].

[]

S. Shaik, A. Shurki. Angew Chem Int Ed 38 (1999) 586 - 625.

[]

R. McWeeny. Rev Mod Phys 35 (1960) 335.

[]

V.I. Pupyshev, N.F. Stepanov. in Modern Problems of Physical Chemistry. v. 11, Moscow State University Publishers, Moscow, 1979 [in Russian].

[]

A.M. Tokmachev, A.L. Tchougréeff, and I.A. Misurkin. Int J Quant Chem 85 (2001) 109-117.

[]

I.B. Bersuker. Electronic structure and properties of transition metal complexes; Nauka: Moscow, 1986 [in Russian].

[]

B. Hay. Coord Chem Rev 126 (1993) 177.

[]

A.V. Soudackov, A.L. Tchougréeff, I.A. Misurkin. Theor Chim Acta 83 (1992) 389 - 416.

[]

A.L. Tchougréeff. Chem Phys Reports 16 (1997) 1035.

[]

A.L. Tchougréeff. Phys Chem Chem Phys 1 (1999) 1051 - 1060; A.M. Tokmachev, A.L. Tchougréeff, and I.A. Misurkin. J Mol Struct THEOCHEM 506 (2000) 17 - 34.

[]

R. McWeeny, B.T. Sutcliffe. Methods of molecular quantum mechanics. AP, London (1969).

[]

A.A. Levin, P.N. Dyachkov. Electronic structure and transformations of heteroligand molecules, Moscow, Nauka, 1990 [in Russian]; A.A. Levin, P.N. Dyachkov. Heteroligand molecular systems: bonding, shapes and isomer stabilities. Taylor and Francis. NY, London. 2002.

[]

L. Di Sipio, E. Tondello, G. De Michelis, L. Oleari. Chem Phys Lett 11 (1971) 287-289.

[]

A.M. Tokmachev, A.L. Tchougréeff. To be published.

[]

L. Pauling. The Nature of the Chemical Bond, Cornell University Press: Ithaca, NY, 1960.

[]

W. J. Mortier, K. Van Genechten, J. Gasteiger. J Am Chem Soc 107 (1985) 829-835; J. Gasteiger, H. Saller. Angew Chem Int Ed Engl 24 (1985) 687-689; L.G. Hammarstrom, T. Liljefors, J. Gasteiger. J Comp Chem 9 (1988) 424-440.