A.L. Tchougréeff and A.M. Tokmachev
Karpov Institute of Physical Chemistry
10 Vorontsovo pole, Moscow, 105064, Russia
and
Center for Computational Chemistry at
the Keldysh Institute for Applied Mathematics of RAS
4 Miusskaya pl., Moscow, 125047, Russia

Deductive Molecular Mechanics of sp³ Carbon Atom

Abstract

Introduction

Molecular mechanics (MM) [] is normally considered as a purely empirical scheme based on representation of the total energy as a sum of contributions (force fields) which are explicit functions of molecular geometry. The form of these force fields is usually taken as one or two first non-zero terms in the Taylor expansion for the energy. The parameters of the force fields are adjusted to reproduce the experimental (or sometimes quantum chemical) data on the heats of formation, equilibrium geometries, and sometimes also vibration frequencies etc. By construction the MM methodology is inapplicable to essentially quantum situations and to the molecules where electronic delocalization or/and correlation are important. Nevertheless the MM remains extremely useful tool of computational chemistry due to very low computational costs and high quality of results on molecular geometry and heats of formation which typically exceeds even that obtained by the high-level quantum chemical methods.

Despite its wide applicability MM approach lacks any reliable theoretical substantiation. The necessity of such substantiation is caused not only by purely theoretical discomfort but also by strong methodological demand since the hybrid quantum mechanical/molecular mechanical (QM/MM) schemes combining the advantages of the QM and MM approaches have been actively developing during last decade. Definition of junction between quantum and classical subsystems on the solid theoretical basis requires some QM description underlying the MM one []. The problem of substantiation of the MM is also close to that of theoretical understanding of organic stereochemistry. Explanatory concepts in this area [] are somehow guessed but are not derived from any reliable QM description of molecular electronic structure and thus deserve theoretical substantiation. The same is true for additive systematics for the heats of formation of organic molecules [].

There were some attempts reported in the literature to understand the applicability of classical additive schemes from the QM principles. The authors of Ref. [] applied the PCILO method [] as the underlying QM scheme for additive systematics. At the same time the non-variational nature of local one-electron states and amplitudes of two-electron configurations employed in Ref. [] has not allowed to derive the analytical form of the main force fields - those of stretching and bending. Recently, we proposed an alternative approach leading to a generic MM scheme derived from the QM description []. It is based on a semiempirical implementation of the trial wave function taken in the form of antisymmetrized product of strictly localized geminals (APSLG) []. This implementation possesses some attractive features: correct behaviour of the wave function at all intrabond separations; variational determination of basis one-electron states; reliable results on heats of formation and molecular structures; linear scaling of computational costs. It allowed us [] to construct a series of reliable non-iterative (MM-like) schemes. At the same time these schemes are purely numerical. Here we present an abridged analytical derivation of generic MM scheme for the case of sp³ carbon atom with special attention paid to explicit form of the force fields and to demonstration of their sources.

The paper is organized as follows: in the next Section we briefly review the APSLG implementation used; then the bond energy functions are considered; in the next Section we thoroughly analyze the bending force field as derived from the APSLG energy functional; then the consequences of adjustment the one-electron states due to geometry variation are considered and, finally, the numerical results are given and discussed.

APSLG-MINDO/3 method

The APSLG trial wave function has the form:

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

(0)

where each bond geminal g_m⁺ is constructed as a linear combination of three singlet two-electron configurations (two ionic and one covalent - or Heitler-London) with variable amplitudes:

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

(0)

subject to normalization condition imposed u_m²+v_m²+2w_m² = 1. In the case of lone pair geminal only one contribution survives:

gm+ = rma+rmb+.

(0)

The basis one-electron states | r_m ñ and |l_m ñ assigned to the m-th geminal are taken as strictly local hybrid orbitals (HOs), i.e. obtained by a unitary SO(4) transformation of four atomic orbitals for each ''heavy'' (non-hydrogen) atom:

tms+ = å i Î A  hmiAais+.

(0)

The 4×4 hybridization matrices h^A are determined variationally. Geminal amplitudes Eq. (2) are in their turn determined by diagonalizing the effective bond Hamiltonians.

The electronic Hamiltonian used is that of the MINDO/3 type [] with resonance parameters b^AB slightly re-adjusted []. The core-core interaction in the MINDO/3 scheme is not purely Coulomb one but is modified by a short-range repulsion term:

Enn = 1 2 å a ¹ b  Cab, Cab = Za Zb (gab+Dab),

(0)

where g_ab is the two-center Coulomb integral.

Bond energy, equilibrium bond length and bond stretching constant

