Content-type: text/HTML

A.M. Tokmacheva1  and A.L. Tchougréeffb
a Institute of Inorganic Chemistry, RWTH Aachen
Prof.-Pirlet-Str. 1, 52074 Aachen, Germany
b Karpov Institute of Physical Chemistry
Vorontsovo pole 10, 105064 Moscow, Russia

Transferability of Parameters of Strictly Local Geminals' Wave Function and Possibility of Sequential Derivation of Molecular Mechanics

Transferability of Parameters of Strictly Local Geminals' Wave Function and Possibility of Sequential Derivation of Molecular Mechanics

Dedicated to Prof. N.F. Stepanov
on the occasion of his 70-th birthday.

Abstract

The problem of substantiation of molecular mechanics (MM) remains actual due to growing popularity of hybrid quantum/classical (QM/MM) schemes. Recently proposed deductive molecular mechanics (DMM) seems to be a natural tool to derive mechanistic models of molecular energy (classical force fields) from a suitable quantum mechanical (QM) description of molecular structure. It is based on an assumption that the trial wave function underlying the MM description is one of the antisymmetrized product of strictly local geminals (SLG). A proof of transferability of electronic structure parameters (ESPs) in this approximation is an essential component of a logical framework for the transition from the QM to an MM description since it allows constructing expressions for potential energy surfaces by proper consideration of the response of the ESPs to the variations of geometry parameters. In the present paper the ESPs defining density matrix elements and basis one-electron states (hybrid orbitals - HOs) in the SLG approximation are formally considered. The transferability of the density mareix elements with respect to the parameters of molecular electronic structure and the linear response relations for the HOs are proven to take place under very nonrestrictive conditions. Special attention is paid to numerical estimates of the ESPs' features giving an ëxperimental" support to approximate expressions for the molecular energy.

Introduction

Molecular mechanics (MM) [,] is currently a versatile and very popular tool in laboratory and industrial practice. It remains the most practical way of analyzing potential energy surfaces (PES) for large molecules even though linear scaling O(N) methods of quantum chemistry are becoming available [,]. The nature of the MM as a combination of totally empirical classical force fields allows to realize its main drawbacks (for example, inapplicability to highly correlated systems) and advantages (fast evaluation of the total energy and its gradient and high accuracy of molecular geometries obtained) and to characterize them as the most preferable according to the "quality/cost" criterion but with the field of application limited to certain combination of properties, processes, and classes of molecules.

Despite its long history and a wealth of successful applications the MM approach remains not substantiated theoretically. It is clear that the molecular PES can be expanded up to the second order in nuclear displacements near the equilibrium geometry and this fact is often considered as a general substantiation of the possibility to use the MM-type expansions for the total energy []. However, it is very questionable since (i) the transferability of the elements of the dynamic (second derivatives) matrix between molecules is not proven, and (ii) the geometry variables (bond lengths and valence angles) used in the definitions of the force fields do not diagonalize the dynamic matrix even approximately, whereas the standard MM PESs are by construction diagonal in the bond-length/valence-angle representation.

Substantiation of MM is not merely an academic question. The last years have demonstrated a growing interest [,,,] to hybrid quantum mechanical/molecular mechanical (QM/MM) schemes where different parts of the entire molecule are treated by different methods. The relatively small part of the molecule (reactive center) is treated by an adequate quantum mechanical (QM) method while the rest (inert environment) is described by classical force fields. Thus constructed QM/MM methods combine the universality of the QM description and the low cost of the MM that makes them an optimal tool for analysis of the structure of large molecules including those of biological importance. Though the QM/MM methods shift the limits for numerical applications they require more insight into the electronic structure underlying the MM schemes since it is not clear how to construct the junction between quantically and classically treated subsystems. Though many ad hoc recipes for the junction construction are proposed in the literature (see, for example, reviews [,,]), the problem has not been yet solved. Meanwhile, the absence of sequential derivation for the intersubsystem junction leads to numerous artefacts in QM/MM schemes reviewed in Ref. []. The general formal solution of this problem proposed in Refs. [,] is based on the assumption that a generic MM scheme can be obtained from a suitable QM method by a series of approximations. Thus the prospect of development of effective QM/MM schemes constitutes the practical aspect of the general task of the MM substantiation. An important prerequisite for the announced derivation of MM from QM is the choice of the quantum chemical method underlying the MM description. By taking a proper form of the trial wave function having some common features (or resemblance) with the MM schemes we can significantly simplify the task. The representation of the total energy of the covalently bound molecule in the current MM schemes is tightly connected to the concept of two-center bonds. The energy is taken as a sum of intrabond, bending, torsion and non-bonding contributions. The specific forms of force fields depend on the implementation and quite a lot of schemes with systems of parameters fitted to reproduce various characteristics of different classes of molecules have been proposed [,,,,,,]. In the course of the evolvement of the MM approach increasingly sophisticated contributions are added to the original simplistic picture [] which allow to extend the MM approach to more and more complex and not transparently tractable classes of the molecules [,,]. Thus an acceptable substantiation of the MM should not be reduced to deriving or validating any specific MM scheme. By contrast, a generic mechanistic picture must be obtained on the basis of the adequate QM description.

The standard QM methods in that or another manner based on the delocalized picture of molecular electronic structure provided by the SCF approximation are not quite suitable for the purpose of the MM substantiation. First, the general SCF energy expression contains contributions of the same form for all pairs of atoms and separation of the energy into bonding and non-bonding contributions is not built in to the SCF methods. Second, even applying different localization techniques intended to reconcile the SCF picture of the electronic structure with the chemist's view relying upon the bonds concept leads not to very much progress due to "tails" of the localized orbitals which are difficult for handling and absolutely non-transferable (it is also impossible to prove theoretically the transferability of the MO LCAO coefficients or, even, of the SCF density matrix elements).

The above discussion allows to specify some criteria for selecting a QM method suitable to serve as a starting point for the MM derivation: discrimination and different treatment of bonding and non-bonding contributions to the energy, strict locality of one-electron states, variational determination of the electronic structure parameters (ESPs). Failure to fulfil any of these criteria leads to a necessity to take whichever of the mentioned features as assumptions (as for example, in Ref. [], where the analysis was based on the PCILO method []).

The unsuitabilty of the SCF approximation as a basis for the additivity concept in general (and thus for the MM in particular) was clear even on very early stage of the theory development. Yet then (see []) it had been proposed to use a geminal based description for substantiating the additive methods and the bond energy concept. However, the one-electron carrier spaces to be used for geminal consctuction had not been sufficiently specified and the entire geminal-based scheme had not been explored to a due extent by adequate numerical studies. Only recently we proposed [,] a semiempirical method which satisfies the suitability criteria formulated above. It is based on the trial wave function in the form of antisymmetrized product of strictly local geminals (SLG). Variational determination of strictly local one-electron states and of the density matrix elements provides the flexibility of the wave function necessary to describe the electronic structure and properties of organic compounds. Numerical estimates based on the MINDO/3 parameterization of the semiempirical Hamiltonian [,] and three (MNDO, AM1, and PM3) parameterizations of the NDDO family [] have shown that the method proposed supercedes the SCF one in description of the heats of formation, molecular geometries, and ionization potentials of organic molecules. Incidentally, it provides a good basis for description of molecular electronic structure in chemical terms. The method is well suited for description of large molecules since it belongs to a class of O(N) methods and can be used for constructing hybrid schemes []. Also the approximate transferability of the ESPs of the SLG approach was numerically demonstrated in Ref. [], where it was shown that the energy calculated with some ESPs "fixed" and characteristic for a particular class of atoms and for bonds is close to that obtained by direct minimization of the SLG energy. More refined treatment was proposed in Refs. [,], where linear response relations for the form of one-electron states were used to derive explicit expressions for the angular dependence of the energy in the case of sp3 hybridized carbon and nitrogen atoms. The numerical estimates of the MM force field parameters obtained in [] on the basis of the analytical expressions for the constants of force fields are close to those accepted in the MM. Special attention was paid to off-diagonal force fields and to possible sources of the angular dependence of the energy [] as well as those of piramidalization potential in nitrogen-containing compounds []. The development of these ideas has led to the formulation of the "deductive molecular mechanics" [], which represents atoms by their hybridization tetrahedra (see below) of different shapes adjusted according to geometry and valence state variations rather than by harmonically interacting point masses ("balls-and-springs") used in standard MM. All these results are obtained within the assumptions of (i) the perfect transferability of ESPs characterizing chemical bonds and (ii) of the validity of linear response relations for hybridization with respect to geometry variations. In the present paper we address the problem of transferability of the density related ESPs, explore the precision and the validity limits for the linear response formulae for the shapes of hybridization tetrahedra and consider the possibility of recovering the standard MM description from the deductive molecular mechanics.

The paper is organized as follows: in Section we briefly review the main features of the SLG method (with the MINDO/3 semiempirical parameterization) relevant to constructing the MM description. In Section the ESPs characterizing the density distribution in chemical bonds and lone pairs as they appear in the SLG method are analyzed and the features assuring their transferability are singled out. In Section we briefly discuss the structure of the hybridization manifold and give expressions for the variations of the hybridization related ESPs in response to geometry changes. Results of numerical experiments are given in these Sections to support our theoretical derivation of transferability and of linear response relations. In Section we discuss different possible approximate descriptions of the ESPs and energy leading to various versions of the deductive MM theory. In Section we consider the possibility of deriving the standard MM picture from the DMM by projecting out excessive (from the MM point of view) variables of the latter and give theoretical expressions for the force fields' parameters of the standard MM. Finally we discuss the relation between the transferability of the density related ESPs and that of the MM force fields.