The APSLG-MINDO/3 energy can be re-written in the form close to that of the MM (with different bonding and non-bonding contributions, 1-3 bond interaction contributions) as a sum of increments:

EA = å tm Î A  [2UmtPmtt+(tmtm|tmtm)TmGmtt]+2 å \Sb tktm¢ Î A k < m\endSb gtktm¢TkPkttPmt¢t¢, ERmLmbond = 2gRmLm[Gmrl-2PmrrPmll]-4brmlmRmLmPmrl, EABnonbond = QAQBgAB.

(0)

Here Umt is the matrix element of attraction of electron to its own core, brmlmRmLm is the intrabond resonance integral, QA is the Mulliken charge on the atom A, and gtktm¢Tk = 2(tktk|tm¢tm¢)Tk-(tktm¢|tm¢tk)Tk is the reduced Coulomb integral.

The energy expression Eq. (6) depends on the electronic structure parameters (ESPs) of two types (i) the elements of density matrices

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)

where t and t^¢ are either r or l, and (ii) the matrices transforming the AO's to the | t ñ one-electron states which describe the hybridization. (The energy depends on the latter through molecular integrals). It was shown in Ref. [] by analysis of the bond related contributions to the energy that there exist two parameters (functions of bond interatomic separation) for each bond z_m^-1 and m_m describing, respectively, intrabond correlation and bond asymmetry. It was shown that for all chemical bonds these parameters are small. Therefore, it is possible to describe transferability of the elements of density matrices by expression of them as power series on these parameters. It was shown [] that the contribution of the zeroth order in m_m can be written as:

Gmtt¢ = 1 4 æ ç è 1-tt¢ 1 G(zm) ö ÷ ø , Pmtt = 1 2 , Pmrl = zm 2G(zm) ,

(0)

where G(z_m) = Ö{1+z_m²}. At the equilibrium intrabond separations the limit z_m® ¥ is valid. In this limit

Gmtt¢ = 1 4 æ ç è 1-tt¢ 1 zm ö ÷ ø , Pmrl = 1 2 æ ç è 1- 1 32zm2 ö ÷ ø .

(0)

Taking the contributions of up to second order in m_m the density matrix elements become slightly modified []:

Gmrl = Gm0rl é ê ë 1+mm2 (2G(zm)+1)(1-G(zm)) 2(G(zm)+1)2 ù ú û , Pmtt = Pm0tt é ê ë 1+tmm G(zm)-1 G(zm)+1 ù ú û , Pmrl = Pm0rl é ê ë 1+mm2 2G(zm)+1-G2(zm) 2(G(zm)+1)2 ù ú û ,

(0)

where subscript 0 corresponds to the estimates by Eq. (8). The formulae Eqs. (9) and (10) allow to conclude that the bond order (2P_m^rl) is transferable up to second order with respect to both parameters z_m^-1 and m_m; the bond covalency (2G_m^rl) is transferable up to second order with respect to m_m and up to first order with respect to z_m^-1; the bond polarity (P_m^rr-P_m^ll) is transferable up to first order with respect to both z_m^-1 and m_m. The transferability of bond orders will be used in the following considerations. The limit z_m® ¥, m_m® 0 can be considered as a good approximation in the case of non-polar bonds (C-H, N-H, C-C, N-C etc.) for the respective bond lengths close to equilibrium. It gives the following estimates:

Pmtt¢ = 1 2 , Gmtt¢ = 1 4 .

(0)

This result can be also obtained from the SCF consideration of symmetric bond. Such picture can be termed as FA (fixed geminal amplitudes) to contrast it to the more refined picture of Eqs. (8), ( 9), and (10) which can be termed as TA (tuned geminal amplitudes). These notations are in the spirit of those of Ref. QMMMIJQC.

From the general APSLG energy expression the equilibrium bond length can be deduced by equating to zero the first derivative of the bond-related energy. We consider quite reliable case of non-polar bond Eq. (8). The derivatives of the gemina related ESPs exactly compensate each other and the bond energy minimum equation is:

¶E ¶q = ZRmZLm ¶DRmLm ¶q -2 zm G(zm) ¶brmlmRmLm ¶q - 1 2 æ ç è 1- 1 G(zm) ö ÷ ø ¶gRmLm ¶q = 0.

(0)

In the limit z_m® ¥ we recover the equilibrium geometry condition for the FA picture. Note that if q is the intrabond distance the partial derivatives of b_rmlm^R_mL_m and g_RmLm are positive which leads to attraction of the bonded atoms. The short-range repulsion is only due to special form of the core-core interaction.

The same concepts can be used to determine the elasticity constant for the bond stretching by taking the second derivative of the energy with respect to the bond length. In the FA picture we get:

kRmLm = æ ç è ZRmZLm d2DRmLm drRmLm2 -2 ¶2brmlmRmLm ¶rRmLm2 - 1 2 d2gRmLm drRmLm2 ö ÷ ø rRmLm0  .

(0)

Numerical results on this constant will be given below.

Hybridization and bending force field

This Section describes in details the role of hybridization in construction of the MM force fields. In general we consider a kind of a linear response relation between variations of nuclear coordinates and ESPs characterizing HOs. Its general form can be readily obtained by expanding the energy up to the second order with respect to both atomic coordinates q and the ESPs x with subsequent minimization:

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

(0)

It leads to the second order energy expression:

E = E0+ 1 2 ( q-q0| ÑqÑqE| q-q0) - - 1 2 ( q-q0| ÑqÑxE( ÑxÑxE) -1ÑxÑqE| q-q0) .

(0)

To derive the bending force field we need to know the main source of angular dependence of the energy. In the APSLG approximation it is obviously determined by the structure of HOs. In Ref. [] the one-center hybridization dependent contribution to the energy has been analyzed. It turned out that the one-center energy is hybridization dependent only if subtle polarization and correlation effects are taken into account. This contradicts to the fact that the correct hybridization can be reproduced numerically even by methods without intraatomic electron-electron repulsion []. Therefore, we must admit that the main source of hybridization is the resonance energy, which is in agreement with earlier concepts in this area [].

Hybridization description

To consider the variations of the resonance energy we studied first the mathematical structure of hybridization. Each of HOs centered at a given atom is a linear combitation of the s- and p-AOs. One can note, however, that when a molecule is rotated as a whole the expansion coefficients at the s-AO are invariant, whereas the coefficiens at the p-AOs transform as if they were components of a 3-vector. Thus the HOs can be represented as object combining a scalar and a vector parts: (s_m,[(v)\vec]_m) []. This mathematical object is a normalized quaternion.

The ends of vectors [(v)\vec]_m form a hybridization tetrahedron containing full information about hybridization. As it is mentioned above the specific hybridization is given by an SO(4) matrix. These matrices form the six-parametric group and the most direct way to parameterize the latter is to represent them as products of six Jacobi matrices performing rotations in two-dimensional subspaces by the angles w_sx^A, w_sy^A, w_sz^A, w_yz^A, w_xz^A, and w_xy^A (subscripts refer to the two-dimensional subspace). The first triple of angles corresponds to formation of hybridization tetrahedron (a pseudorotation) while the second triple of angles corresponds to rotation of the hybridization tetrahedron as a whole (a quasirotation, where prefix quasi signifies that no physical body actually rotates). The SO(4) group generates all possible hybridizations thus being a dynamical group for the system of all possible HOs at a given atom [].

The above parametrization of the SO(4) group by the Jacobi matrices is inconvenient for analytical considerations since the matrix elements are clumsy combinations of trigonometrical functions. Therefore, we propose more workable parameterization of the SO(4) group by a pair of normalized quaternions similar to the quaternion parameterization of the SO(3) group of 3-dimensional rotations []. The quaternions (q and p) correspond to fictitious rotations defined by triples of angles:

wg± = eabgwab±wsg.

(0)

As in the case of quaternion representation of rotations the components of the quaternions at hand are given by:

q0 = cos w+ 2 ,q1 = w+x w+ sin w+ 2 ,q2 = w+y w+ sin w+ 2 ,q3 = w+z w+ sin w+ 2 , p0 = cos w- 2 ,p1 = w-x w- sin w- 2 ,p2 = w-y w- sin w- 2 ,p3 = w-z w- sin w- 2 .

(0)

The SO(4) matrix in terms of the quaternion components (the derivation will be published elsewhere) has the form:

h = æ ç ç ç ç è q0p0+q1p1+q2p2+q3p3 q0p1-q1p0-q2p3+q3p2 -q0p1+q1p0-q2p3+q3p2 q0p0+q1p1-q2p2-q3p3 -q0p2+q1p3+q2p0-q3p1 q0p3+q1p2+q2p1+q3p0 -q0p3-q1p2+q2p1+q3p0 -q0p2+q1p3-q2p0+q3p1 q0p2+q1p3-q2p0-q3p1 q0p3-q1p2+q2p1-q3p0 -q0p3+q1p2+q2p1-q3p0 q0p2+q1p3+q2p0+q3p1 q0p0-q1p1+q2p2-q3p3 -q0p1-q1p0+q2p3+q3p2 q0p1+q1p0+q2p3+q3p2 q0p0-q1p1-q2p2+q3p3 ö ÷ ÷ ÷ ÷ ø .