SLG description of molecular electronic structure

Constructing the SLG trial wave function according to [,] requires the following moves. First, the one-electron basis of the strictly local hybrid orbitals (HOs) must be contructed. These orbitals are obtained by an orthogonal transformation of the s and p AOs for each "heavy" (non-hydrogen) atom. These transformations are represented by 4×4 orthogonal matrices hA Î O(4) for each heavy atom A. All the HOs are assigned either to respective two-electron chemical bonds or to electron lone pairs. Each chemical bond refers to two such HOs - | r ñ and | l ñ (right- and left-end ones, respectively). Each lone pair is formed by one HO only (a right one for the sake of definiteness).

Chemical bonds and lone pairs are described by singlet two-electron functions - geminals [] taken in the form originally proposed by Weinbaum []. With use of the second quantization notation they are written as:
gm+ = umrma+rmb++vmlma+lmb++wm(rma+lmb++lma+rmb+)
(0)
for chemical bonds and:
gm+ = rma+rmb+
(0)
for lone pairs. The normalization condition imposed on the amplitudes reads:
um2+vm2+2wm2 = 1.
(0)
The geminals thus defined in the carrier space spanned by so constructed HOs are termed to be strictly local geminals (SLG).

The wave function of electrons in the molecule is then taken as the antisymmetrized product of the geminals given by Eqs. (1), (2):
| F ñ =
Õ
m 
gm+| 0 ñ .
(0)

The SLG energy is a function of the intrabond matrix elements of spinless one- and two-electron density matrices:
Pmtt¢ = á 0| gmtms+tms¢gm+| 0 ñGmtt¢ = á 0| gmtmb+tma¢+tma¢tmbgm+| 0 ñ ,
Pmrr = um2+wm2Pmll = vm2+wm2Pmrl = Pmlr = (um+vm)wm,
Gmrr = um2Gmll = vm2Gmrl = Gmlr = wm2,
(0)
where t and t¢ ( = ±1) correspond to the right (+1) and to the left (-1) ends of the bond and to the Fermi operators r and l, respectively. Certain energy contributions further reduce to Coulomb interaction of effective atomic charges residing on the atoms. Applying the Mulliken scheme for the charges to the SLG wave function yields them in the form:
QA = 2
å
tm Î A 
Pmtt-ZA.
(0)

If the Hamiltonian of the MINDO/3 form [] in the HO basis is used with the SLG trial wave function the total energy can be written in a form somewhat close to the MM energy with interactions between bonded and non-bonded atoms treated in different ways and closely relates it to that given in Ref. [] in the context of analysis of a variety of additive schemes of molecular energy:
E =
å
A 
EA+
å
A < B 
EAB,
EA =
å
tm Î A 
[2UmtPmtt+(tmtm|tmtm)TmGmtt]+2
å
[(tktm¢ Î A) || (k < m)] 
gtktm¢TkPkttPmt¢t¢,
ERmLmbond = 2gRmLm[Gmrl-2PmrrPmll]- 4brmlmRmLmPmrl,
EABnonbond = QAQBgAB + ZAZBD AB,
(0)
where Tm refers to the atom bearing the HO tm; DAB is MINDO/3 specific core-core repulsion; the reduced Coulomb integral is:
gtktm¢Tk = 2(tktk|tm¢tm¢)Tk-(tktm¢|tm¢tk)Tk.
(0)
and the resonance (electron hopping) one is:
btm1tm2¢AB =
å
i Î A 

å
j Î B 
hm1iAhm2jBbijAB.
(0)
In this energy expression it is assumed that the nonbonding contribution appears for all pairs of atoms whereas the bonding only for those with an actual bond.

Molecular integrals entering the above expressions depend on molecular geometry and on the orthogonal matrices hA for all non-hydrogen atoms A. The expressions for the matrix elements in the HO basis are given in Ref. [].

The energy expression Eq. (7) depends on the amplitudes of bond geminals through the values of Pmtt¢ and Gmtt¢, and on the form of HOs through the molecular integrals. The amplitudes of two ionic (um, vm) and one covalent (Heitler-London type, Ö2wm) configurations in Eq. (1) for each geminal are determined with use of the variational principle as well as the matrices hA of transformation of the AO basis to the HO one.

This comprises the essence of the semiempirical SLG method.

Electronic structure parameters in the SLG picture

In the previous Section we reviewed briefly the semiempirical implementation of the SLG method for analysis of electronic structure and expressed the total molecular energy in the form Eq. (7) which allows the representation of the molecular PES as a sum of local increments. These increments depend on the ESPs of two classes (i) those defining the hybridization of atomic basis sets and (ii) the intrabond density matrix elements. In this Section we concentrate upon the proof of transferability of the electron density matrix elements, related to the geminal amplitudes, and on the structure of the hybridization manifold as it appears in the SLG approximation.

Density matrix elements and their transferability

As it is reported in Section 1 the energy in the SLG approximation is a function of one- and two-electron density matrices. Their matrix elements are in turn expressed through the geminal amplitudes, appearing while diagonalizing the effective bond Hamiltonians. Thus any analysis of the properties of the density ESPs starts from description of the latter.

Effective bond Hamiltonians