(0)

Due to algebraic structure the SO(4) matrix H close to the unity matrix also represents small variation of HOs in a vicinity of a given set of HOs. The first order contribution for example is:

d(1)H = æ ç ç ç ç è 0 p1-q1 p2-q2 p3-q3 q1-p1 0 -q3-p3 q2+p2 q2-p2 q3+p3 0 -q1-p1 q3-p3 -q2-p2 q1+p1 0 ö ÷ ÷ ÷ ÷ ø ,

(0)

and analogously the second order contribution can be extracted. The first order correction to the HOs in the quaternion form (s,[(v)\vec]) which appear when small quasi- and pseudorotations d[(w)\vec] _l and d[(w)\vec] _b are applied to the system of HOs at a given atom has a simple form:

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

(0)

where × stands for the vector product of 3-vectors.

Valence angles and bending energy

It is known for almost half a century [] that in general terms stereochemistry (form of the coordination polyhedron) is determined by relation between the two-center energy which favours population of excited and ionized states of an atom under consideration and the excitation and ionization energies themselves which tend to keep an atom in its ground (unhybridized) neutral state. In the standard MM framework it is assumed that there is a characteristic angle to be maintained for each pair of incident chemical bonds (triple of atoms). These angles are chosen as totally independent parameters. At the same time it is known that only five of six angles (in the case of substituted carbon) are really independent []. This contradiction is important and should be taken into account when one tries to derive the MM scheme. Also the parameters characterizing the bending energy in the standard form

1 2 kABC(q-q0)2

(0)

are considered to be independent parameters. In the framework of the APSLG derived MM-like scheme the system of HOs at a given atom (hybridization tetrahedron) has more rigid structure than the coordination tetrahedron because the form of hybridization tetrahedron is determined by only three parameters (pseudorotation angles [(w)\vec] _b). It leads to important purely algebraic consequences which will be considered below.

We want to demonstrate the general correspondence between the structure of the APSLG wave function and the corresponding energy functional on one hand and the MM force fields on the other hand. To do this we consider a characteristic example - the sp³ carbon atom. For the sake of simplicity we give explicit formulae only for the simple case of hydride (methane) while the generalization of the results to more complex cases is straightforward and will be published elsewhere. We characterize the stereochemistry of the carbon atom (the coordination tetrahedron) by a set of four unit vectors [(e)\vec]_RmLm which have the carbon atom as their origin and point to the hydrogen atoms. In the diatomic coordinate frame (DCF) with the z axis directed along the [(e)\vec]_RmLm vector the resonance integrals can be written as:

brmlmRmLm = bssRmLmhmsRmhmsLm+bzsRmLmhmzRmhmsLm

(0)

or

brmlmCH = bssCHsmC+bzsCHvmzC,

(0)

where

vmzTm = ( ® v   Tm m  , ® e   RmLm  )

(0)

(in the case of hydrogen atom the contribution of the s-orbital always equals to unity).

Now we obtain the general equilibrium conditions for the hybridization tetrahedron. Taking the first order variation of the resonance integrals Eq. (23) and making use of formulae Eq. (20) we get:

d(1)brmlmCH = -bssCH(d ® w   b  , ® v   C m  )+ +bzsCH[smC(d ® w   b  , ® e   RmLm  )+(d ® w   l  × ® v   , ® e   RmLm  )].

(0)

From this equation the expression for pseudotorque [(N)\vec] and quasitorque [(K)\vec] can be easily obtained as a sum of contributions to the resonance energy for all four bonds which are linear in d[(w)\vec] _b and d[(w)\vec] _l , respectively:

® N   = 4 å m  Pmrl{bssCH ® v   C m  -bzsCHsmC ® e   RmLm  }, ® K   = 4 å m  PmrlbzsCH ® e   RmLm  × ® v   Rm m  .

(0)

The equilibrium condition for the HOs is then:

® N   = ® K   = ® 0   .

(0)

These conditions govern both the form and orientation of the hybridization tetrahedron. The solution of these equations in the case of symmetric hydride is obvious because if we direct all vectors [(v)\vec]_m^C along vectors [(e)\vec]_RmLm both equations become satisfied. We note that in the first equation two different contributions (from HOs [(v)\vec]_m^C and vectors [(e)\vec]_RmLm) vanish independently and that the second contribution is always zero when directions of HOs coincide with directions of bonds. Tetrahedral coordination and hybridization correspond to the energy minimum in the symmetric case. Corrections to this characteristic case can be meanwhile considered as perturbations.

Let us consider now a situation when parameters [(w)\vec] _b are fixed by some reasons. It can be termed as FO (fixed orbitals) picture following the notations of Ref. [] because the form of hybridization tetrahedron is fixed in this case. It worths mentioning that if the form of the hybridization is fixed ([(w)\vec] _b), the one-center contributions to the energy are precisely [(w)\vec] _l independent irrespective to the above assumptions concerning the geminal-related ESPs. The minimum of the resonance energy which is the only orientation dependent contribution is achieved for some values of vectors [(e)\vec]_RmLm (in the symmetric hydride case they are directed to the vortices of tetrahedron). 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:

( ® e   RmLm  )¢ = ® e   RmLm  +d ® j   m  × ® e   RmLm  + 1 2 æ è d ® j   m  Äd ® j   m  - Ád ® j   2 m  ö ø ® e   RmLm  = = ® e   RmLm  +d ® j   m  × ® e   RmLm  + 1 2 æ è d ® j   m  (d ® j   m  , ® e   RmLm  )-d ® j   2 m  ® e   RmLm  ö ø .

(0)

These vectors must be inserted in Eq. (24) and the elasticity constant can be obtained by extracting the second order contribution in vectors d[(j)\vec] _m. In the case of hydride:

d[(j)\vec] m[(j)\vec] m¢(2)E = -2dmm¢PmrlbzsCH× × ì í î ( ® e   RmLm  ,d ® j   m  )(d ® j   m  , ® v   C m  )-d ® j   2 m  ( ® v   C m  , ® e   RmLm  ) ü ý þ .

(0)

This equation is quite remarkable since it shows that in the FO picture there are no contributions to the bending which can be attributed to interbond interaction. The bending force field is produced by energies of separate chemical bonds.

Typically in the MM framework the increment from the bending is considered as a quadratic function of valence angles. The formula for bending Eq. ( 29) can be re-written in this form. Variation of the valence angle q_mm^¢ with m < m^¢ results in rotations of the involved vectors around the axis orthogonal to the both coordination tetrahedron vectors:

d ® j   m  = - dqmm¢ 2 ® e   RmLm  × ® e   Rm¢Lm¢  ê ê ® e   RmLm  × ® e   Rm¢Lm¢  ê ê ; d ® j   m¢  = -d ® j   m  .

(0)

After substitution of this expression to the second order expansion Eq. ( 29) and after significant simplifications based on vector algebra we obtain that the bending force field constant can be written as:

kHCH = bzsCH{Pmrl( ® v   C m  , ® e   RmLm  )+Pm¢rl( ® v   C m¢  , ® e   Rm¢Lm¢  )},

(0)

which is a sum of two separate single bond contributions.

According to Eq. (15) some modification of the bending expression Eq. (31) is expected due to adjustment of hybridization tetrahedron to the geometry (coordination polyhedron) deformation. In the FO picture, which we are now considering, this modification can be only due to overall quasirotation of the hybridization tetrahedron. All the angular deformations produced by the four rotation vectors d[(j)\vec] _m form an 8-dimensional space (two components for each vector since the component collinear to the polyhedron vector can be set to zero). At the same the space of quasirotations is 3-dimensional and thus cannot accomodate all possible angular deformations. From algebraic point of view the linear mapping Eq. (14) in this case:

d ® w   l  = Ad ® j   = -( Ñ[(w)\vec] l2E) -1 å m  ( Ñ[(w)\vec] lÑ[(j)\vec] mE) d ® j   m

(0)

maps an 8-dimensional space to a 3-dimensional one, that means after using:

dimkerA+dimimA = 8,

(0)

that the transformation A has at least a five-dimensional kernel i.e. a five-dimensional subspace of the rotations space spanned by { d[(j)\vec] _m| m = 1¸4} , which results in a zero rotation of the system of the HOs. This kernel can be determined as a subspace orthogonal to the image of the operator A. The structure of the image subspace can be easily figured out since it is obvious that the rotations of a molecule as a whole

F = å m  d ® j   m

(0)

lead to the equivalent quasirotations of the set of HOs. This result can be directly applied to finding the correction to Eq. (31). As one can see, for whatever values of dq_mm^¢ the resulting deformations are all orthogonal to the overall rotation of the system as a whole (correspond to a kernel of operator A Eq. (32)), since

d ® j   m¢  +d ® j   m  = 0,

(0)

i.e. any deformation of valence angles cannot result in quasirotation of HOs and Eq. (31) is not corrected due to adjustment of HOs in the FO picture and thus the only source of the bending energy is the angular dependence of the resonance integrals b_rmlm^CH for the bonds involved in the bending.