Within the original SLG approach [,] the geminals are characterized by the amplitudes (see Eq. (1) um, vm, and Ö2wm = zm, which simplifies the normalization condition Eq. (3) for the amplitudes to: um2+vm2+zm2 = 1. To determine them the effective Hamiltonians for each bond geminal are constructed. The optimal values of these amplitudes are the solutions of the eigenvector problem (see also []):
æ
ç
ç
ç
è
Am
Dm
0
Dm
Cm
Dm
0
Dm
Bm
ö
÷
÷
÷
ø
æ
ç
ç
ç
è
um
zm
vm
ö
÷
÷
÷
ø
= em æ
ç
ç
ç
è
um
zm
vm
ö
÷
÷
÷
ø
,
(0)
corresponding to its lowest eigenvalue.

The matrix elements of the effective bond Hamiltonians are defined as (with the MINDO/3 Hamiltonian):
Am = 2Umr+(rmrm|rmrm)Rm-4gRmLmPmll+
+2
å
B ¹ Rm 
gRmBQB+2
å
[(tm1 Î Rm) || (m1 ¹ m)] 
grmtm1RmPm1tt;
Bm = 2Uml+(lmlm|lmlm)Lm-4gRmLmPmrr+
+2
å
B ¹ Lm 
gLmBQB+2
å
[(tm1 Î Lm) || (m1 ¹ m)] 
glmtm1LmPm1tt;
Cm = 1
2
(Am+Bm)-DgmDm = -Ö2brmlmRmLm,
(0)
where
Dgm = gm-gRmLmgm = 1
2

å
t Î {r,l} 
(tmtm|tmtm)Tm.
(0)

The calculations of Refs. [,] performed on organic compounds of different classes (alkanes, alcohols, amines etc.) have demonstrated a remarkable stability of all the geminal related ESPs. The values of the polarity Prrm - Pllm do not exceed 0.07 by absolute value for the compounds containing carbon, nitrogen, and hydrogen atoms (for the situation with oxygen and fluorine see below). Also the ionicity u2m + v2m for a rich variety of bonds has a stable value about 0.4. The bond orders 2Prlm all acquire values between 0.92 and 1.0. These features though not completely unexpected, since the transferability of the parameters of the single bonds in organic compounds is well known experimentally, require a theoretical explanation.

Pseudospin representation and the perturbative estimates of the bond-geminal ESPs

In order to provide the required explanation we notice that the effective Hamiltonians for the bond geminals can be represented as a sum of the unperturbed part which when diagonalized yields an invariant, i.e. exactly transferable, values of the ESPs and of a perturbation responsible for specificity of different chemical compositions and environments.

Pseudospin operator of the bond geminal  

Let us introduce a pseudospin [^(t)]m operator corresponding to the pseudospin value tm = 1. The matrices of its components in the basis of the configurations defining the geminal are given by:
^
t
 

zm 
= æ
ç
ç
ç
è
1
0
0
0
0
0
0
0
-1
ö
÷
÷
÷
ø
,
^
t
 

+m 
= æ
ç
ç
ç
è
0
Ö2
0
0
0
Ö2
0
0
0
ö
÷
÷
÷
ø
,
^
t
 

-m 
= æ
è
^
t
 

+m 
ö
ø
f
 
.
(0)
The configurations with á [^(t)]zm ñ = ±1 are the ionic ones with both electrons located on the same end of chemical bond (right or left, respectively). In terms of the pseudospin operator averages the quantities in Eq. (5) are expressed as follows:
Pmtt = 1
2
(1+t
^
t
 

zm 

), Pmrl = 1
2

^
t
 

+m 

Pmlr = 1
2

^
t
 

-m 

,
Gmtt = 1
2
(
^
t
 
2
zm 

+t
^
t
 

zm 

), Gmrl = Gmlr = 1
2
(1-
^
t
 
2
zm 

).
(0)
The effective bond Hamiltonians can be rewritten in terms of the pseudospin operators introduced just above. Indeed, the effective Hamiltonian Eq. (10) for each of the bond-representing geminals has the form:
Hmeff = H0meff+DIm ^
t
 
2
zm 
+DPm ^
t
 

zm 
,
H0meff = æ
ç
ç
ç
è
Cm
Dm
0
Dm
Cm
Dm
0
Dm
Cm
ö
÷
÷
÷
ø
,
DIm = 1
2
(Am+Bm)-Cm = DgmDPm = 1
2
(Am-Bm),
(0)
where the part of the effective Hamiltonian proportional to [^(t)]zm2 is responsible for reproducing the relative contribution of the ionic and covalent configurations to the bond geminal, and the contribution proportional to [^(t)]zm relates to the asymmetry (polarity) of the bond. The geminal amplitudes obtained by diagonalizing the unperturbed bond Hamiltonians H0meff and the density or [^(t)]-type ESPs thus obtained are perfectly invariant:
u0m = v0m = w0m = 1
2
æ
ç
è
z0m = 1
Ö2
ö
÷
ø
P0mtt¢ = 1
2
G0mtt¢ = 1
4
,

^
t
 

zm 


0 
= 0, 
^
t
 
2
zm 


0 
= 1
2

^
t
 

+m 


0 
=
^
t
 

-m 


0 
= 1
(0)
with respect to the sorts of atoms and to the geometry parameters though the unperturbed effective Hamiltonians H0meff themselves are, of course, not. The invariant values of the ESPs are also rather close to the exact SLG values obtained in Refs. [,].

Perturbative estimate of ESPs with respect to noncorrelated bare Hamiltonian  

One can try to estimate the optimal values of the ESPs specific for each bond and molecule perturbatively by using the linear response approximation []. According to the latter the response d á A ñ of a quantity described by the operator A to the perturbation lB of the Hamiltonian (where l is the parameter characterizing the intensity of the perturbation) has the form:
d á A ñ = l á á A;B ñ ñ ,
(0)
where the zero frequency Green's function in the energy domain is given by the relation:
á á A;B ñ ñ = 2
å
i ¹ 0 
á 0| A| i ñ (e0-ei)-1 á i| B| 0 ñ .
(0)
Inserting the amplitudes of the geminals for the ground and the two excited states of the bond and two excitation energies (which are equal to Ö2| Dm| and 2Ö2| Dm| ) for the effective Hamiltonian H0meff to the general expression for the Green's functions Eq. (18) one can immediately check that only the diagonal (with respect to A and B running over the relevant components of the pseudospin operator and its squared z-component) Green's functions are nonvanishing. Therefore, the first-order responses of the averages of the pseudospin operator components can be written as:

^
t
 

zm 

= d
^
t
 

zm 

= Bm-Am
4brmlmRmLm
,
d
^
t
 
2
zm 

= - Dgm
8brmlmRmLm
,
d
^
t
 

+m 
+ ^
t
 

-m 

= 0.
(0)
The last expression demonstrates that at least in the linear response approximation the bond order does not change, (i.e. is invariant even for the different atoms forming the bond and geometries). This result, though suggests certain stability of the bond orders, in fact is a consequence of the SCF type of the two-electron wave function implicitly assumed by the decomposition of the bond effective Hamiltonian Eq. (15). The unperturbed part thus defined yields the symmetric ground state with the total weight of the ionic configurations equal to that of the covalent one. This coincides with the result characteristic for the SCF approximate treatment of the symmetric bond. It is obvious that the formulae for the bond-order variation and the two-electron density matrix elements (and the SCF approximation itself) are not valid at larger interatomic separations (the denominators in Eq. (19) become too small) though the exact solution of the SLG problem has a correct asymptotic behaviour for single bonds even at infinite interatomic separations. This attractive feature of the SLG model is lost in the perturbative treatment based on the Hamiltonian separation Eq. (15) as it is in the SCF approximation as well.

Perturbation of the density matrix elements for correlated ground state  

In order to overcome the above failure of the perturbative estimation of the two-electron density and of the bond orders let us consider a symmetric bond. This would correspond to a different decomposition of the effective bond Hamiltonian than that of Eq. (15). We assume that the contribution to the effective bond Hamiltonian which is proportional to [^(t)]zm2 is included into the unperturbed (zero order) Hamiltonian. The problem then reduces to a 2×2 matrix diagonalization. The ESPs, as they appear from solution of this problem, are:
Gmtt¢ = 1
4
æ
ç
è
1-tt¢ 1
G(zm)
ö
÷
ø
Pmtt = 1
2
Pmrl = zm
2G(zm)
.
(0)
where
zm = 4brmlmRmLm/DgmG(zm) =
Ö
 

1+zm2
 
(0)
Small interatomic separations characteristic for the real bonds correspond to the limit zm >> 1 and the ESPs have the asymptotic behaviour as:
G0mtt¢ = 1
4
æ
ç
è
1-tt¢ 1
zm
ö
÷
ø
Pmrl = 1
2
æ
ç
è
1- 1
2zm2
ö
÷
ø
.
(0)

Now, when the total ionic contribution to the geminal is calculated exactly (variationally), the bond polarity can be estimated perturbatively in the linear response approximation, but with the correlated ground state of the symmetric effective bond Hamiltonian taken for evaluating the Green's function. It can be conveniently done with use of a dimensionless asymmetry parameter:
mm = Bm-Am
DgmG(zm)
.
(0)
Inserting the corresponding ground state to the definition of the Green's function á á [^(t)]zm;[^(t)]zm ñ ñ one obtains:

^
t
 

zm 

= mm G(zm)-1
G(zm)+1
(0)
for the polarity of the bond between atoms with the fixed hybridizations. It vanishes for infinite interatomic separation as it should be for the exact wave function. The bond ionicity and the bond order are not affected in the linear response approximation.

The expression for the bond polarity coincides with that of Eq. (24) even if the second-order perturbation correction to the wave function is used (i.e., the contribution to the bond polarity proportional to mm2 is absent). At the same time the second-order corrected bond ionicity and bond order have the following form:

^
t
 
2
zm 

=
^
t
 
2
zm 


c 
é
ê
ë
1+mm2 2G(zm)+1
2(G(zm)+1)
ù
ú
û
,

^
t
 

+m 
+ ^
t
 

-m 

=
^
t
 

+m 
+ ^
t
 

-m 


c 
é
ê
ë
1+mm2 2G(zm)+1-G2(zm)
2(G(zm)+1)2
ù
ú
û
,
(0)
where the ESPs with the subscript c correspond to the estimates by Eq. (20).

Lone pairs

Another archetypical form of two-electron group is the lone pair. As it is mentioned above the lone pair is described by a degenerate geminal containing the contribution of only one ionic configuration. For the sake of definiteness we set it to be the right-end ionic configuration of the corresponding degenerate bond (the amplitude um becomes equal to unity, see Eq. (2)). The ESPs related to the lone pair can be readily evaluated:

^
t
 

zm 

= d
^
t
 

zm 

= 1, d
^
t
 
2
zm 

= 1
2
,
Pmrl = Pmlr = 0, Gmrr = 1, Gmll = Gmrl = Gmlr = 0.
(0)
These quantities are perfectly invariant and transferable from one molecule to another and basically characterize (within the accepted approximation, of course) the qualitative difference between the atoms of different chemical elements: by the number of lone pairs they bear.

Numerical estimates of density ESPs' transferability

The above analytical results must be controlled by numerical estimates in order to get a feeling of the real sense of the "first" and ßecond" orders. Table represents the results of calculations on the ESPs á [^(t)]zm ñ , á [^(t)]zm2 ñ , and á [^(t)]+m ñ by the SLG method (Eq. (10)) and by the approximate formulae Eqs. (19), (20), (22), (24), and (25) for some characteristic bonds in small molecules. The results show that in the case of bonds with small polarity all the formulae perform very well. The most precise approximations Eqs. (24) and (25) give results which perfectly coincide with the exact (SLG-MINDO/3) ones even for very polar O-H and F-H bonds. Also estimates according to the asymptotic (zm >> 1) formulae Eq. (20) give reasonable results for the ESPs of the bonds in not too polar molecules at their equilibrium geometries. The main source of stability of the bond order values is the validity of the above limit which in its turn takes place due to the fact that the difference between one- and two-center electron-electron repulsion integrals (Dgm) at interatomic separations characteristic for chemical bonding is much smaller than the resonance interaction at the same distance. The data of Table illustrate the difference between that which may be called MM atom types. For example, the primary C-H bonds in the ethane and propane molecules have very similar ESPs at the same time somewhat differing from those for the secondary C-H bonds in the propane molecule.

Further analysis of the quantities mm allows to single out two types of factors loaded upon this parameter: those related to the bond itself (which are again hybridization dependent) and the rest describing the environment of the bond. These factors contribute additively:
mm = m0m+m1m,
(0)
where the intrinsic bond-related part is:
m0m = 1
DgmG(zm)
[2(Uml-Umr)+(lmlm|lmlm)Lm-(rmrm|rmrm)Rm+
+
å
[(tm1 Î Lm) || (m1 ¹ m)] 
glmtm1Lmnm1-
å
[(tm1 Î Rm) || (m1 ¹ m)] 
grmtm1Rmnm1],
(0)
where nm equals to 1 for an incident chemical bond and to 2 for a lone pair. This contribution is characteristic for the pair of atoms RmLm with given ratios of the s- and p-weights in the HOs |rmñ and |rmñ ascribed to the bond at hand and clearly depending on the chemical nature of these atoms through specific values of atomic parameters and numbers of lone pairs.

The contribution to the bond asymmetry coming from the environment of the bond is:
m1m = 1
DgmG(zm)
[2
å
B ¹ Lm 
gLmBQB-2
å
B ¹ Rm 
gRmBQB+4gRmLm
^
t
 

zm 

+
+2
å
[(tm1 Î Lm) || (m1 ¹ m)] 
glmtm1Lmtm1
^
t
 

zm1 

-2
å
[(tm1 Î Rm) || (m1 ¹ m)] 
grmtm1Rmtm1
^
t
 

zm1 

].
(0)
In the molecules where all bonds are weakly polar one can expect that the external part m1m is small. Also, in the molecules containing many polar bonds, the effect of randomly distributed effective atomic charges almost vanishes by this leading to small values of external contributions m1m. Thus the only situation when one could expect the environment to affect the characteristics of the otherwise transferable bond is that when the bond under consideration falls to a close vicinity with a few strongly charged atoms arranged in such a way that their fields sum up to a nonzero overall field acting along the bond. That is, clearly, one of the situations which elaborated MM parameterizations mark as a special one, requiring specific values of parameters. The separation of mm according to Eq. (27) allows to obtain transferable, environment independent approximating functions for the ESPs by substituting into Eqs. (24) and (25) the parameter m0m instead of mm. The numerical results are given in Table . It shows that this approach leads to the approximate ESPs perfectly coinciding with those obtained by the SLG method itself that demonstrates the applicability of such a scheme. The comparison of Tables and allows to single out the effects of the environment. For example, the primary C-H bonds in the ethane and propane molecules result in coinciding values of the ESPs estimated by using m0. Therefore, the small difference between the ESPs of these bonds obtained by the SLG method is totally caused by the slightly different environment, i.e. by the m1 values which are equal to 0.003 and 0.006, respectively. The small magnitude of the deviations between the SLG and approximately estimated ESPs can be rationalized by the smallness of their effect on the total energy of the molecule. For example, even in the case of rather polar water molecule the approximation of the bond ESPs by their values from Table leads to increase of the total energy by only 0.014 kcal/mol as compared to the exact SLG calculation.

Mathematical description for hybridization in SLG picture

In the framework of the SLG scheme the structure of one-electron basis states is defined by orthogonal transformations of AOs for each atom with an sp-valence shell. The energy expression Eq. (7) is the function of the parameters defining these transformations. The 4×4 O(4) matrix hA of transformation from the AO to the HO basis set on the atom A depends on six angular variables. Three of them (pseudorotation angles [(w)\vec] b = (wsx,wsy,wsz) with subscripts indicating pairs of basis AOs mixed by the corresponding 2×2 Jacobi rotations) define the structure of the HOs (s-/p-mixing and relative directions of the HOs) while other three (quasirotation angles [(w)\vec] l = (wyz,-wxz,wyz)) define the SO(3) matrix performing rotation of the set of four HOs as a whole (prefix quasi refers to the fact that no physical body rotates under its action, only the system of HO's). Generally, the transformation of orbitals caused by a pseudorotation forms a set of HOs which is known as hybridization pattern (like sp3, sp2 etc.) which is more or less stable, while the set of quasirotation angles is totally non-transferable, depends on the relative placement of bonded atoms and, obviously, is governed by the resonance contribution to the energy since only the latter depends on the directions of the HOs.

The mathematical description of hybridization is based on employing the algebraic group structure of the hybrids' manifold. Due to the latter any small variation of HOs in a vicinity of a given set of HOs represented by a 4×4 orthogonal matrix h can be expressed with use of the SO(4) matrix H close to the unity matrix:
H = I+d(1)H+d(2)H,
h¢ = Hh » h+d(1)h+d(2)h,
d(1)h = d(1)Hhd(2)h = d(2)Hh.
(0)
The general form of matrix H is analyzed in [] where its expansion upto the second order in the vicinity of the unity matrix is given. Also simple expressions for the first-order variation of the structure of the HOs due to small quasi- and pseudorotations d[(w)\vec] l and d[(w)\vec] b applied to the set of HOs at a given atom are derived:
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 and we use a quaternion (s,[(v)\vec]) representation (see [] and below) of the HOs where the coefficient s at the s-AO in the HO LCAO expansion is a scalar and the coefficients at the p-AOs form the vector part [(v)\vec] of the quaternion. Since we know the dependence of the molecular integrals (and, therefore, energy) on the HO LCAO coefficients we find according to Eq. (31) their dependence on the ESPs [(w)\vec] b and [(w)\vec] l defining hybridization.

Deductive molecular mechanics and possibility of derivation of standard MM

Basics of deductive molecular mechanics

The analysis of the properties of the ESPs pertinent to the SLG approximation performed in Sections 0.1 and 0.1 allows to rewrite the energy Eq. (7) as follows:
E =
å
m 
é
ê
ë
æ
ç
è
2Um+ Dgm
2
-2brmlmRmLm ö
÷
ø
+ 1
2

å
k < m 

å
tt¢ Î {r,l} 
dTkTm¢gtktm¢Tk ù
ú
û
+
(a)
+
å
A < B 
ZAZBD AB+
(b)
+
å
m 

å
t Î {r,l} 
é
ê
ë
t
^
t
 

zm 

(Umt+ 1
2
(tmtm|tmtm)Tm) +
(c)
+ 1
2

å
k < m 

å
tt¢ Î {r,l} 
dTkTm¢gtktm¢Tk æ
è
t
^
t
 

zk 

+t¢
^
t
 

zm 

ö
ø
ù
ú
û
+
(d)
+
å
m 
é
ê
ë

å
t Î {r,l} 
1
2
(tmtm|tmtm)Tm- gRmLm ù
ú
û
d
^
t
 
2
zm 

+
(e)
+
å
A < B 
QAQBgAB +
å
m 
gRmLm
^
t
 

zm 

2
 
+
(f)
+ 1
2

å
k < m 

å
tt¢ Î {r,l} 
dTkTm¢gtktm¢Tk tt¢
^
t
 

zk 


^
t
 

zm 

+
(g)
-
å
m 
brmlmRmLmd
^
t
 

+m 
+ ^
t
 

-m 

(h)
(0)
This representation allows the following family of approximate treatments for the energy. In the case the geminal amplitude related ESPs are fixed at their transferable values - the corresponding approximation is termed as the FA i.e. the fixed one - the energy Eq. (32) reduces to the lines (a) and (b), which yields an expression dependent on the molecular geometry and the hybridization ESPs only. The other lines in Eq. (32) reappear if corrections to the amplitude related ESPs are taken into account. This corresponds to the amplitude tuning in response to the geometry variations or environment details and this family of approximations is thus termed as the TA, i.e. the tuned amplitudes approximation. The simplest one in this family is that which retains only the terms linear in á [^(t)]zm ñ ~ mm (lines (c) and (d)) thus allowing the bond asymmetry (polar bonds and nonvanishing effective charges). Including the terms linear in zm -1 provides the possibility to take into account the variations of two-electron density matrices in response to geometry and environment variations (line (e)), but only including the quadratic terms yields corrections to the bond orders (line (h)). Further in this paper we mean by the TA approximation namely its simplest (mm-linear) version. Since the corrections to the bond orders appear only in the second order in small parameters mm and zm -1 the FA and the simplest TA approximations may have some value as it can be seen from Table by comparing columns (10) and (19), (20).

Further components of the description are those related to the HOs. The latter enter into the theory through the Hamiltonian matrix elements in the HO basis. The matrix elements entering Eqs. (7), (32) are either invariant with respect to basis transformations (the interatomic Coulomb interaction gAB) or can be uniquely expressed through contributions of s-AO to the HOs (the one-center matrix elements). The only class of molecular integrals depending on the whole structure of the HOs (including directions) is that of the resonance integrals. As we mentioned in Section 0.1 each spx-HO can be considered as a normalized quaternion (s,[(v)\vec]). Following [] we represent the entire system of HOs at any given atom by four vector parts [(v)\vec] m of the corresponding orthonormal quaternions. Even this representation is superfluous since only six Jacobi angles suffice to describe the system of HOs of each given atom completely. Nevertheless, usage of the vector parts is visual. If the latter are assumed to have the corresponding nucleus as their common origin the tetrahedral shape thus obtained contains (with an excess) all necessary information about the system of HOs of the given atom. In [] such a construct was called the hybridization tetrahedron of the heavy atom at hand. Using the hybridization tetrahedra as elements of the theoretical construct allows further discrimination of possible approximations. Due to the mentioned dependencies of the molecular integrals on the Jacobi angles both the FA and TA approximations to the energy Eq. (32) depend on the relative orientation of the hybridization tetrahedra through the bond resonance integrals brmlmRmLm. The resonance integrals depend also on the shapes (relative weights of the s- and p-contributions to the HOs which ultimately define the interhybrid angles) of the hybridization tetrahedra. All other terms in Eq. (32) in the FA and TA approximations depend only on the shapes of the hybridization tetrahedra. This leads to the possibility to either fix the relative weights of the s- and p-orbitals (FO i.e. fixed orbitals approximation) at spn (n = 1 ¸3) or any other allowable values and by this fix the shapes of the hybridization tetrahedra which thus become interacting rigid bodies or to allow the relative weights of the s- and p-orbitals to be tuned thus leading to the TO - tuned orbitals - picture of the flexible hybridization tetrahedra. Whichever combination of the FA or TA treatments for the density matrix elements on one hand with the FO or TO treatments for the HOs on the other hand results in a representation of the molecular energy Eq. (32) as such of the system of tetrahedral bodies (rigid or flexible) whose interactions and self energies depend on distances between their centers, their shapes and relative orientations. For example, the energy variation due to small pseudo- and quasirotations of the hybridization tetrahedron in the vicinity of the equilibrium for sp3-hybridized carbon atom in the FA approximation is given by a diagonal quadratic form:
dww(2)E = 2 æ
è
æ
è
d ®
w
 
Rm
b 
|GbbRmRm| d ®
w
 
Rm
b 
ö
ø
+ æ
è
d ®
w
 
Rm
l 
|GllRmRm| d ®
w
 
Rm
l 
ö
ø
ö
ø
,
GllRmRm = 4
Ö3
æ
è
bzsRmLmsmLm-bzzRmLm
Ö
 

1-(smLm) 2
 
ö
ø
Á,
GbbRmRm = 2[ æ
ç
è
bssRmLm+ 1
Ö3
bzsRmLm ö
÷
ø
smLm-
- æ
ç
è
bszRmLm+ 1
Ö3
bzzRmLm ö
÷
ø

Ö
 

1-( smLm) 2
 
]Á,
(0)
Three equal eigenvalues (diagonal elements) correspond to variations of [(w)\vec] b while three other equal eigenvalues correspond to those of [(w)\vec] l. Following the general reasoning [] this can be treated as a mechanistic model for molecular energy, for which the name "deductive molecular mechanics" has been coined in []. These constructs apparently satisfy the conditions formulated in the Introduction for a QM derived mechanistic model of molecular PES. However, there are several questions which have to be clarified. First, we have to see what might be the relation between the DMM and the standard MM known in the literature. Second, the precision of the the FA, TA, FO, and TO approximate treatments must be evaluated. Third, the interrelation between the ESPs and the molecular geometry must be established. The rest of the paper is devoted to investigating these three topics.

Relation between deductive and standard MM

General setting

The content of the deductive molecular mechanics as formulated in Ref. [] and above is a description of the molecular energy in the form of Eq. (32) as a function of shapes and mutual orientations of the hybridization tetrahedra and of geometry parameters. On the other hand the standard MM can be qualified as a scheme directly parameterizing the molecular energy as a function of molecular geometry only. From this point of view the Jacobi angles variables [(w)\vec] b, [(w)\vec] l describing the shapes and orientations of hybridization tetrahedra are superfluous and must be excluded. This can be done by finding the response of the corresponding ESPs to the variations of bond lengths and valence angles with use of linear response relations between different subsets of variables pertinent to the DMM picture. To do so let us consider a minimum of the energy with respect to both geometry and the ESPs. In the vicinity of a minimum the energy can be expanded upto second order with respect to nuclear displacements q and variations of the ESPs x:
E = E0+ 1
2
( x-x0| ÑxÑxE| x-x0) +
+( x-x0| ÑxÑqE| q-q0) +
+( q-q0| ÑqÑxE| x-x0) +
+ 1
2
( q-q0| ÑqÑqE| q-q0) ,
(0)
where linear terms disappear due to minimum conditions. For the sake of definiteness we restrict ourselves by the FA picture. Then the only remaining ESPs are the Jacobi angles describing the shapes and orientations of the hybridization tetrahedra. Minimization of the energy Eq. (34) with respect to x for a given value of q leads to basic linear response relation between the ESPs and the geometry distortions:
x-x0 = -( ÑxÑxE) -1ÑxÑqE(q-q0).
(0)
The scheme of this type was used in Ref. [] to demonstrate the transferability of scaling factors typically applied to reach an agreement between the calculated (ab initio) IR frequencies and the experimental ones.

Linear response relations for hybridization ESPs

The main use of the formulae Eq. (35) is for exclusion of the angular variables describing the hybridization tetrahedra from the mechanistic picture. Now we estimate the precision of the linear response relations (Eq. (35)) between geometry and hybridization variations themselves by numerical study of elongation of one C-H bond and deformations of valence angles. We consider tetrahedral methane molecule as a reference (its parameters then correspond to subscript 0 in Eqs. (34), (35)). First of all, we notice that the ÑxÑxE matrix further simplifies for methane since sLmm = 1 and, therefore, simple analytical expressions become possible. Also we remark that the FA approximation is adequate here since, for example, even very large elongation of one C-H bond by 0.1 Å leads to changes of the bond geminal amplitudes u,v, and w not exceeding 0.003. The same applies to the averages of the pseudospin ([^(t)]) operators.

Linear response of hybridization to bond elongation  

Let us consider first the relation between hybridization and elongation of the C-H bond. For this end we need the mixed second order derivatives coupling the bond stretching with the hybridization ESPs. For every C-H bond in methane we can introduce diatomic coordinate frame with the z axis directed along the bond and express the resonance integral as:
bCHmrmlm = bssCHsm+bzsCHvmz,
(0)
where the subscript m enumerates the C-H bonds. In the case of methane sm = [1/2] and vmz = [(Ö3)/2]. Changing the the bondlength causes the response of the vector d[(w)\vec] l to be exactly zero (vector product of collinear vectors) since the directions of the chemical bonds and the HOs coincide in the reference structure of methane. At the same time the response of the shape of the hybridization tetrahedron represented by the vector d[(w)\vec] b is nonvanishing and can be written as []:
d ®
w
 

b 
= - Ö3
4
· Ö3qssCH-qzsCH
Ö3bssCH+bzsCH
®
e
 

m 
drm,
(0)
where the derivative of bmn with respect to the interatomic distance is qmn, [(e)\vec]m is the unit vector directed along the CHm bond and drm is the variation of the length of the same C-H bond. The denominator in this expression corresponds to the eigenvalue of the ÑxÑxE matrix referring to [(w)\vec] b, while the numerator correponds to the block of the ÑxÑqE matrix where q is the bondlength rm and x is [(w)\vec] b.

Formula Eq. (37) gives the analytical expression for the coupling coefficient between the bond elongation (in Å ) multiplied by unit vector of this bond direction and changes of pseudorotation angles d[(w)\vec] b in methane (in radians). Its numerical value C1 is 0.2764 rad·Å-1. This distortion corresponds to the following form for the matrix of small transformation of the whole set of HOs (matrix H in Eq. (30)):
æ
ç
ç
ç
ç
è
1
-d
-d
-d
d
1
0
0
d
0
1
0
d
0
0
1
ö
÷
÷
÷
÷
ø
(0)
where the linear response estimate for d is C1·[(dr)/(Ö3)]. This form of the transformation matrix is perfectly numerically reproduced both in the FA and TA pictures. Fig. 1 represents the relative deviation of d estimated in the linear response approximation with respect to the exact value obtained by the energy minimization (both values are obtained within the FA picture). This figure shows that the linear response is good even for large distortions (deviation from linearity is only about 1% for stretching of 0.05 Å). The approximately linear dependence of the error itself on the variation of the bond length certifies that the second order estimate should perfectly describe the d parameter obtained by the energy minimization.

Linear response of hybridization to valence angle deformation  

The linear response relations between the molecular shape and the shape of hybridization tetrahedron are rather tricky due to complex structure of the hybridization manifold. The molecular shape can be formally characterized by unit vectors with origin at an atom considered, taken as a center, and pointing to those bonded to the central one. In the case of methane the deformations of thus defined coordination polyhedron are small rotations of unit vectors [(e)\vec]m directed from the carbon atom to hydrogen atoms. Their small rotations d[(j)\vec] m form an 8-dimensional space which decomposes to a direct sum of two subspaces: one 3-dimensional corresponding to rotations of the molecule as a whole and another 5-dimensional corresponding to independent variations of valence angles. The former one is precisely mapped on the 3-dimensional space of quasirotations d[(w)\vec] l while the latter (5-dimensional) one is mapped on the 3-dimensional space of pseudorotations d[(w)\vec] b corresponding to changes of the shape of the hybridization tetrahedron []. Due to very general theorems of linear algebra [] there exists at least a two-dimensional kernel in the space of deformations of molecular shape which maps to zero deformation of hybridization tetrahedron. In [] the term "hybridization incompatible" has been coined for the deformations from this kernel. The structure of deformations laying in the kernel of the mapping is quite simple: they are produced by equal variations of opposite (spiro) valence angles. In contrast, the variations which correspond to increase of one valence angle by dc and decrease of its spiro counterpart by the same value fall into "coimage" of this mapping i.e. to the subspace which one-to-one maps to the space of pseudorotations d[(w)\vec] b. The deformations in the coimage can be called "hybridization compatible". It is clear that only these latter variations should be considered. It is also clear that any variation of the valence angle is a sum of equal amounts of hybridization compatible and hybridization incompatible deformations. The denominator in the linear response relation Eq. (35) is the same as for Eq. (37) while the relevant block of the ÑxÑqE matrix (with q taken as a difference of two opposite valence angles) is proportional to bzsCH. Applying the linear response technique to the "hybridization compatible" variation of two spiro valence angles we obtain:
d ®
w
 

b 
= - bzsCH
Ö2(Ö3bssCH+bzsCH)
dc ®
k
 
,
(0)
where [(k)\vec] is the ort directed along the symmetry axis preserved during this distortion [].

The coupling coefficient between the change of the pseudorotation vector and totally hybridization compatible deformation of valence angles can be easily found by Eq. (39). Its numerical value C2 is -0.20734 for the equilibrium interatomic distance. The considered distortion produces the following HO transformation matrix (matrix H in Eq. (30)):
æ
ç
ç
ç
ç
è
1
0
0
-d
0
1
0
0
0
0
1
0
d
0
0
1
ö
÷
÷
÷
÷
ø
,
(0)
where the linear response estimate for d is C2·dc. The change of angle between HOs (i.e. between the vector parts [(v)\vec]m of the quaternions representing the HOs) is -2Ö2d/3 for small values of d. The above form of the HO transformation matrix is perfectly confirmed by our numerical calculations performed within the FA picture even for very large distortions which is a consequence of the mathematical structure of the hybridization manifold described above. Fig. 2 represents the relative deviation between the linear response estimate of the d parameter in Eq. (40) and its value coming from the energy minimization procedure. These data show that the linear response estimate performs very well even for improbably large distortions (the deviation from the linear responce estimate is smaller than 0.25% for the distortion of 0.3 rad (about 17°)).

The smallness of the coupling coefficient C2 even for the totally hybridization compatible deformations allows to qualitatively understand certain features of the electronic structure of cyclopropane as it appears in the SLG approach: a very large distortion of the C-C-C valence angle from the tetrahedral to 60° one leads only to a relatively small distortion of the corresponding interhybrid angle. We model this process by strongly deforming the methane molecule. Simple estimate is based on Eqs. (39), (40) and runs as follows. The valence angle deformation when going from methane to cyclopropane is of 49.5° (=109.5°-60°); only one half of it is hybridization compatible; after multiplying by C2 this yields the value of the interhybrid angle between the HOs corresponding to the üntouched" C-H bonds of 114.5° (i.e. the angle variation amounts only 5°). From the energy minimum condition for hydrides it follows that the C-H bonds indeed must follow the directions of the HOs. Numerical experiments performed with use of the SLG-MINDO/3 method show that if one of the H-C-H valence angles is fixed at the cyclopropane value of 60° the energy minimum corresponds to its spiro counterpart of 115°. These results can be directly compared with the experimental H-C-H valence angle in the cyclopropane molecule which equals to 115.1°.

Analogous estimate can be applied to cyclobutane. In this case we consider the distorted methane molecule with angle 90°. The response of the HOs to the deformation is proportional in our model to deviation of the valence angle from the tetrahedral one. The deviation of the C-C-C angle from the tetrahedral one in cyclobutane (19.5°) amounts 40% of that in cyclopropane. Therefore, we can expect that about the same ratio will be observed for the deviations of the H-C-H valence angle from the tetrahedral one in the cyclobutane and cyclopropane molecules. In fact, this ratio in the SLG-MINDO/3 numerical experiment is about 39%.

Estimates of parameters of the standard MM force fields based on DMM

Announced transition from the DMM model of molecular PES to a model dependent on molecular geometry is formally obtained by inserting Eq. (35) to Eq. (34) which yields
E = E0+ 1
2
( q-q0| ÑqÑqE| q-q0) -
- 1
2
( q-q0| ÑqÑxE( ÑxÑxE) -1ÑxÑqE| q-q0) .
(0)
It apparently consists of two contributions: (i) the leading one - the second derivatives of the energy Eq. (32) with respect to the geormetry parameters q, and (ii) a smaller correction appearing as a result of projecting out the ESPs related to the hybridization tetrahedra. Due to the SLG form of the wave function the energy expression Eq. (32) is naturally represented as a sum of atom and bond increments. Staying for the sake of simplicity within the FA picture we can conclude that the only geometry dependent contribution is that proportional to the resonance integrals of respective bonds. These contributions depend namely on the natural nuclear coordinates []: bond lengths and valence angles entering the definition of the standard MM force fields. This makes sensible to consider the geometry dependence of individual bond energies. For a pair of singly bonded atoms the corresponding energy is:
ERm+ELm+ERmLmbond+ERmLmnonbond
(0)
which for a symmetric bond may be written as:
E = 2(Um-gRmLm)+ gm
2
æ
ç
è
1- 1
G(zm)
ö
÷
ø
+
+ gRmLm
2
æ
ç
è
1+ 1
G(zm)
ö
÷
ø
-2brmlmRmLm zm
G(zm)
+ ZRmZLmDRmLm,
(0)
where Um is the mean arithmetic value of the one-center electron core attraction parameters Umr and Uml.

We consider in more details the energy curve for the C-H bond. The curve corresponding to the sp3 hybridization of carbon atom and to the symmetric TA picture (Eq. (20)) is given by Fig. 3. It has correct qualitative behaviour for all interatomic separations. The minimum depth on this curve is approximately -0.23 a.u. and can be considered as the "pure" energy of the C-H bond. It is generally accepted in the literature that the energy of C-H bond is approximately 0.15 a.u. At the same time the latter value is thermodynamic one while our value is obtained by extracting the contributions to the energy intrinsic to this bond and excluding the interaction between the bonds. The difference between the thermodynamic value for the bond energy and that obtained from the SLG energy in the FA picture can be explicitly written in quite simple form:
1/4[Us(C)+3Up(C)+3gtt¢C+6DHH-EA(C)],
(0)
where EA(C) is the energy of non-hybridized carbon atom. Adding this value to the minimum of the energy curve for the C-H bond gives the value close to the thermodynamic one since the heats of formation are well reproduced within the SLG-MINDO/3 method.

It is interesting to compare the form of the bond energy curve Fig. 3 with the Morse potential. To this purpose we tried to approximate the curve of Fig. 3 by the Morse function D0[1-exp(-a(r-re)/re)]2 by minimizing the area between two curves in the interval from 0.72 Å to 2.50 Å . With the parameters D0 and re fixed at the values equal to the minimum depth and position on the curve (0.2295 a.u. and 1.078 Å ) the optimal value of parameter a is then 2.306 but with these parameters two curves are in fact quite different (the area between curves is almost 11% of area between the bond energy curve and the abscissa). If we optimize all three parameters of the Morse curve they become slightly modified D0=0.2333 a.u., re=1.045 Å , and a=2.295. This reduces the area between the curves by 30%. It should be concluded that the energy profile in the TA approximation is not partcularly well reproduced by whatever Morse curve.

In order to estimate the parameters of harmonic force fields we consider the symmetric correlated single bond, where the energy can be obtained without any reference to its environment. In our case the derivative of the bond energy with respect to a geometry parameter q has the form:
Em
q
= ZRmZLm DRmLm
q
-2 zm
G(zm)
brmlmRmLm
q
- 1
2
æ
ç
è
1- 1
G(zm)
ö
÷
ø
gRmLm
q
,
(0)
where the derivatives of different ESPs with respect to geometry variables exactly cancel each other so that the final expression for the energy derivative acquires the form expected from the perturbative analysis of the ESPs given above. If q is the interatomic distance setting the derivative to zero yields the equiation for determining the minimum position. Results are given in Table . In the limit zm >> 1 we recover the equilibrium geometry condition for the FA picture. The meaning of other notation in the Table is following. The TAsymm refers to a TA estimate for the symmetric bond (bond asymmetry - polarity - terms omitted), while TApert refers to the perturbative inclusion of bond asymmetry effects to the TA picture with use of the m0 parameter. All estimates are quite reasonable. At the same time the latter one looks more promising because it perfectly well corresponds to the equilibrium interatomic separation in methane and other hydrocarbons.

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)
We see from Table that at least for one of the ëxperimental" estimates for the stiffness of the C-H bond [] the agreement is quite acceptable. The deviation from other cited values may be understood since the bond stretch parameters fitted by other authors in the context of the structure oriented MM schemes are implicitly loaded with the average effects of surrounding atoms. This disagreement cannot be ascribed to the effects of the projected out Jacobi angles. Indeed, the off-diagonal constant coupling stretching for two incident C-H bonds in the methane molecule can be written as:
Koff = 1
4Ö3
(Ö3qssCH-qzsCH)2
Ö3bssCH+bzsCH
(0)
The only reason why this term appears is the deformation of the carbon hybridization tetrahedron Eq. (37) effectively coupling stretchings of two C-H bonds. Its estimated magnitude is only 0.120 mdyn/Å which is in an agreement with its IR estimate 0.03 mdyn/Å [] upto order of magnitude. This estimate however establishes the scale of the corresponding effects. One can see that the DMM specific corrections to the bare estimates of the harmonic stretching constants Eq. (46) are small. Nevertheless in the cases when the bare harmonic constant vanishes (like the considered off-diagonal constant does) the DMM correction gives an estimate correct to the order of magnitude and allows to solve the question on the presence of the off-diagonal terms on purely theoretical basis.

Analogous treatment of the energy terms quadratic in valence angles' deformations yields the bare estimate for the harmonic bending constant in the form:
kHCH = Ö3
2
bzsCH,
(0)
From Table one can see that the elasticity constants for bending force fields are in good agreement with the values accepted in the literature.

Discussion

In the previous Section we provided the exclusion of the angular variables characterizing the shapes and orientations of the hybridization tetrahedra from the mechanistic DMM model of molecular PES. This results in a model announced in the Introduction, which is similar to the standard MM models but is obtained by the sequential derivation from the QM (SLG) model of molecular electronic structure. As it is mentioned the transferability of the ESPs characterizing chemical bonds in molecules and linear response relations for hybridization ESPs are main components of deriving MM theory of molecular PESs from corresponding QM theory. Both these features have been mathematically derived and numerically checked in Section 1.

Despite its long history the very term "transferability" remains somewhat vaguely defined synonym of äll the best" in parameterization schemes, referring largely to their capacity to be used without change for any molecule in a sufficiently wide class of similar ones. From quantitative point of view this concept have got some attention in two related areas. First we mention the estimates of transferability given in Ref. [] where that of the semiempirical quantum chemical parameters has been related to the fact that the corresponding quantities remain the same for all molecules of similar structure upto the second order with respect to overlap integrals between AOs residing at neighbour atoms. That allowed to define the transferability for the quantum chemical parametrs (ultimately, for the Hamiltonian matrix elements) as invariance of some quantity to a given order of precision with respect to a small parameter. Analogously in Ref. [] the problem of constructing transferable dynamical matrices in relation to analysis of vibrational spectra has been considered. The stability of the dynamical matrix was analyzed with respect to small parameter of relative mass variation under the isotope substitution in a series of related molecules.

The importance of the transferability of the geminals has been pointed out yet in []. It was stated that the assumption of the transferability of the geminal amplitudes is a prerequisite for that of the bond energy. However, in [] the geminal transferability had not been shown and the authors concentrated on the statements equivalent to the transferability of the MM bond stretching force fields. The transferability of the MM force fields must be considered as an important chracteristic of this approach. The reasons to treat that or another force field as a transferable between two specific molecules or classes of molecules are either purely pragmatic or this question is solved on the school-wise grounds []. As far as we know in the literature there were no attempts reported to prove this property of the force fields from any general point of view. Here we present a step towards quantitative analysis of the transferability of the MM force fields by proving the transferability of the density matrix elements - equivalent to that of the geminal amplitudes, but more directly related to the energy. Under the assumptions given by Eqs. (15), (16) the averages of the pseudospin operators (and thus all the bond ESPs) are invariant in that sense that they do not depend on environment of the bond under consideration and even on particular composition of the bond, i.e. on the nature and the hybridization of the atoms the bond connects. This corresponds to the FA picture (see above and Ref. []). It is important that the invariance (at the established level of precision) of the density matrix elements can be proven only for the basis of the variationally determined HOs - a specific characteristic of the SLG approach [,]. In the basis of AOs the density matrix elements are not invariant even approximately. Though the approximation sufficient to obtain formally these invariant results (the SCF approximation) is very crude it, nevertheless, breakes only at large interatomic separations which normally are not covered by any MM-like approximation. This result allows to pose further questions: to what extent the density ESPs' invariance may stand further improvements of the description and whether it is possible to relate the invariance of the density matrix elements with the transferability of the MM force fields. To answer these questions we notice that the invariant values of ESPs can be improved by perturbative corrections (the TA picture) reflecting all diversity of chemical compositions and environments the bond may occur in. Nevertheless, all the variety of perturbations is characterized by two small dimensionless parameters: zm-1 Eq. (21) and mm Eq. (23). Both parameters depend on the atoms connected by the bond, their separation, and their hybridization. The perturbative treatment allows to estimate the precision of transferability. For example, using Eqs. (22), (24), and (25) we conclude that the bond order is the quantity transferable upto second order with respect to both zm-1 and mm; the ionicity (the total weight of the ionic configurations) is transferabile upto second order with respect to mm and upto first order with respect to zm-1; the bond polarity is transferable upto first order with respect to both zm-1 and mm. The second order transferability of bond orders explains to certain extent the success of the concept of ßingle bond" suitable for a large variety of chemical bonds. Note that the second order transferability takes place for the bond orders also in the case when we employ the SLG bond wave function with the correct asymptotic behavior despite the fact that the transferable numerical value itself is obtained from the SCF wave function which does not have the correct asymptotic behavior. Within this picture all specific characteristics of the force field are loaded into parameters of the (effective) Hamiltonian, which are numbers specific either for a given atom in certain hybridization state or for a pair of such hybrid states of atoms - ends of the bond. Although the density ESPs can be considered as constants independent on any details of molecular composition or geometry, the force fileds which are basically sums of products of ESPs by matrix elements of molecular Hamiltonian are geometry dependent and composition specific. The force fileds thus obtained are expected to be the same for the same composition of the bond and to depend weakly (to the extent of the variance of the mm1 parameters) on the environment. These properties are basically much more than necessary for substantiating whichever MM-like description.

Conclusion

In the present paper we discussed the problem of deriving the MM representation of the molecular PES from a relevant QM description. Using the SLG wave function we analyzed the ESPs related to bond geminals and to hybridization tetrahedra. It was shown that the bond-related parameters can be represented as functions of parameters of the MINDO/3 Hamiltonian in the HO basis, transferable from one molecule to another. The functional form of the ESPs found is valid at arbitrary interatomic separations. At the interatomic separations close to the equilibrium bond lengths characteristic for the MM-like treatments two approximations, both suitable for substantiation of the ESPs transferability were considered. One is the fixed geminal amplitudes approximation which results in perfectly transferable numbers referring to the ESPs in question. Another, more exact is the tuned geminal amplitudes approximation which takes into account small corrections to the invariant ESPs. Two small parameters characterizing specificity of the bond and effects of its environment were introduced. By this the whole manifold of quantum chemical parameters defining the effective bond Hamiltonian boils down to only two relevant parameters zm -1 and mm. The presence of such only two-dimensional manifold and smallness of the parameters for a wide range of bonds in quite different environments essentially explains the transferability of the density related ESPs. This allows for a family of mechanistic models describing molecular PESs in terms of hybridization tetrahedra with interactions dependent on distances between their centers and on mutual orientations. Linear response relations for variation of hybridization parameters due to elongation of chemical bonds or specific changes of valence angles are considered. They allow to exclude the angular variables describing the shapes and orientations of hybridization tetrahedra and to represent the molecular energy in both the FA and TA approximations as that of the system of interacting point masses [] ("balls-and-springs" picture) depending on the molecular geometry only. This energy has a form of a sum of local (bond) increments corresponding to the force fields of the standard MM. The estimates of the paramemters of these force fields coming from the analytical expressions are compared with those obtained in numerical experiments showing the high accuracy of analytical estimates. The reasons for this possibility are both the transferability of the bond-related ESPs and the linear response relations for the hybridization tetrahedra are numerically tested for their precision.

Acknowledgements

This work was performed with financial support of the RFBR through the grants 02-03-32087, 04-03-32146, and 04-03-32206. It has been completed during the stay of A.M.T. in the RWTH, Aachen in the frame of the Alexander von Humboldt Postdoctoral Fellowship which is gratefully acknowledged as is the kind hospitality of Prof. R. Dronskowski. A.L.T. gratefully acknowledges valuable discussions with Profs. N.F. Stepanov, A.A. Levin, and I. Mayer.

References

[]
V.G. Dashevskii, Conformations of Organic Molecules, Khimiya, Moscow, 1974 [in Russian].

[]
U. Burkert and N. Allinger, Molecular Mechanics, ACS, Washington, DC, 1982.

[]
S. Goedecker, Rev. Mod. Phys., 71 (1999) 1085.

[]
S.Y. Wu and C.S. Jayanthi, Physics Reports, 358 (2002) 1.

[]
A. Warshel and M. Levitt, J. Mol. Biol., 103 (1976) 227.

[]
J. Gao, In: Reviews in Computational Chemistry, VCH Publishers, Inc., New York, 1996.

[]
P. Sherwood, In: Modern Methods and Algorithms of Quantum Chemistry, Proceedings, Second edition, Ed. J. Grotendorst, NIC Series, 3, John von Neumann Institute for Computing, Jülich, 2000.

[]
A.L. Tchougréeff and A.M. Tokmachev, In: Advanced Topics in Theoretical Chemical Physics, Eds. J. Maruani, R. Lefebvre, and E. Brändas, Kluwer, Dordrecht, 2003.

[]
A.L. Tchougréeff, Phys. Chem. Chem. Phys., 1 (1999) 1051.

[]
A.M. Tokmachev, A.L. Tchougréeff, and I.A. Misurkin, J. Mol. Struct. (Theochem), 506 (2000) 17.

[]
A.M. Tokmachev and A.L. Tchougréeff, Int. J. Quantum Chem., 84 (2001) 39.

[]
N.L. Allinger, J. Am. Chem. Soc., 99 (1977) 8127.

[]
N.L. Allinger, Y.H. Yuh, and J.-H. Lii, J. Am. Chem. Soc., 111 (1989) 8551.

[]
W.L. Jorgensen and J. Tirado-Rives, J. Am. Chem. Soc., 110 (1988) 1657.

[]
B.R. Brooks, R.E. Bruccoleri, B.D. Olafson, D.J. States, S. Swaminathan, and M. Karplus, J. Comp. Chem., 4 (1983) 187.

[]
W.D. Cornell, P. Ciepak, C.I. Bayly, I.R. Gould, K. Merz, D.M. Ferguson, D.C. Spellmeyer, T. Fox, J.W. Caldwell, and P. Kollman, J. Am. Chem. Soc., 117 (1995) 5179.

[]
T.A. Halgren, J. Comp. Chem., 17 (1996) 490.

[]
N.L. Allinger, K. Chen, and J.-H. Lii, J. Comp. Chem., 17 (1996) 642.

[]
P.Comba and T.W. Hambley, Molecular Modeling, VCH Publishers, Weinheim et al., 1995.

[]
J. Kozelka, In: Metal Ions in Biological Systems, Vol. 33, Eds. A. Sigel and H. Sigel, Marcel Dekker Inc, New York et al., 1996.

[]
I.V. Pletnev and V.L. Melnikov, Russ. Chem. Bull., 46 (1997) 1278.

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

[]
S. Diner, J.-P. Malrieu, and P. Claverie, Theor. Chim. Acta, 13 (1969) 1.

[]
T.L. Allen and H. Shull, J. Chem. Phys., 35 (1961) 1644.

[]
A.M. Tokmachev and A.L. Tchougréeff, J. Comp. Chem., 22 (2001) 752.

[]
A.M. Tokmachev and A.L. Tchougréeff, J. Phys. Chem. A, 107 (2003) 358.

[]
A.M. Tokmachev and A.L. Tchougréeff, Int. J. Quantum Chem., 85 (2001) 109.

[]
A.M. Tokmachev and A.L. Tchougréeff, Int. J. Quantum Chem., 88 (2002) 403.

[]
A.L. Tchougréeff and A.M. Tokmachev, Int. J. Quantum Chem., 96 (2004) 175.

[]
A.L. Tchougréeff and A.M. Tokmachev, Int. J. Quantum Chem., 100 (2004) to appear.

[]
A.L. Tchougréeff, J. Mol. Struct. (Theochem), 630 (2003) 243.

[]
V.A. Fock, M.G. Veselov, and M.I. Petrashen', J. Exp. Theor. Phys., 10 (1940) 723 [in Russian].

[]
S. Weinbaum, J. Chem. Phys., 1 (1933) 593.

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

[]
J.M. Parks and R.G. Parr, J. Chem. Phys., 28 (1958) 335.

[]
A.A. Abrikosov, L.P. Gor'kov, and I.Ye. Dzyaloshinskii, Quantum Field Theoretical Methods in Statistical Physics, Pergamon, Oxford, 1965.

[]
H. Poincaré, La Science et l'Hypotése, Flammarion, Paris, 1923 [in French].

[]
V.I. Pupyshev, Y.N. Panchenko, C.W. Bock, and G. Pongor, J. Chem. Phys., 94 (1991) 1247.

[]
O. Ermer and S. Lifson, J. Am. Chem. Soc., 95 (1973) 4121.

[]
S. Chang, D. McNally, S. Shary-Tehrany, M.J. Hickey, and R.H. Boyd, J. Am. Chem. Soc., 92 (1970)3109.

[]
R.G. Pearson, Symmetry Rules for Chemical Reactions, Wiley, New York, 1976.

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

[]
A.I. Kostrikin and Yu.I. Manin, Linear Algebra and Geometry, Nauka, Moscow, 1986 [in Russian].

[]
I. Fisher-Hjalmars, In: Modern Quantum Chemistry. Istanbul Lectures. Vol. 1, Ed. O. Sinanoglu, Academic Press, New York et al., 1965.

[]
N.F. Stepanov, G.S. Koptev, V.I. Pupyshev, and Yu.N. Panchenko, In: Modern Problems of Physical Chemistry, Vol. 11, MSU Publ., Moscow, 1979 [in Russian].

[]
A.Y. Meyer, In: Theoretical Models of Chemical Bonding, Part 1, Ed. Z.B. Maksi\'c, Springer, Heidelberg, 1989.

Table 0: Estimates of parameters of bond geminals
MoleculeBondz-1m á [^(t)]zm ñ á [^(t)]zm2 ñ á [^(t)]+m ñ
(10)(19) (24)(10)(19) (20)(25) (10)(20) (22)(25)
H2 H-H0.1270.0000.0000.0000.0000.4370.4360.4370.4370.9920.9920.9920.992
CH4 C-H0.1810.0930.0650.0950.0650.4140.4100.4110.4140.9820.9840.9840.982
NH3 N-H0.1720.0970.0690.0980.0690.4190.4140.4150.4190.9830.9860.9850.983
H2O O-H0.1610.3410.2440.3450.2470.4640.4190.4200.4660.9590.9870.9870.959
HF F-H0.2850.5400.3090.5620.3080.4510.3570.3630.4570.9250.9620.9590.929
C2H6 C-C0.1930.0000.0000.0000.0000.4050.4040.4050.4050.9820.9820.9810.982
C-H0.1810.0710.0500.0720.0500.4130.4100.4110.4130.9830.9840.9840.983
C-C0.1930.0160.0110.0160.0110.4050.4040.4050.4050.9820.9820.9810.982
C3H8 C1-H0.1810.0740.0520.0750.0520.4130.4100.4110.4130.9830.9840.9840.983
C2-H0.1800.0510.0350.0510.0350.4120.4100.4110.4120.9840.9840.9840.984
Cyclo- C-C0.1990.0000.0000.0000.0000.4020.4000.4020.4020.9810.9810.9800.981
propane C-H0.1800.0840.0590.0860.0590.4140.4100.4110.4140.9830.9840.9840.983
N2H4 N-N0.2900.0000.0000.0000.0000.3610.3550.3610.3610.9600.9600.9580.960
N-H0.1760.0920.0650.0930.0650.4160.4120.4130.4160.9830.9850.9840.983
N-C0.2190.0250.0160.0260.0160.3930.3900.3930.3930.9770.9770.9760.977
CH3NH2N-H0.1730.0980.0690.0990.0690.4190.4120.4150.4190.9830.9850.9850.983
C-H0.1870.0780.0530.0790.0530.4100.4060.4080.4100.9820.9830.9820.982
O-C0.2180.2760.1780.2820.1790.4200.3910.3940.4210.9640.9770.9760.964
CH3OH O-H0.1580.3280.2360.3320.2400.4620.4210.4220.4640.9620.9880.9880.961
C-H0.1820.0770.0530.0780.0530.4130.4090.4110.4130.9830.9840.9840.983
CH3F F-C0.3550.4170.2120.4420.2080.3820.3220.3330.3830.9300.9420.9370.932
C-H0.1800.0890.0620.0910.0620.4140.4100.4110.4140.9820.9840.9840.982

Table 0: Parameters of bond geminals estimated with m0
MoleculeBondm0m1 á [^(t)]zm ñ á [^(t)]zm2 ñ á [^(t)]+m ñ
(10)(24)(10)(25)(10)(25)
H2 H-H0.0000.0000.0000.0000.4370.4370.9920.992
CH4 C-H0.0750.0180.0650.0520.4140.4130.9820.983
NH3 N-H0.0910.0060.0690.0650.4190.4190.9830.984
H2O O-H0.343-0.0020.2440.2490.4640.4660.9590.959
HF F-H0.5400.0000.3090.3080.4510.4570.9250.929
C2H6 C-C0.0000.0000.0000.0000.4050.4050.9820.982
C-H0.0680.0030.0500.0470.4130.4130.9830.983
C-C0.0070.0090.0110.0050.4050.4050.9820.982
C3H8 C1-H0.0680.0060.0520.0470.4130.4130.9830.983
C2-H0.061-0.0100.0350.0420.4120.4130.9840.983
Cyclo- C-C0.0000.0000.0000.0000.4020.4020.9810.981
propane C-H0.092-0.0070.0590.0640.4140.4150.9830.982
N2H4 N-N0.0000.0000.0000.0000.3610.3610.9600.960
N-H0.094-0.0020.0650.0660.4160.4170.9830.983
N-C0.033-0.0070.0160.0210.3930.3930.9770.977
CH3NH2N-H0.102-0.0050.0690.0730.4190.4190.9830.983
C-H0.078-0.0010.0530.0540.4100.4100.9820.982
O-C0.297-0.0210.1780.1930.4200.4250.9640.962
CH3OH O-H0.346-0.0180.2360.2530.4620.4690.9620.958
C-H0.093-0.0160.0530.0650.4130.4140.9830.982
CH3F F-C0.431-0.0150.2120.2150.3820.3870.9300.932
C-H0.092-0.0030.0620.0640.4140.4150.9820.982

Table 0: DMM based estimates of the MM force field parameters as compared to those accepted in some standard MM parameterizations
r0CH kCH kHCH
Å mdyn/Å mdyn/deg
FA: 1.069 8.30 0.509
TAsymm: 1.078 7.77
TApert: 1.096 7.17
Standard MM:
[]: 1.113 []: 4.5 ¸ 4.7 []: 0.549
[]: 1.105 []: 5.31 []: 0.508
[]: 1.090 []: 7.90 []: 0.493

Figure

Figure 3: Deviation of the linear response estimate of the d angle Eq. (38) from the precise one for elongation of the C-H bond in methane

Figure

Figure 3: Deviation of the linear response estimate of the d angle Eq. (40) from the precise one for angular distortion of methane.

Figure

Figure 3: C-H bond energy profile for methane molecule as obtained within the symmetric TAFO approximation.


Footnotes:

1AvH Postdoctoral Fellow, on leave from the Karpov Institute of Physical Chemistry, Moscow, Russia


File translated from TEX by TTH, version 2.67.
On 18 Aug 2004, 13:36.