Adjustment of form of hybridization tetrahedron

Let us now consider the picture which allows to adjust both orientation and form of hybridization tetrahedron. It can be termed as TO (tuned orbitals) in accordance with notations of Ref. []. To this end we need to systematically estimate the second derivatives of the energy with respect to variables q and x. Generally, the second order correction to the energy due to variation of HOs can be written as:

dww(2)E = -4 å m  Pmrldww(2)brmlmCH, where dww(2)brmlmCH = bssCHd(2)smC+bzsCH( ® e   RmLm  ,d(2) ® v   C m  )

(0)

This expression significantly simplifies for a symmetric case, where all s_m^C are the same ( [1/2]), the directions of all HOs coincide with those of the bonds, i.e., [(v)\vec]_m^C = [(Ö3)/2][(e)\vec]_RmLm , and the bond orders for all bonds and the resonance integrals in the corresponding DCFs coincide. Taking these relations into account we show [] that the second order energy expansion with respect to small quasi- and pseudorotation angles for the symmetric hydride (methane) takes the form:

dww(2)E = 4Prl é ê ë 2 Ö3 bzsCHd ® w   2 l  + æ ç è bssCH+ 1 Ö3 bzsCH ö ÷ ø d ® w   2 b  ù ú û ,

(0)

which is diagonal with two eigenvalues which are triply degenerate.

To write them down the necessary mixed second derivatives we consider the contribution to the energy of the second order with respect to the angles and the small variations of the resonance integrals:

dbw(2)brmlmCH = -dbssCHm(d ® w   b  , ® v   C m  )+dbzsCHmsmC(d ® w   b  , ® e   RmLm  ) -dbzsCHm(d ® w   l  , ® e   RmLm  × ® v   C m  ).

(0)

Now we consider the deformation of the hybridization tetrahedron due to angular distortions of molecular geometry. The required mixed second derivative can be easily written from Eqs. (28) and (38):

Ñ[(w)\vec] bÑ[(j)\vec] mE = 4PmrlbzsCHsmC( ® e   RmLm  ×) Ñ[(w)\vec] lÑ[(j)\vec] mE = 4PmrlbzsCH( ® e   RmLm  Ä ® v   C m  -( ® v   C m  , ® e   RmLm  )Á).

(0)

The relation Eq. (14) in this case can be written as

d ® w   = -( Ñ[(w)\vec] 2E)-1 å m  ( Ñ[(w)\vec] Ñ[(j)\vec]mE) d ® j   m  ,

(0)

where a 6×6 matrix of the second derivatives with respect to angles coming from Eq. (37) is meant by Ñ_[(w)\vec] ²E. This linear operator maps 8-dimensional space of angular molecular deformation to a 6-dimensional space of angular deformations of the HOs. Thus there exists at least a 2-dimensional subspace of angular deformations of a tetrahedron which affect neither orientation nor the form of the HOs. These deformations can be called ''hybridization incompatible''. To single them out we notice that as previously the 3-dimensional subspace of overall rotations maps to the 3-dimensional subspace of pure HOs rotations. The Ñ_[(w)\vec] ²E matrix is nondegenerate as well as its reverse matrix. Therefore, the mixed second derivative matrix has a kernel mapping some of the coordination tetrahedron deformations to zero. It means that the 5-dimensional subspace formed by independent distortions of valence angles form a 2-dimensional kernel which maps to the zero deformation of the HOs and a 3-dimensional image which maps to the 3-dimensional space of the HO deformations d[(w)\vec] _b. We perform direct calculation of these subspaces for the symmetric case with the coordination tetrahedron vectors pointing to the octants. It can be obtained the following picture of the tetrahedron deformation: if a valence angle increases by certain amount dq_1m then the opposite (spiro) angle decreases by the same magnitude. Clearly in a tetrahedron one can chose three pairs of spiro-angles for m = 2¸4. On the other hand it can be shown that the symmetric combinations of the spiro angles (when they increase or decrease simultaneously) all fall into the kernel of the å\limits_m( Ñ_{[(w)\vec] b}Ñ_{[(j)\vec] m}E) operator thus forming the 2-dimensional subspace of the hybridization incompatible deformations of the tetrahedron.

In the case of the hybridization compatible rotations of the vectors the result of action of the operator reads:

å m  Ñ[(w)\vec] bÑ[(j)\vec]mE ì í î d ® j   m  ü ý þ = 4PrlbzsRL   æ  ú Ö 2 3   æ è dc12 ® k   +dc13 ® j   +dc14 ® i   ö ø

(0)

and, therefore, the reaction of the form of hybridization tetrahedron on the angular distortions of molecular geometry can be calculated as:

d ® w   b  = - bzsCH Ö2(Ö3bssCH+bzsCH) æ è dc12 ® k   +dc13 ® j   +dc14 ® i   ö ø ,

(0)

provided the parameters dc_1m describe the hybridization compatible deformations of the coordination tetrahedron. The above expression has a remarkable feature: the variation of the hybridization angles falls back as compared to the variation of the corresponding valence angles since the coefficient relating the variation of geometry angle to the hybridization angle is less than unity.

The same considerations also apply to the bond stretching where variation of the resonance parameters can be presented as

dbmnCHm = gmnCHdrCHm.

(0)

It leads to the mixed second derivatives of the form:

¶2brmlmCHm ¶ ® w   b  ¶rCHm = -gssCH ® v   C m  +gzsCHsmC ® e   RmLm  , ¶2brmlmCHm ¶ ® w   l  ¶rCHm = -gzsCH ® e   RmLm  × ® v   C m  .

(0)

Therefore, in the symmetric case additional quasirotation is absent since the vectors [(v)\vec]_m^C and [(e)\vec]_RmLm are collinear. The responce of the form of the hybridization tetrahedron can however be written as:

d ® w   b  = - Ö3 2 · gssCH ® v   C m  -gzsCHsmC ® e   RmLm  Ö3bssCH+bzsCH drCHm.

(0)

It is well known that solution of reverse problem in molecular vibrational spectroscopy leads to off-diagonal terms in the potential energy function [] while they are usually absent in the MM schemes. The proposed scheme allows to solve the question of the off-diagonal terms on purely theoretical basis. So, coupling of stretching for two incident C-H bonds in the methane molecule can be written as:

¶2E ¶rCHm¶rCHm¢ = 1 2Ö3 Prl (Ö3gssCH-gzsCH)2 Ö3bssCH+bzsCH

(0)

and is positive. The same moves result in off-diagonal terms connecting bond stretching with angular distortions at the same atom.

Results and Discussion

In this paper we starting from the semiempirical APSLG description of molecular electronic structure produced explicit formulae for different force fields of the sp³ hybridized carbon atom. Here we discuss the limits of their applicability. First of all we should mention that some numerical analysis of schemes either fixing or tuning different ESPs is performed in Ref. []. In this paper the construction is performed in the linear response approximation. So we need estimates of numerical reliability of this model. It can be shown that even relatively large elongation of C-H bond of 0.10 Å leads to very small changes of geminal amplitudes not exceeding 0.003. The difference between estimates of d[(w)\vec] _b for this distortion in the FA setting numerically and by formula Eq. (45) is less than 0.2%. In the case of methane distortions conserving the S₄ symmetry the HOs remain unchanged even for angular deformations corresponding to valence angle 60^°. These deformations are hybridization incompatible and this conclusion derived in the linear approximation remains true even for the above large distortion. The next process studied is the totally hybridization compatible deformation of the coordination tetrahedron. In this case even for large distortion dq₁₂ = 6^° the difference between numerical and theoretical (by Eq. (42)) estimates of d[(w)\vec] _b is less than 0.03%. Therefore, the use of linear response formulae is valid.

The formulae given above allow to obtain the parameters of the MM force fields. We give here the numerical estimates. The parameter r⁰ (''equilibrium'' C-H bond length) in the FA picture equals to 1.069 Å . The symmetric TA picture Eq. (8) gives the value 1.078 Å , which is closer to the C-H bond length in the methane molecule (1.094 Å ). At the same time both these values seem quite reasonable. In the case of second derivative of the bond energy function we cannot expect a good coincidence with the same parameters adopted in the MM schemes for two reasons. First of all the MM elasticity constants have some implicit contribution of other factors. For example, the structure-related elasticity constant values implicitly refer to a certain range of nonvalence interaction energies. On the other hand we must admit that the elasticty constant obtained in the present MM derivation strongly depends on the form of specific core-core repulsion D_CH adopted in the MINDO/3 method. Within the original method this repulsion term is parameterized only to reproduce the equilibrium bond lengths (first derivatives of the energy) but not the second derivatives of the energy (MINDO/3 method is not parameterized to reproduce the IR frequencies). Indeed, we obtain the value 8.30 mdin/Å in the framework of the FA picture and 7.77 mdin/Å in the framework of the symmetric TA picture. The standard MM schemes usually use this elasticity constant in the range from 4.5 to 4.7 mdin/Å MM2,EAS,CFF3,Boyd. The constant obtained from the solution of reverse spectral problem in the assumption of harmonic potential is equal to 5.31 mdin/Å []. At the same time Ref. [] gives the value of 7.9 mdin/Å for this constant.

The estimates of other constants based on analysis of the resonance contribution to the energy seem to be more reliable because the bending force field derived does not depend on any derivatives of molecular integrals on geometric parameters but is determined by coarse features of electronic structure like angle between HO and direction of chemical bond. Eq. (31) allows to estimate the bending elasticity constant for the H-C-H angle. It equals to 0.509 mdin/deg. and is quite close to the constants used in the MM force fields (for example, 0.549 mdin/deg. in Ref. [] and 0.508 mdin/deg. in Ref. []) and to the value 0.493 mdin/deg. obtained by solving the reverse spectral problem for methane Gribovetal. It shows valdity of the stereochemistry considerations based on the resonance contribution to the energy and confirms our conclusions about the main source of the bending energy. The explicit formula Eq. ( 46) allows to estimate the off-diagonal constant h coupling the stretching of two incident C-H bonds. The numerical value of this constant obtained from our consideration is 0.120 mdin/Å . This value is approximately three times larger than that estimated in Ref. Gribovetal in the harmonic approximation. It is not surprising since for such a small constant its value strongly depends on the force field used.

Conclusion

In this paper we proposed a way of analytical construction of the principal MM force fields - those of stretching and bending. To this end we performed a formal analysis of the energy expression written within the semiempirical implementation of the APSLG method and thoroughly elaborated different types of approximations to its electronic structure parameters either matrix elements of electron density or hybridization parameters. Both types of parameters can be either fixed or tuned. Here we developed a linear response version of the theory with fixed geminal ESPs (FA) and applied it to very characteristic example - methane molecule. We constructed bond energy functions and using the parametrization of the hybridization manifold by a pair of quaternions gave explicit formulae for different force fields including off-diagonal ones which depend on the details of hybridization response to the small geometry variations. The general theoretical consideration is confirmed by numerical analysis of force field parameters.

Acknowledgements

This work has been performed with partial 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 Tops oe A/S is acknowledged as well. ALT gratefully acknowledges valuable discussions with Prof. A.A. Levin, Dr. I.A. Misurkin, and Dr. I.V. Pletnev.

References

[]

Burkert, U.; Allinger, N. Molecular Mechanics; ACS: Washington, DC, 1982.

[]

Tokmachev, A. M.; Tchougréeff, A. L. Int J Quantum Chem 2001, 84, 39.

[]

Gillespie, R. J.; Hargittai, I. The VSEPR Model of Molecular Geometry; Allyn & Bacon, 1991.

[]

Cohen, N.; Benson, S. W. Chem Rev 1993, 93, 2419.

[]

Daudey, J.-P.; Malrieu, J.-P.; Rojas, O. In Localization and Delocalization in Quantum Chemistry, Vol. 1, Atoms and Molecules in the Ground State; Chalvet, O.; Daudel, R.; Diner, S.; Malrieu, J.-P., Eds.; D. Reidel Publ.: Dordrecht et al., 1975.

[]

Diner, S.; Malrieu, J.-P.; Claverie, P. Theoret Chim Acta 1969, 13, 1.

[]

Tokmachev, A. M.; Tchougréeff, A. L. Int J Quantum Chem 2002, 88, 403.

[]

Tokmachev, A. M.; Tchougréeff, A. L. J Comp Chem 2001, 22, 752.

[]

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

[]

Tchougréeff, A. L.; Tokmachev A. M.; in preparation.

[]

Hoffmann, R. J Chem Phys 1963, 39, 1397.

[]

Coulson, C. A. Valence; Oxford University Press: Oxford, 1961.

[]

Elliot, J. P.; Dawber, P. G. Symmetry in Physics, Vols. 1, 2; Macmillan Press Ltd: London, 1979.

[]

Volkenstein, M. V.; Gribov, L. A.; Elyashevich, M. A.; Stepanov, B. I. Molecular Virbations; Nauka: Moscow, 1972 [in Russian].

[]

Allinger, N. L. J Am Chem Soc 1977, 99, 8127.

[]

Engler, E. M.; Andose, J. D.; Schleyer, P. v. R. J Am Chem Soc 1973, 95, 8005.

[]

Ermer, O.; Lifson, S. J Am Chem Soc 1973, 95, 4121.

[]

Chang, S.; McNally, D.; Shary-Tehrany, S.; Hickey, S. M. J.; Boyd, R. H. J Am Chem Soc 1970, 92, 3109.

[]

Pearson, R. G. Symmetry Rules for Chemical Reactions; Wiley-Interscience: New York et al., 1976.