Towards a possible ab initio molecular mechanics.
I. Stability of density matrix elements^{1}
May 3, 2007
A.L. Tchougréeff,
A.M. Tokmachev
I. Mayer
Abstract
Molecular structures are characterized by a large degree of additivity and
transferability of various intramolecular interactions determining the shape
and energetics of molecules. This property is constantly utilized by the
different molecular mechanics" (MM) schemes which allow
one to obtain quite reliable molecular geometries and relative energy values
for a wide range of molecular systems, which is especially remarkabke in the
light of the simplicity of the assumptions made. Numerous MM schemes
presented in the literature use different sets of parameters (force fields
etc.), which are adjusted empirically. The known success of MM
models poses two important questions: first, one wishes to understand why do
they work at all, and, second, one would like to develop schemes, in which
the parameters of MM can be determined theoretically. Such an analysis could
also give some deeper insight permitting to predict whether the given MM
scheme is expected to be successful if applied to the class of problems
actually at hand. From a more pragmatic point of view, having a bridge
connecting a quantum mechanical (QM) description of molecular structure with
a classical moldel (MM) can help to improve hybrid QM/MM methods by
providing a systematic derivation of the form of the junction between the
subsystems treated by MM and QM, respectively and by giving a priori
estimates of the junction parameters.
Our previous studies based on the semiempirical quantum chemical
Hamiltonians of MINDO and MNDO types permitted us to draw the conclusion
that the success of MM can indeed be understood on the basis of quantum
mechanics. The analysis of the geminaltype wave functions constructed by
making use of oriented hybrid orbitals showed that the parameters of the
wave function and the energy contributions of the individual bonds are
indeed well transferable, in good agreement with the simple chemical picture
of the systems studied.
It seems to be desirable to develop a similar analysis also at the ab initio level of the theory, because that would connect together the
chemical and physical description of molecules and explain the observed
transferability of molecular interactions more rigorously than semiempirical
theories can do it. On the other hand this treatment might be also useful
for developing QM/MM junctions in the cases when the QM part of the systems
are described at the ab initio level. We report here results of the
first step of this analysis: the geminal parameters expressed in a
symmetrically orthogonalized optimized minimal basis set exhibit a degree of
stability similar to that observed in the semiempirical case.
Introduction
Chemistry is dominated by the idea of transferability either explicitly
formulated or implicitly assumed. The concept of isolated chemical bonds 
which are to a great extent independent of their surrounding  allowed the
chemists to formulate a consistent approach which was remarkably successful
already in the XIX century in interpreting and predicting different
structural and reactivity features of various compounds, especially organic
ones. All these achievements were reached not only without any relation to
the quantum mechanical basics of the molecular electronic structure but even
without the knowledge of electrons themselves, not to speak about the
detailed understanding of the laws governing their motion. The modern
incarnation of that more than a century old approach is represented by the
various molecular mechanics" (MM) methods providing direct
parameterization schemes for molecular potential energy as function of the
coordinates of nuclei of the atoms composing the molecule Dashevskii,BurkertAllinger1982. Although the standard paradigm in the MM
area is that there is no direct relationship between the quantum theory of
molecular structure and the MM force fields, the entire MM theory is tacitly
based upon the concept of the isolated twocenter chemical bonds, which is
of course a crude representation of a basically quantum mechanical nature of
the ground state manyelectron wave function of the molecule (see below).
Meanwhile, modern quantum chemistry gives very little for understanding of
the fundamental transferability properties mentioned above. At best, they
are reproduced in the calculations (numerical experiments) perfomed at
different levels of sophistication, so that these important features remain
experimental facts still missing whatever theoretical explanation. However,
as stressed by Coulson, any explanation why must be
given in terms of concepts which are regarded as adequate or suitable. So
the explanation must not be that the electronic computer shows that D(HF) >> D(HH), since this is not an explanation at all,
but merely a confirmation of experiment" []. Also a
theoretical understanding of the origin of transferability observed in
experiments of both real and numerical types could be of a siginificant
practical use in the context of the hybrid QM/MM methods, since it would
potentially provide the form of the junction between the parts of a complex
molecular system treated on the ab initio QM and the classical MM
levels, respectively.
In the recent series of papers [,,] an attempt has been
undertaken to reconcile the quantum mechanical description with the idea of
transferable bonds. That was done in the framework of a semiempirical
quantum chemical model of molecular structure, in which one uses implicitly
orthogonalized valence AOs and the Hamitonian is parametrized against
experimental data. Technical details of these papers will be referred to
throughout the present paper. But here we just want to give some very
general notes concerning the overall structure of the theory which may
provide an explanation for the bonds' transferability and represent a basis
for a further development of a classical picture of the molecular structure,
compatible with its quantum mechanical description.
First of all we notice that the transferability of bonds must be a feature
pertinent to the wave function, which is built up as an (antisymmetrized)
product of welldistinguishable twoelectron components. This can be the
only basis of the semiobservability of the chemical bonds, as defined by
Ruedenberg [], which is tacitly assumed within the entire
MM picture and more generally  in chemistry itself, where applies. It is
important to mention here that for the molecular systems whose electronic
structure does not fit into the picture of twocenter twoelectron bonds,
the models of the MM force field type cannot be easily developed either. It
is noteworthy that the molecular Hamiltonian is a rather nontransferable
entity: those corresponding to different molecules differ by singular terms
of electronnuclear attraction in the direct physical space; in this respect
they are not close in any reasonable mathematical sense. Similar
consideration can be applied to the matrix form of the electronic
Hamiltonian: those for different molecules act in spaces of different
dimensionality and their matrix elements are geometry dependent. Thus there
are no a priori signs of the transferabilty of the Hamiltonian as
far as the description in the direct physical space is concerned. (By using
a refined analysis, some transferable elements of the Hamiltonian could be
identified in the framework of a special mixed" second quantization formalism [].) By contrast, formal
characteristics of the transferability, which are pertinent to the wave
functions, have been established in the framework of the semiempirical
theory (see Refs. [,,] and below).
The minimal requirement to a subsystem within a molecular system to be
distinguishable or observable is that it has a stable composition. In this
respect we may ask that to what extent the numbers of particles residing in
different parts of molecular system are well defined, and are, therefore,
good quantum numbers. To find this out we have to consider the dispersion
(or meansquare deviation) of the numbers of particles in the parts.
Let us introduce operators [^N]_{M} of the numbers of electrons in some
parts M of the molecular system such that:
where [^N] is the operator of the number of electrons in the entire
molecular system under consideration. For the normalized wave function Y we get:
áY 
^ N

M

Yñ = n_{M} , "M . 
 (0) 
As we mentioned above, the physical distinguishability or observability of a
part of the system is formally expressed as vanishing fluctuation
(dispersion) of the number of electrons in it. The dispersion of the number
of electrons in the Mth part:
áY( 
^ N

M

n_{M}) ^{2}Yñ 
 (0) 
must be as small as possible.
The minimum condition for the overall dispersion

min
 
å
M

áY( 
^ N

M

n_{M})^{2}Yñ 
 (0) 
can be trivially satisfied by the function Y taken as the
antisymmetrized product of group functions (GF):
where each of the group functions g_{M}^{+} is an eigenfunction of the
corresponding electron number operator [^N]_{M}:

^ N

M

g_{M}^{+} 0
ñ = n_{M}g_{M}^{+}0
ñ . 
 (0) 
This implies that n_{M} is an integer number. The above condition can be
satisfied if one assumes that the entire space L of oneelectron states of
the molecular system is divided into disjunct subspaces L_{M}:


L_{M}ÇL_{N} = {0}, M ¹ N. 



 (0) 
which are also assumed mutually orthogonal. Each of the subspaces L_{M} serves as a carrier space for the Mth electron group. To build up
the operators g_{M}^{+}, an orthogonal basis set { c_{m}} of vectors (functions) is introduced in each of the subspaces
L_{M}. If we define the operator m^{+} as one creating an
electron in the oneelectron state c_{m}, then the operators g_{M}^{+} are constructed as expansions:
g_{M}^{+} = 
å
{ m_{i}}

C_{{ mi}} 
N_{M} Õ
\substack i = 1 m_{i} Î { M}

m_{i}^{+} 
 (0) 
with the amplitudes C_{{ mi} } to be determined
variationally. (Each {m_{i}} represents a selection of n_{M} elements
of the local basis {c_{m}}.) Within this construct the hitherto
vaguely defined parts of the system" become the carrier
subspaces. The number of particles operator for each group or carrier
subspace can be defined as:

^ N

M

= 
å
m Î { M}

m^{+}m^{}, 
 (0) 
the average occupancies of the parts n_{M} become the integer numbers of
electrons residing in each of the carrier subspaces, and the dispersions Eq.
(3) vanish as each group fuction g_{M}^{+} is an eigenfunction of
the respective operator [^N]_{M} with the eigenvalue n_{M}.
While the overall dispersion Eq. (4) vanishes for a wave
function represented by a single GF product Eq. (5). This is not the
case for the exact molecular wave function, which may be always represented
as a linear combination of several functions with a GF structure, but with
an arbitrary distribution of electrons between the groups. This physically
corresponds to accounting the charge transfer between the subsystems. If the
selection of the groups (break up of the entire orbital space into carrier
subspaces and assigning integer numbers of electrons to them) was
appropriate, then the dispersion Eq. (4) for the exact wave function
will be small.
This concept can be viewed alternatively as a separation of the Hamiltonian [^H] into its ßignificant" part [^H]_{0}:
which commutes with all the number of particles operators introduced above:

é ë


^ N

M

, 
^ H

0

ù û

= 0, "M 
 (0) 
and its noncommuting part [^W] playing the rôle of an ïnteraction"
or "perturbation" defined with respect to [^H]_{0}. This perturbation [^W] must be weak. Provided the set of the operators [^N]_{M} is
optimized with respect to this criterion, the entire molecular system
naturally falls out into subsystems. Under this condition the magnitude of
the dispersions Eq. (4) for the exact wave function is
controlled by the same smallness parameters as the convergence of the
perturbation series with the perturbation [^W]. The matrix elements of [^W] are nonvanishing only betwen the approximate ground state Y
(or other eigenvectors of [^H]_{0}) and those other states Y^{¢} which differ from Y by the
the numbers of electrons in some subsystems. Then the distinguishability
means that
where DE_{CT} stands for the characteristic energy of a
charge transfer state.
Our following considerations are based on the trial wave function in the GF
form. This form of the wave function assures a simplified structure of the
one and twoelectron density matrices and, as a consequence, a simplified
expression for the molecular energy, which can be written as a sum of
intragroup contributions and Coulomb and exchange interactions between
electron groups. Selecting an appropriate form of the electronic wave
function for each of the groups allows to formalize the chemical
observations on the transferability of certain parts of molecules and to
establish transferability of certain types of electronic structure variables
(ESVs) as well. The simplest form of the wave function relevant to our
purpose of reconstructing a classical bonding picture in organic compounds
is based on the antisymmetrized product of twoelectron functions (geminals)
[,]. This is particularly relevant in the case of
usual saturated organic molecules, the electronic structure of which can be
described by twocenter, twoelectron bonds and lone electron pairs 
basically following Lewis. In the semiempirical context each carrier space
could be specified by a pair of strictly local hybrid orbitals (HOs)
centered on the atoms connected by the bond at hand. This form of the trial
wave function in combination with the semiempirical Hamiltonians allowed us
to derive a mechanistic picture for the molecular energy Tch62,Tch64,Tch73 and to show a path to the junction between QM and MM
parts in the hybrid QM/MM methods. It has been demonstrated []
that the matrix elements of one and twoelectron density matrices are
stable  and thus transferable for analogous bonds  up to certain degree
in two small parameters representing intrabond electron correlation and
polarity. Remarkably, the offdiagonal matrix element of the oneelectron
density matrix  the Coulson bond order  is stable (transferable) within
the second order in both these small parameters. It is noteworthy that the
transferability of the density matrix elements takes place namely in the
basis of the strictly local HOs. Although it is possible that the
transferability of the density matrix elements would take place in some
other basis as well, this is definitely not the case for an arbitrary basis (e.g., for the original AO basis).
The HOs used for constructing the geminals turn out to be very sensitive to
the variations of the molecular Hamiltonian irrespective to the origin of
these variations: either the geometry changes or the changes of the chemical
composition of the molecular system. The actual form and orientation of the
HOs spanning the carrier space for the respective bond geminal takes upon
itself a major part of the variations in the Hamiltonian, whereas the
density matrix elements in the basis of HOs remain relatively stable as
noted above.
In the present paper we start the long way towards an ab initio MM
theory. The semiempirical approach to deductive MM serves us as an example
and we try to extend it taking into account the complications arising in the
ab initio domain. Our first objective is to study the stability
(transferability) properties of the density matrix elements. The rest of the
paper addresses this issue: we devise a quantum chemical theory underlying
the tentative ab initio MM and use it to establish the
transferability of the density matrix elements on the solid theoretical
basis supported by numerical experiments.
Strictly local geminals formalism for the electronic wave function
The general and formal construction of the classical picture of chemical
bonding starts from selecting an adequate form of the trial wave function to
be used in the derivations and numerical experiments. The deductive MM
method based on the semiempirical Hamiltonians employs a geminalbased wave
function. An important restriction on the structure of the bonding geminals
is that each of them is constructed in a carrier space spanned by two
strictly local HOs  r_{m}
ñ and  l_{m}
ñ (four spinorbitals):
g_{m}^{+} = u_{m}r_{ma}^{+}r_{mb}^{+}+v_{m}l_{ma}^{+}l_{mb}^{+}+w_{m}( r_{ma}^{+}l_{mb}^{+}+l_{ma}^{+}r_{mb}^{+}) , 
 (0) 
where the amplitudes of configurations u_{m}, v_{m}, and w_{m} are
subject to the normalization condition:
u_{m}^{2}+v_{m}^{2}+2w_{m}^{2} = 1 
 (0) 
and should be determined variationally. In the case of a lone electron pair
the corresponding carrier space is spanned by one HO and the geminal is
trivial:
g_{m}^{+} = r_{ma}^{+}r_{mb}^{+}. 
 (0) 
This type of geminals based on the HOs can be considered as strictly local
geminals (SLG) [,]. The overall wave function of electrons in
the molecule is then taken as the antisymmetrized product of the geminals,
in accordance with Eq. (5).
The operators r_{ms}^{+} and l_{ms}^{+}
create an electron with the spin projection s in oneelectron states
centered on the "right" and on the "left" ends of the bond, respectively.
These operators are special cases of the general operators m_{i}^{+} of Eq. (8) for the case of twoelectron twoorbital
groups (geminals). The structure of these operators should reflect the
strictly local character of the geminals. These orbitals are obtained by
transformations involving only orbitals of the same atom:
t_{ms}^{+} = 
å
a Î A

h_{tma}(A)a_{s}^{+}; t = r,l. 
 (0) 
In this formula a_{s}^{+} stand for the operators creating
electrons in the (orthogonal) basis AOs centered at the atom A; the
matrices h(A) have the SO(4) structure and transform valence s and porbitals into a set of four orthogonal HOs. The independent parameters of
these matrices produce another subset of the ESVs, which describe the actual
hybridization related ones and should also be determined variationally.
This picture is natural for semiempirical settings, but it cannot be
directly transferred to the ab initio context where a series of
complications arises. First of all, the dimensionality of the atomic
subspace is larger than four due to explicit presence of the core shells for
the heavy atoms (second row and higher), even if a minimal basi set is used.
Another complication comes from the fact that the standard ab initio
methods employ nonorthogonal AOs basis sets and one cannot form orthogonal
HOs simply by performing intraatomic orthogonal transformations of them.
Nevertheless, we can propose a procedure for obtaining approprate
oneelectron states, which can be naturally assigned to geminals. We perform
this procedure in two steps. The first step is an orthogonalization of the
initial set of AOs (including orbitals assigned to different atoms).
Obviously there is some arbitrariness in this step. We employ the symmetric
(Löwdin) orthogonalization to obtain the orthonormal set of orbitals,
which is as close to the initial one as only possible. (It also conserves
the assignment of the orbitals to atoms.) The second step is an intraatomic
unitary transformation of orbitals. It is performed for each atom
separately:
t_{ms}^{+} = 
å
a Î A

H_{tma}(A)a_{s}^{+}, 
 (0) 
where the operators a_{s}^{+} correspond to the
orthonormalized set of AOs and the transformations H(A) is not limited to
four orbitals. The whole transformation is a product of the
(atomic)blockdiagonal transformation matrix with blocks given by the
matrices H(A) and the symmetrical orthonormalization matrix S^{1/2}. It
is noteworthy that the order of the two steps can be interchanged, i.e. the whole procedure can be viewed as a formation of nonorthogonal HOs
by the transformations H(A) with subsequent orthogonalization of all HOs.
The assignment of HOs obtained by the transformation Eq. (17) to
twoelectron groups is not yet specified. It is possible to distribute all
the orbitals between geminals or to use one orbital as a core and a subset
of four orbitals as valence HOs assigned to geminals. At the same time, our
ultimate goal of deriving a classical (MM) picture makes preferable to
recover hybridization structure Eq. (16) from the general
transformation H(A) because it allows to construct the most important
valence bond configurations and to relate the structure of oneelectron
states to the stereochemistry of an atom. We impose a set of restrictions on
the structure of the matrix H(A) and consider it as a product of
transformations. First we form new shells (sets of basis AOs; an sshell
is given by one function while a pshell is given by three functions with
the same radial part). The optimal sets of s and pfunctions are
obtained by two independent orthogonal transformations of the basis s and
pshells, respectively. Among these shells one sshell (AO) corresponds
to the core, while the remaining s and pshells can be combined to
produce some spshells (basis subsets containing one sAO and three pAOs with correct transformation properties). In our subsequent calculations
we use only one of these spshells. At the second step we apply some
hybridization SO(4) transformation mixing one s and three porbitals
for each spshell. Therefore, we produce the same setting for studying the
transferability of the density matrix elements as we had in the
semiempirical case. The structure of the HOs and the algorithms for
determination of the hybridizationrelated ESVs are important for the
construction of ab initio MM and these questions will be addressed
elsewhere. Here we concentrate on the density matrix elements in the basis
of these HOs.
The use of only the minimal number of basis orbitals at this stage of
development may be justified by the known fact that even the use of STONG type minimal basis sets at the RHF level of theory represents
an often surprisingly good ``model chemistry". Furthermore, the a
posteriori analyses of large scale calculations usually also lead to as
many significantly populated effective core and valence orbitals as is the
number of functions in the minimal basis sets [,].
Biorthogonal treatment for the geminals with overlap
General biorthogonal quasienergy formula for geminals
Overlap of the orbitals assigned to the different bonds makes it difficult
to calculate expectation values  in particular the total energy  already
at the single determinant level of theory,
However, the situation there is much simpler because any single determinant
wave function can be expressed (up to an insignificant normalization
constant) by using orthonormalized oneelectron orbitals, which is not the
case for geminals.
In the single determinant case one usually calculates the energy after
performing an orthogonalization of the orbitals. However, it is easy to show
that the energy of a single determinant wave function Y_{SD}ñ
can be given also as
E = 
áY_{SD}Y_{SD}ñ

= á 
~ Y

SD

 
^ H

Y_{SD}ñ 
 (0) 
where [(Y)\tilde]_{SD}ñ is the single determinant built up of the
orbitals {[(j)\tilde]_{i}} representing the biorthogonal
counterparts of the nonorthogonal orbitals {j_{i}} contained in
the determinant Y_{SD}ñ:

~ j

i

= 
å
j

(S^{1})_{ji}j_{j} 
 (0) 
Here s is the overlap matrix of the occupied orbitals
having the elements S_{ij} = áj_{i}j_{j}ñ.
(When writing equation (18), we have taken into account that á[(Y)\tilde]_{SD}Y_{SD}ñ = 1, owing to the biorthogonality of the
two sets of orbitals.)
No such strict relationship holds for the geminal wave functions. However,
at least at the equilibrium geometries, each geminal is expected to be not
too far from its respective HartreeFock approximation. Therefore, for
qualitative considerations it may be of meaning to consider the biorthogonal
quasienergy of the geminal wave function Y_{G}ñ:
We note that the equality [E\tilde] = E would hold also if Y_{G}ñ were the exact (full CI) wave function.
The idea of using the biorthogonal quasienergy for qualitative
considerations of geminals, describing different chemical bonds is not new;
for instance it has been proposed in the book of Surján [].
However, there one cannot find all the necessary quantities, and the
expressions are in a form not very suitable to be used in our case.
Therefore we have decided to present here the result explicitly.
It is easy to see, that if one introduces the biorthogonal transformation of
the basis orbitals l_{m}ñ,r_{m}ñ then the geminals, [g\tilde]_{m}ñ formed with the same coefficients u,v and w of
the latter will be ßtrongly biorthogonal" to the original ones (provided
that the orbitals within same geminal have been orthogonalized).
As a conseqence, the equality á[(Y)\tilde]_{G}Y_{G}ñ = 1
holds, oneelectron matrix elements appear only between orbitals assigned to
the same geminal and twoelectron ones involve at most two geminals. A
somewhat lengthy but essentially straightforward derivation gave the result:




= á 
~ Y

G

 
^ H

Y_{G}ñ\notag 
 (0)  


n å
i = 1


ì í
î

2(u_{i}^{2}+w_{i}^{2})á 
~ r

i

 
^ h

r_{i}ñ+2w_{i}(u_{i}+v_{i})(á 
~ r

i

 
^ h

l_{i}ñ+á 
~ l

i

 
^ h

r_{i}ñ) 
 (0)  

+2(v_{i}^{2}+w_{i}^{2}))á 
~ l

i

 
^ h

l_{i}ñ+u_{i}^{2}[ 
~ r

i


~ r

i

r_{i}r_{i}]+2u_{i}w_{i}[ 
~ r

i


~ r

i

l_{i}r_{i}]+u_{i}v_{i}[ 
~ r

i


~ r

i

l_{i}l_{i}] \notag 
 (0)  

+2u_{i}w_{i}[ 
~ r

i


~ r

i

r_{i}r_{i}]+2w_{i}^{2}[ 
~ r

_{i} 
~ l

i

l_{i}r_{i}]+2w_{i}^{2}[ 
~ r

i


~ l

_{i}r_{i}l_{i}]+2u_{i}w_{i}[ 
~ r

i


~ l

i

l_{i}l_{i}] \notag 
 (0)  

+ +u_{i}v_{i}[ 
~ l

i


~ l

i

r_{i}r_{i}]+2v_{i}w_{i}[ 
~ l

i


~ l

i

l_{i}r_{i}]+v_{i}^{2}[ 
~ l

i


~ l

_{i}l_{i}l_{i}] 
ü ý
þ

+2 
å
i < j


~ E

Cx ij

\notag 
 (0) 

Here [E\tilde]_{ij}^{Cx} describes the Coulomb and exchange interactions
between the ith and jth geminals in the biorthogonal framework; it is
given by the following expression


(u_{i}^{2}+w_{i}^{2})(u_{j}^{2}+w_{j}^{2}) 
æ è

2[ 
~ r

i


~ r

j

r_{i}r_{j}][ 
~ r

i


~ r

_{j}r_{j}r_{i}] 
ö ø

\notag 
 (0)  

w_{i}(u_{i}+v_{i})(u_{j}^{2}+w_{j}^{2}) 
æ è

2[ 
~ r

i


~ r

_{j}l_{i}r_{j}][ 
~ r

i


~ r

j

r_{j}l_{i}] 
ö ø

\notag 
 (0)  

(u_{i}^{2}+w_{i}^{2})w_{j}(u_{j}+v_{j}) 
æ è

2[ 
~ r

i


~ r

_{j}r_{i}l_{j}][ 
~ r

i


~ r

j

l_{j}r_{i}] 
ö ø

\notag 
 (0)  

w_{i}(u_{i}+v_{i})w_{j}(u_{j}+v_{j}) 
æ è

2[ 
~ r

i


~ r

_{j}l_{i}l_{j}][ 
~ r

i


~ r

j

l_{j}l_{i}] 
ö ø

\notag 
 (0)  

(u_{i}^{2}+w_{i}^{2})w_{j}(u_{j}+v_{j}) 
æ è

2[ 
~ r

i


~ l

_{j}r_{i}r_{j}][ 
~ r

i


~ l

j

r_{j}r_{i}] 
ö ø

\notag 
 (0)  

w_{i}(u_{i}+v_{i})w_{j}(u_{j}+v_{j}) 
æ è

2[ 
~ r

i


~ l

_{j}l_{i}r_{j}][ 
~ r

i


~ l

j

r_{j}l_{i}] 
ö ø

\notag 
 (0)  

(u_{i}^{2}+w_{i}^{2})(u_{j}^{2}+w_{j}^{2}) 
æ è

2[ 
~ r

i


~ l

_{j}r_{i}l_{j}][ 
~ r

i


~ l

j

l_{j}r_{i}] 
ö ø


 (0)  

w_{i}(u_{i}+v_{i})(v_{j}^{2}+w_{j}^{2}) 
æ è

2[ 
~ r

i


~ l

_{j}l_{i}l_{j}][ 
~ r

i


~ l

j

l_{j}l_{i}] 
ö ø

\notag 
 (0)  

w_{i}(u_{i}+v_{i})(u_{j}^{2}+w_{j}^{2}) 
æ è

2[ 
~ l

i


~ r

_{j}r_{i}r_{j}][ 
~ l

i


~ r

j

r_{j}r_{i}] 
ö ø

\notag 
 (0)  

(v_{i}^{2}+w_{i}^{2})(u_{j}^{2}+w_{j}^{2}) 
æ è

2[ 
~ l

i


~ r

_{j}l_{i}r_{j}][ 
~ l

i


~ r

j

r_{j}l_{i}] 
ö ø

\notag 
 (0)  

w_{i}(u_{i}+v_{i})w_{j}(u_{j}+v_{j}) 
æ è

2[ 
~ l

i


~ r

_{j}r_{i}l_{j}][ 
~ l

i


~ r

j

l_{j}r_{i}] 
ö ø

\notag 
 (0)  

(v_{i}^{2}+w_{i}^{2})w_{j}(u_{j}+v_{j}) 
æ è

2[ 
~ l

i


~ r

_{j}l_{i}l_{j}][ 
~ l

i


~ r

j

l_{j}l_{i}] 
ö ø

\notag 
 (0)  

w_{i}(u_{i}+v_{i})w_{j}(u_{j}+v_{j}) 
æ è

2[ 
~ l

i


~ l

_{j}r_{i}r_{j}][ 
~ l

i


~ l

j

r_{j}r_{i}] 
ö ø

\notag 
 (0)  

(v_{i}^{2}+w_{i}^{2})w_{j}(u_{j}+v_{j}) 
æ è

2[ 
~ l

i


~ l

_{j}l_{i}r_{j}][ 
~ l

i


~ l

j

r_{j}l_{i}] 
ö ø

\notag 
 (0)  

w_{i}(u_{i}+v_{i})(u_{j}^{2}+v_{j}^{2}) 
æ è

2[ 
~ l

i


~ l

_{j}r_{i}l_{j}][ 
~ l

i


~ l

j

l_{j}r_{i}] 
ö ø

\notag 
 (0)  

(v_{i}^{2}+w_{i}^{2})(v_{j}^{2}+w_{j}^{2}) 
æ è

2[ 
~ l

i


~ l

_{j}l_{i}l_{j}][ 
~ l

i


~ l

j

l_{j}l_{i}] 
ö ø

\notag 
 (0) 

One may introduce the biorthogonal Coulomb and exchangetype operators by
the definitions of their action on an arbitrary orbital j as

^

a_{j}b_{j}

j(1) = 
ó õ


r_{12}

dv_{2}j(1) \notag 
 (0)   (0)  
^

a_{j}b_{j}

j(1) = 
ó õ


r_{12}

dv_{2} b_{j}(1) \notag 
 (0) 

Using them, one may define the effective (biorthogonal) potential due to the
electrons in the jth geminal as


(u_{j}^{2} + w_{j}^{2}) 
æ ç
è

2 
^

_{rjrj}  
^

r_{j}r_{j}

ö ÷
ø

+(v_{j}^{2} + w_{j}^{2}) 
æ ç
è

2 
^

l_{j}l_{j}

 
^

l_{j}l_{j}

ö ÷
ø

\notag 
 (0)   (0)  

w_{j}(u_{j}+v_{j}) 
æ ç
è

2 
^

r_{j}l_{j}

 
^

r_{j}l_{j}

+2 
^

l_{j}r_{j}

 
^

l_{j}r_{j}

ö ÷
ø

\notag 
 (0) 

and then the effective (biorthogonal) core for the ith geminal will be

^

eff
i

= 
^ h

+ 
j ¹ i


^

eff
j

=  
1 2

D 
å
a


Z_{a} r_{a}

+ 
j ¹ i


^

eff
j


 (0) 

Then the biorthogonal quasienergy of the whole system becomes



n å
i = 1


ì í
î

2(u_{i}^{2}+w_{i}^{2})á 
~ r

i

 
^

eff
i

r_{i}ñ+ 2w_{i}(u_{i}+v_{i})(á 
~ r

i

 
^

^{eff}_{i}l_{i}ñ+á 
~ l

i

 
^

eff
i

r_{i}ñ) \notag 
 (0)  

+ 2(v_{i}^{2}+w_{i}^{2}))á 
~ l

i

 
^

eff
i

l_{i}ñ+ u_{i}^{2}[ 
~ r

i


~ r

i

r_{i}r_{i}]+ 2u_{i}w_{i}[ 
~ r

i


~ r

i

l_{i}r_{i}]+u_{i}v_{i}[ 
~ r

i


~ r

i

l_{i}l_{i}] \notag 
 (0)  

+ 2u_{i}w_{i}[ 
~ r

i


~ r

i

r_{i}r_{i}]+2w_{i}^{2}[ 
~ r

i


~ l

i

l_{i}r_{i}]+2w_{i}^{2}[ 
~ r

i


~ l

i

r_{i}l_{i}]+ 2u_{i}w_{i}[ 
~ r

i


~ l

i

l_{i}l_{i}] \notag 
 (0)  

+ +u_{i}v_{i}[ 
~ l

i


~ l

i

r_{i}r_{i}]+2v_{i}w_{i}[ 
~ l

i


~ l

i

l_{i}r_{i}] + v_{i}^{2}[ 
~ l

i


~ l

i

l_{i}l_{i}] 
ü ý
þ


 (0) 

Considerations based on excluding effects of BSSEtype
In the present study we are interested in the factors influencing the
geminal coefficients and the resulting density matrix elements, rather than
in the actual value of the total energy. Assuming that the orbitals r_{i},
l_{i} belonging to the same geminal are orthogonalized and that they are
optimized by using a sufficiently large local basis, then one may use the
conventional (not biorthogonal) integrals when determining the
geminal coefficients .(According to our experience gained in the field of
intermolecular interactions [], the energy should
preferably be calculated without such a substitution.) That is the case
because the difference between an integral containing the biorthogonal
orbital in the bra" and that containing there the original
one can be connected with the finiteness of the basis  these effects cause
the known basis set superposition error" (BSSE) in the
calculation of intermolecular interactions. They are known to diminish as
the basis set improves.
For that reason we are going here to show that the difference between
integrals containing the biorthogonal and original orbitals is a basis
finiteness effect. As the derivation is not particularly related to the
geminal structure of wave functions, we shall use simplified notations m,n, [(n)\tilde] to denote the individual orbitals encountered. Also, the
derivation will be pertinent to some ``units" A  which in our case can be a
geminal.
Let us consider the difference of ïntraunit" (m,n Î A) integrals
and introduce the projector [^P]^{A} on the subspace of the basis
orbitals assigned to the given unit. If the orbitals within each
unit are orthonormalized, then
and we may write


á 
~ n

 
^ h

mñán 
^ h

mñ = \notag 
 (0)  

á 
~ n

( 
^ P

A

+1 
^ P

A

) 
^ h

mñán 
^ h

mñ\notag 
 (0)  

á 
~ n

 
^ P

A


^ h

mñán 
^ h

mñ+á 
~ n

(1 
^ P

A

) 
^ h

mñ 
 (0)  

á 
~ n

 
å
r Î A

rñár 
^ h

mñán 
^ h

mñ+á 
~ n

(1 
^ P

A

) 
^ h

mñ\notag 
 (0)  


å
r Î A

d_{nr}ár 
^ h

mñán 
^ h

mñ+á 
~ n

(1 
^ P

A

) 
^ h

mñ\notag 
 (0)  

á 
~ n

(1 
^ P

A

) 
^ h

mñ\notag 
 (0) 

which would vanish if the function [^h]mñ could exactly be
expanded with the aim of the basis orbitals assigned to group A  thus it
is connected with the BSSE effects. Analogous considerations apply for the
twoelectron integrals, too.
Stability of the density matrix elements
In the present Section we prove that under very nonrestrictive conditions
the matrix elements of the one and twoelectron density matrices satisfy
the stability (transferability) conditions also for the ab initio
effective bond Hamiltonians, which are defined by the chemical composition
and geometry of the molecular system. We show that for a broad class of
molecules the latter two characteristics do not affect to a significant
degree the values of the density matrices in the basis of the symmetrically
orthogonalized HOs.
Effective bond Hamiltonians
Within the original SLG approach [,] the geminals are
characterized by the amplitudes u_{m}, v_{m}, and w_{m} (see Eq. (13)). We introduce a new notation z_{m} = Ö2w_{m} simplifying
the normalization condition for the amplitudes to: u_{m}^{2}+v_{m}^{2}+z_{m}^{2} = 1. The effective Hamiltonians for the
geminals representing chemical bonds are expressed in terms of the molecular
integrals in the HO basis. Each geminal is expanded over three singlet
twoelectron basis configurations and the optimal values of their amplitudes
are given by the solutions of the eigenvector problems (see also ParksParr):

æ ç ç
ç è



ö ÷ ÷
÷ ø


æ ç ç
ç è



ö ÷ ÷
÷ ø

= e_{m} 
æ ç ç
ç è



ö ÷ ÷
÷ ø

, 
 (0) 
corresponding to their lowest eigenvalues. The actual expressions for the
matrix elements of the effective bond Hamiltonians defined in the basis of
orthogonal (and real) HOs are (subscript m is omitted for brevity):

 
 

= (h_{rr}+h_{ll})+[rlrl]+[llrr], 

 
 
 


 (0) 
where h_{tt¢} stands for the matrix element of the oneelectron
part of the effective Hamiltonian for the bond under consideration and [tt^{¢}t^{¢¢}t^{¢¢¢}] are the
twoelectron Coulomb repulsion integrals written according to the [1212]
convention. The difference with the semiempirical case is that F ¹ 0 and D ¹ G, in general.
It is constructive to perform the transformation to the basis of normalized
symmetric and antisymmetric combinations of the ionic configurations [1/(Ö2)](r_{ma}^{+}r_{mb}^{+}± l_{ma}^{+}l_{mb}^{+}). Then the
effective Hamiltonian can be written as:
where the first state in the list is the symmetric combination of the ionic
configurations, the second state is the covalent configuration, and the last
one is the antisymmetric combination of the ionic configurations. This form
of the effective Hamiltonian allows to extract the most important composite
parameters determining the structure of the bond. A key parameter is the
difference of the energies of the symmetric ionic and covalent states:
Dg = 
1 2

(R+L)+FC = 
1 2

( [rrrr]+[llll]) [rlrl] > 0. 
 (0) 
This parameter equals to the difference between the average of the
intrahybrid Coulomb repulsions and their interhybrid counterpart. In the
simple symmetric case it would correspond to the difference of the oneand
twocenter Coulomb integrals g_{11} and g_{12}. Another
important parameter is b:
b =  
1 2Ö2

(D+G) = h_{lr} 
1 2

([rlrr]+[rlll]) > 0, 

describing electron hopping between the ends of the bond. The other two
characteristics represent the bond's asymmetry:


= 2(h_{rr}h_{ll})+[rrrr][llll], 

 


 (0) 
The sign of these matrix elements in the case of polar bonds depends on the
arbitrary assignment of the "right" and "left" labels to atoms. If the
"right" atom is more electronegative than the "left" one then RL < 0 and, typically, DG > 0. We use this
convention in the following.
It is convenient to use some dimensionless parameters instead of the matrix
elements. We use Dg as a measure of all the interactions and
introduce a set of dimensionless parameters:
All the parameters introduced are nonnegative and the last two parameters
do not appear in the semiempirical case. They allow to rewrite the
effective Hamiltonian Eq. (63) in a more convenient form:

^ H

= C 
^ I

+Dg 
æ ç ç
ç è


ö ÷ ÷
÷ ø

. 
 (0) 
It is remarkable that the matrix in the parentheses determines the structure
of the geminal irrespective to the average value of energy. In the important
case of symmetric bonds (or symmetrized effective Hamiltonians) all the
variety of possible solutions (ground state wave functions) is parameterized
by a single parameter z depending on the given chemical composition
of the bond, the shape and orientation of the HOs and the bond length.
Already this observation brings one to the conclusion that some degree of
the transferability of geminal amplitudes is can be expected, because the
multidimensional manifold of the matrix elements of the molecular electronic
Hamiltonians boils down to a single parameter z, which by only this
reason can be the same for very different molecules. Below we, however,
shall show even more: namely, that the significant components of the
electronic structure description in the SLG approximation  the elements of
the one and twoelectron density matrices in the HOs basis are described by
the formulae which further reduce the possible variance of the the
quantities of interest. This result may look trivial because the structure
of the symmetric geminal actually depends on only one parameter but the
point is that the parameter z has a simple and explicit expression
through the parameters of the problem.
Ground state of the effective symmetrized bond Hamiltonian and
the stability of the density
The actual solution of the eigenvector problem for the effective symmetrized
Hamiltonian (with the antisymmetric ionic state excluded) reduces to
diagonalizing 2×2 matrix
depending on the parameter z (see. Eq. (70)). We denote
the amplitudes of the symmetric ionic and covalent configurations in the
ground state as sinj and cosj, respectively.
Diagonalizing the Hamiltonian Eq. (70) yields (in the symmetric
case) the following solution for the ground state:
or in a detailed form,



sinj Ö2

= 
1 2


æ ú
Ö


, \notag 
 (0)  

 (0)  


cosj Ö2

= 
1 2


æ ú
Ö


, \notag 
 (0) 

where the notation
is introduced for the hereinafter ubiquitous square root.
The amplitudes of configurations enter the expressions for any observables
(including the energy) not by themselves, but only in some combinations:
namely, the elements of the one and twoelectron density matrices written
in the basis of the HOs spanning the carrier space for the bond geminal
under consideration. Those for the mth bond (geminal) are by definition:


= 
å
s

á 0 g_{m}t_{ms}^{+}t_{ms}^{¢}g_{m}^{+}0
ñ , 


= 
å
s

á 0 g_{m}t_{ms}^{+}t_{ms}^{¢+}t_{ms}^{¢¢}t_{ms}^{¢¢¢}g_{m}^{+}0
ñ . 



 (0) 
All the necessary matrix elements can be readily determined as simple
functions of geminal amplitudes u_{m}, v_{m} and w_{m} taking into
account the actual definition of the geminals Eq. (13). This
results in the following:

P_{m}^{rr} = 2( u_{m}^{2}+w_{m}^{2}) , P_{m}^{ll} = 2(v_{m}^{2}+w_{m}^{2}) , P_{m}^{rl} = P_{m}^{lr} = 2(u_{m}+v_{m})w_{m}, 

G_{m}^{rrrr} = 2u_{m}^{2}, G_{m}^{llll} = 2v_{m}^{2}, G_{m}^{rrll} = G_{m}^{lrrl} = 2w_{m}^{2}, 

G_{m}^{lrlr} = G_{m}^{rlrl} = 2u_{m}v_{m},G_{m}^{rrlr} = G_{m}^{rlrr} = 2u_{m}w_{m},G_{m}^{rlll} = G_{m}^{lrll} = 2w_{m}v_{m} 



 (0) 
where we took into account that some of the integrals appear twice in the
expression for the energy for different spin projection indices on which the
density matrix elements are themselves independent. Only the above
quantities each of which is not more complicated than involving one square
root irrationality (for the symmetric or symmetrized bond) enter the
expression for the energy in terms of the Hamltonian matrix elements in the
HOs basis. In the following we shall consider only three density matrix
elements P_{m}^{rl} = 2w_{m}(u_{m}+v_{m}), P_{m}^{rr}P_{m}^{ll} = 2(u_{m}^{2}v_{m}^{2}), G_{m}^{rrll} = 2w_{m}^{2},which naturally reproduce three important
characteristics: Coulsontype bond order (offdiagonal element of the first
order density matrix), bond polarity and bond covalency.
The solution of the symmetric problem results in the following values of the
density matrix elements:
P^{rl} = 
z G(z)

, P^{rr}P^{ll} = 0, G^{rrll} = 
1 2


æ ç
è

1+ 
1 G(z)

ö ÷
ø

. 
 (0) 
As we mentioned above, the very possibility to characterize any symmetrical
bond by a single parameter explains a lot of the observed transferability:
the entire diversity of semiempirical parameters and/or ab initio
basis sets boils down to a single number (Anderson's parameter z).
The definition of z (Eq. (66)) shows that z®0 in the limit of infinite separation, corresponding to the homolitic bond
cleavage. However, the analysis performed in a semiempirical setting Ref.
[] indicated that the interatomic separations characteristic for
the real bonds correspond to values z >> 1. It means that the density
matrix elements have the asymptotic behavior described by:
P^{rl}\asymp 1 
1 2z^{2}

, G^{rrll}\asymp 
1 2


æ ç
è

1+ 
1 z

ö ÷
ø

, 
 (0) 
where z^{1} serves as a small parameter. The most important feature
of the above expression is that the variations of the offdiagonal matrix
element of the oneelectron density matrix (Coulsontype bond order) with
respect to the limiting value of unity are of the second order with respect
to the small parameter z^{1}. This explains the results of the
numerical experiments of Refs. [,] performed on organic
compounds of different classes (alkanes, alcohols, amines etc.) at
the semiempirical level, where the Coulson bond orders P_{m}^{rl} all acquired values between 0.92 and 1.0. (Of course, the asymmetry of bonds
had also been properly taken into account.)
The same behaviour is observed in the ab initio case, too. We are
going to support the present theoretical conclusions by the numerical
results performed in the framework of the fully variational SLG model
described in the previous Section employing the 631G basis set. We repeated
some calculations in larger basis sets 6311G and 6311++G and those
calculations demonstrated that although the numerical values were affected,
all the qualitative conclusions remain valid. Table presents
the results of our calculations for a set of simple test molecules with
characteristic chemical bonds in different environment. First, the Table
presents the numerical values of the parameters determining the structure of
the Hamiltonian Eq. (70). We can see that for all molecules at
their equilibrium geometry structures the values of z are
considerably larger than unity and that the second important is the
asymmetry parameter l (see below). The Table also presents the
estimates of the density matrix elements Eq. (78), which are
compared with their actual values obtained in the ab initio SLG
computations. As expected, the two schemes give exactly the same values for
nonpolar chemical bonds. The agreement is reasonable even for highly polar
bonds. At the same time the difference between the two schemes quickly
increases when the polarity of the bond becomes larger. It means that
further corrections taking into account the bond asymmetry are necessary.
But before we perform a more refined analysis of the density matrix
elements, we consider another model based on the symmetric structure of the
geminal.
Estimates of the density matrix elements from the HeitlerLondon
solution
The nonorthogonality of the oneelectron basis functions represents one of
the most important differences between the ab initio and
semiempirical settings compared in the present paper. In the semiempirical
setting as a rule  though nowadays it is not mandatory  the AOs centered
on different atoms are assumed orthogonal with all effects of
nonorthogonality reloaded to the parameters of the method. On the contrary,
in the ab initio setting the overlap of the different oneelectron
functions is treated explicitly. We have considered above the model based on
the single parameter z_{m}, reflecting the electron correlation
within a symmetrized geminal built up of Löwdinorthogonalized orbitals.
It led to establishing the transferability
properties and the numerical results similar to those obtained in the
semiempircal case []. However, in the ab initio setting
the overlap integral S of the nonorthogonalized HOs ascribed to bond
geminal can be considered as an alternative parameter, the effect of which
on the result may be important. These quantities were not explicitly used in
the above estimates since those were performed in the basis of orthogonal
twoelectron configurations constructed from the orthogonalized HOs. By
virtue of this, the overlap integral affected the matrix elements of the
effective bond Hamiltonian and the parameter z only indirectly.
Nevertheless, we can determine the overlap between the HOs with the
knowledge of the transformation matrices H(A) Eq. (17) and use it
explicitly for constructing geminal structure models.
It is known that the HeitlerLondon (HL) wave function gives a good estimate
for the bondlength and for the energy of the H_{2} molecule. It has the
form of the covalent (homopolar) contribution to the bond geminal, but it is
built up from the original, i.e. nonorthogonal HOs.
Transforming the HL function to the basis of orthogonal HOs results in the
appearance of the ionic terms relative to the orthogonal basis. Let us
assume that the states  r_{m}
ñ and  l_{m}
ñ employed above are the
symmetrically orthonormalized original states  r_{m}^{0}
ñ and  l_{m}^{0}
ñ so
that:

æ ç
è



ö ÷
ø

= S^{1/2} 
æ ç
è



ö ÷
ø

, 
 (0) 
where the 2×2 overlap matrix S contains 1 as diagonal
elements and the overlap S as offdiagonal ones. The original
nonorthonormal states are restored form the orthonormal ones by applying
matrix S^{1/2}. The normalized HL state:

1

 r_{ma}^{0+}l_{mb}^{0+}+l_{ma}^{0+}r_{mb}^{0+}
ñ 
 (0) 
can be readily transformed to the form Eq. (13) taking into
account that this transformation is a direct product of two S^{1/2} transformations for the orbitals. It means that the amplitudes are
determned by the elements of the matrix S:
 r_{ma}^{0+}l_{mb}^{0+}+l_{ma}^{0+}r_{mb}^{0+}
ñ = Sr_{ma}^{+}r_{mb}^{+}+r_{ma}^{+}l_{mb}^{+}+l_{ma}^{+}r_{mb}^{+}+Sl_{ma}^{+}l_{mb}^{+}, 
 (0) 
where both sides have the normsquare 2( 1+S^{2}) .
The symmetric form of the geminal representation allows to extract the
effective value z_{S} of the Anderson's parameter z parameter
as a function of the intrabond overlap, which is equal to:
This parameter determines all the geminal amplitudes and the density
matrices according to the formulae given above. To test this model we
performed a series of calculations similar to those described in the
previous subsection. The numerical results are presented in Table . This Table demonstrates that the values of z determined on
the basis of the overlap are significantly smaller than the values arising
from the effective bond Hamiltoniani, although they often have similar trend
when going from one bond to another. At the same time, the numerical results
within the HL model are relatively poor and this model can be used for the
construction of the mechanistic picture for the potential energy surfaces
only if some (uniform) scaling of the HL based estimates of density matrix
elements is performed.
Bond asymmetry
So far our consideration was limited to the properties of either symmetric
or symmetrized effective bond Hamiltonians. At the same time, the difference
between the estimates obtained in this manner and the actual values of the
density matrices can be significant (although not exceedeing 5% for the
Coulson bond orders and 20% for the twoelectron density matrix elements)
in the case of polar bonds. Qualitatively the symmetrized model obviously
does not reproduce at all the difference between the diagonal oneelectron
density matrix elements P^{rr}P^{ll}. In the semiempirical context Tch73 the bond asymmetry has been taken into account perturbatively. It
was natural, since a single unaccounted parameter l measured the
interaction between the symmetric and antisymmetric states and the entire
effect of the bond antisymmetry could be loaded upon it. Actually, instead
of the parameter l another dimensionless parameter m had been
introduced []:
This parameter is small and it naturally arised from the firstorder
correction to the bond polarity:
P^{rr}P^{ll} = 2m 
G(z)1 G(z)+1

, 
 (0) 
where the factor 2 is due to the spin summation in Eq. (76) and
the rest of the multiplier at m is smaller than 1 and tends to 1
for large z. It was possible to single out two types of contributions
to the parameter m. It can be broken down into a sum of hybridization
dependent component m_{0} corresponding to the bond itself and the rest
m_{1} describing the environment of the bond:
where the first part prevails by orders of magnitude.
In the ab initio context the situation is somewhat more complex.
First of all, the interaction of the symmetric ground state with the excited
antisymmetric one leading to the asymmetry of the bond does not reduce to
the difference of the diagonal matrix elements of the original effective
Hamiltonian (reflected by l), but contains also a contribution
coming from the asymmetry of the offdiagonal Coulomb interaction matrix
elements (expressed by c). Next, the energy denominators entering the
perturbation theory expansions are modified through the exchange matrix
element. Nevertheless, it is possible to generalize the principal stability
(transferability) results obtained previously in the semiempirical context
also to the ab initio setting. If we number the ground, symmetric
excited and antisymmetric ionic states of the symmetrized bond hamiltonian
as states 0, 1 and 2, respectively, then the energies and the matrix
elements of the perturbation can be written (omitting the constant C in the Hamiltonian Eq. (70) and expressing
everything in units of Dg):



1G(z) 2

, e_{1} = 
1+G(z) 2

, e_{2} = 1 
w 2

, \notag 
 (0)  

 (0)  


ccosjlsinj 2

, V_{12} =  
lcosj+csinj 2

. \notag 
 (0) 

Using these notations it is easy to write the perturbation theory
corrections to the wave function and the density matrix elements. The
firstorder correction to the bond polarity becomes:
P^{rr}P^{ll} = 2(mm^{¢}) 
G(z)1 G(z)+1w

, 
 (0) 
where the contribution m^{¢} comes from c and is defined
as:
The actual values of the parameters m and m^{¢} are given in
Table . The parameter m^{¢} is usually much smaller
than m (the only exception is the molecule LiH) where they are of the
same order.
The analysis in the semiempirical case [] demonstrated that the
offdiagonal density matrix elements P^{rl} and G^{rrll} deviate
from their ideal" values 1 and [1/2],
respectively, in the second order with respect to the bond asymmetry
perturbation. A similar result can be obtained in the ab initio
casei, too. However, it is not possible to consider (mm^{¢})
as a single parameter of the bond asymmetry because l and c
enter the matrix elements V_{02} and V_{12} in different combinations.
Nevertheless, the corrections remain quadratic functions with respect to
them (or, equivalently, with respect to m and m^{¢}),
although the expressions become somewhat cumbersome. This result is
important because it shows the degree of the transferability for different
density matrix elements. Moreover, taking into account that the corrections
due to c and w are typically small (and of opposite signs),
the equations derived in the semiempirical domain can be good approximations
to the problem. We performed the calculations of the density matrix elements
employing the secondorder corrections with respect to m and m^{¢} (In fact, it actually represented the firstorder correction as
far as bond polarity is conserned.) Our calculations (Table ) show that this approach produces
surprizingly accurate results even for strongly polar bonds, where the
asymmetrical perturbation is not small.
Nonperturbational estimate of stability (transferability)
The previous analysis of stability (transferability) with respect to the bond
polarity is based on the perturbation estimates. At the same time it is
possible to propose an alternative approach to the problem of establishing
transferability for the density matrix elements. Let us consider the
following parametric form for the geminal amplitudes:
Our previous computations [] show that the values of f and
r should be close to p/4 for slightly polar and slightly
correlated chemical bonds. Therefore, we use a natural transformation:
which leads to a new definition of the geminal amplitudes:



1 2

(cosqsinq)(cosh+sinh), \notag 
 (0)  


1 2

(cosqsinq)(coshsinh), 
 (0)  

 (0) 

where the parameters q and h are expected to be small. These
parameters allow to look on the transferability at a different angle.
Namely, the density matrix elements can be written as:


cos2qcosh\asymp 12q^{2} 
1 2

h^{2},\notag 
 (0)  

(1sin2q)sin2h\asymp 2h4qh = 2h(14q), 
 (0)  


1 2

(1+sin2q)\asymp 
1 2

+q. \notag 
 (0) 

These equations show the orders of the transferability of the density matrix
elements with respect to the presumably small parameters q and h. These relations remind the previously established ones using the
parameters z^{1} and m (for example, the bond order is
transferable up to second order with respect to both parameters).
The transferability properties of the density matrix elements are
demonstrated by the above expressions provided the actual values of the
parameters q and h are indeed small. Thus, it is necessary to
have some reasonable estimates of these characteristics through the
parameters of the problem (effective Hamiltonian). To this end we write the
energy as an approximate function of q and h omitting all the
contributions of the order higher than second with respect to these
parameters. This energy function can be explicitly minimized with respect to
its variables. As a result, we obtain the following expressions of the
sought parameters:



zl^{2}+lcw 2(z^{2}wzl^{2})

, 
 (0)  

 (0)  


zlzcl z^{2}wzl^{2}

. \notag 
 (0) 

We performed the actual calculations of the density matrix elements based on
these estimates. The results are presented in Table . The
estimates are good in the case of bonds with small polarity (like CH) but
the quality of the estimates deteriorates in the case of polar bonds. It
means that the estimates made by making use the perturbation theory are
better for the practical calculations. Of course, the somewhat worse results
obtained in this case characterize rather the quality of the estimates Eq. (104) than the transferability properties of the density matrix
elements.
Discussion
In the present paper we attempted to extend the transferability of density
matrix elements on organic compounds of different classes (alkanes,
alcohols, amines etc.) established in Refs. [,] in
the semiempirical context to the ab initio Hamiltonian. This can be
considered as a prerequisite for developing a mechanistric model of
molecular potential energy surface which in its turn is consistent with ab
initio model of molecular electronic structure. It turns out to be possible
provided tha wave function of the molecular system can be presented as an
antisymmetrized product of strictlylocal geminals. Under this condition we
have demonstrated a remarkable stability of all the geminal related ESVs.
The values of the polarity P_{m}^{rr}P_{m}^{ll} do not exceed 0.07 by
absolute value for the compounds containing carbon, nitrogen, and hydrogen
atoms. Also the ionicity u_{m}^{2}+v_{m}^{2} for a rich variety of bonds
had a stable value about 0.4. These features although were not completely
unexpected, since the transferability of the parameters of the single bonds
in organic compounds is well known experimentally, required a theoretical
explanation. This explanation had been given by the observation that the
values of the density matrix elements have to be close to their transferable
values:
which depend neither on molecular geometry nor on the nature nor on
hybridization of the atoms bonded within a narrow range. The same as we see
now applies to the estimates of the similar quantities in the ab
initio context. The small corrections to the transferable values of the
density matrix elements can be expressed in terms of the effective bond
Hamiltonians characterizing the amount of deviation of correlation and of
asymmetry of the considered bond as compared to the transferable values. In
the semiempirical context the above observations allowed to develop a series
of approximate expressions for the molecular energy functional written for
the SLG trial wave function. The approximate formulae for the density matrix
elements suggest several options: one may neglect all the bond specific
corrections and by this to arrive to the picture where all the density
matrix elements in the basis of orthogonal HOs are fixed at their
transferable values. This is termed as the FA (fixed amplitudes where the
geminal amplitudes u,v,w for all geminals are meant). In this case the
energy becomes the function only of molecular geometry and of the rest of
the ESVs of the molecular system which describe the shape and the
orientation of the HOs spanning carrier subspaces for different geminals
(bonds). Alternatives to the the FA picture are numerous. All of them can be
uniformly termed as TA (tuned amplitudes) approximations. Only the simplest
ones were explored in Refs. [,,,]. They can be
characterized as those neglecting any nontransferability corrections of the
order higher than one in any of the small parameters. It is reasonable to
assume that similar moves can be applied in the ab initio context
thus paving the route of a sequential transition from a strict quantum
chemical description of molecular electronic structure of now accepted
ab initio level to its classical description in terms of bonding
and nonbonding force fields.
Conclusion
In the present paper we explored the transferability features of the
approximate SLG wave function in the ab initio context and showed
that the important elements of molecular electronic structure: elements of
the one and twoelectron density matrices are transferable from one
molecule to another to the degree controlled by small composition and
environment specific parameters. If the mentioned corrections are taken into
account the corresponding matrix elements can be calculated as transferable
functions of the molecular integrals pertinent to each of the bonds.
Acknowledgments
This work has been generously supported by the TÉT grant No. RUS12/2004
and the FASI grant No. 02.434.11.7098 as a joint RussianHungarian research
project, by the Russian Foundation for Basic Research (RFBR) grants Nos.
050790067, 070301128 and by the Hungarian Research Fund (OTKA) grant No. T43558.
A.M.T. acknowledges the financial support from a Marie Curie Fellowship (contract MIF1CT2005021198).
References
 []
 V.G. Dashevskii. Conformations of organic molecules.
1973 [in Russian].
 []
 U. Burkert, N. Allinger. Molecular Mechanics.
ACS Publications, 1982.
 []
 C.A. Coulson. Rev. Mod. Phys., 32 (1963)
170.
 []
 A.L. Tchougréeff, A.M. Tokmachev. Int. J. Quant.
Chem. 96 (2004) 175184.
 []
 A.L. Tchougréeff. J. Mol. Struct. THEOCHEM 630
(2003) 243263.
 []
 A.M. Tokmachev, A.L. Tchougréeff. J. Comp. Chem.
26 (2005) 491505.
 []
 K. Ruedenberg. Rev. Mod. Phys., 34 (1962)
326376.
 []
 I. Mayer. Int. J. Quant. Chem. 23 (1983)
341363.
 []
 V.A. Fock, M.G. Veselov, and M.I. Petrashen'.
J. Exp. Theor. Phys., 10 (1940) 723 [in Russian].
 []
 E. Kapuy. Acta Phys. Acad. Sci. Hung. 11 (1960) 409;
ibid. 12 (1960) 185.
 []
 A.M. Tokmachev, A.L. Tchougréeff. J. Comp. Chem.
22 (2001) 752764.
 []
 A.M. Tokmachev, A.L. Tchougréeff. J. Phys. Chem.
A 107 (2003) 358365.
 []
 I. Mayer. J. Phys. Chem. 100 (1996) 62496257
 []
 R. Ponec and D.L. Cooper Faraday Discuss. 135 (2007)
3142.
 []
 P. R. Surján. Second Quantized Approach to Quantum
Chemistry. SpringerVerlag, Berlin et al. 1989.
 []
 I. Mayer Int. J. Quant. Chem. 70, (1998) 4163.
 []
 J.M. Parks and R.G. Parr. J. Chem. Phys. 28
(1958) 335.
 []
 A.M. Tokmachev, A.L. Tchougréeff. Zh. Fiz. Khim.
73 (1999) 347  357.
 []
 A.L. Tchougréeff, A.M. Tokmachev. Int. J. Quant.
Chem. 100 (2004) 667676.
Table 0: Dimensionless parameters of the Hamiltonian Eq. ( 70) and the estimates of the density matrix elements Eq. ( 78)
Molecule  Bond  z  l  c  w  P^{rl}(SLG)
 P^{rl}(z)  G^{rrll}(SLG)  G^{rrll}(z) 
H_{2}  HH  4.3609  0.0000  0.0000  0.2841  0.9747  0.9747  0.6117  0.6117 
CH_{4}  CH  6.4648  1.0072  0.0333  0.3691  0.9822  0.9882  0.5669
 0.5764 
NH_{3}  NH  6.7235  1.9218  0.1159  0.3148  0.9701  0.9891  0.5432
 0.5736 
H_{2}O  OH  6.5442  2.6164  0.1572  0.2712  0.9523  0.9885  0.5187
 0.5755 
HF  FH  6.1613  3.2140  0.2144  0.2351  0.9293  0.9871  0.4896  0.5801 
LiH  HLi  4.4647  1.5786  0.6056  0.3559  0.9727  0.9758  0.5854  0.6093 
CH_{3}F  CH  6.6071  1.1919  0.0854  0.3549  0.9814  0.9887  0.5626
 0.5748 
 FC  6.2957  3.4383  0.3229  0.2547  0.9265  0.9876  0.4814  0.5784

CH_{3}OH  OC  7.5618  2.6840  0.3101  0.3662  0.9634  0.9914  0.5209
 0.5656 
 OH  6.3887  2.5647  0.1463  0.2613  0.9519  0.9880  0.5203  0.5773

 CH  6.4450  1.1948  0.0611  0.3482  0.9802  0.9882  0.5635  0.5767

 CH  6.5107  1.0941  0.0797  0.3507  0.9822  0.9884  0.5654  0.5759

C_{2}H_{6}  CC  9.1266  0.0000  0.0000  0.5886  0.9941  0.9941  0.5545  0.5545 
 CH  6.4181  1.0056  0.0382  0.3492  0.9822  0.9881  0.5674  0.5770

F_{2}  FF  2.9903  0.0000  0.0000  0.1139  0.9484  0.9484  0.6586  0.6586 
H_{2}O_{2}  OO  4.7984  0.0000  0.0000  0.2297  0.9790  0.9790  0.6020  0.6020 
 OH  6.4239  2.7791  0.1510  0.2650  0.9457  0.9881  0.5111  0.5111

N_{2}H_{4}  NN  7.9519  0.0000  0.0000  0.4477  0.9922  0.9922  0.5624  0.5624 
 NH  6.6531  2.0949  0.1213  0.2995  0.9659  0.9889  0.5378  0.5743

 NH  6.7421  1.9392  0.1375  0.3011  0.9705  0.9892  0.5431  0.5734

Table 0: Intrabond overlap and the HeitlerLondon estimates of the density
matrix elements
Molecule  Bond  S  z(HL)  z  P^{rl}(SLG)  P^{rl}(HL)  G^{rrll}(SLG)  G^{rrll}(HL) 
H_{2}  HH  0.7035  2.7861  4.3609  0.9747  0.9412  0.6118  0.6689

CH_{4}  CH  0.7202  2.9931  6.4648  0.9822  0.9485  0.5669  0.6584

NH_{3}  NH  0.6922  2.6575  6.7235  0.9701  0.9359  0.5432  0.6761

H_{2}O  OH  0.6465  2.2213  6.5442  0.9523  0.9119  0.5187  0.7052

HF  FH  0.5934  1.8317  6.1613  0.9293  0.8777  0.4896  0.7396 
LiH  HLi  0.6691  2.4227  4.4647  0.9727  0.9244  0.5854  0.6908 
CH_{3}F  CH  0.7154  2.9312  6.6071  0.9815  0.9464  0.5626  0.6614

 FC  0.5533  1.5951  6.2957  0.9265  0.8473  0.4814  0.7656 
CH_{3}OH  OC  0.6142  1.9723  7.5618  0.9634  0.8919  0.5209  0.7261

 OH  0.6415  2.1803  6.3887  0.9519  0.9090  0.5203  0.7084 
 CH  0.7128  2.8975  6.4450  0.9802  0.9453  0.5635  0.6631 
 CH  0.7118  2.8857  6.5107  0.9822  0.9449  0.5654  0.6637 
C_{2}H_{6}  CC  0.6975  2.7169  9.1266  0.9941  0.9385  0.5545  0.6727 
 CH  0.7128  2.8982  6.4181  0.9822  0.9453  0.5674  0.6631 
F_{2}  FF  0.3728  0.8661  2.9903  0.9484  0.6547  0.6586  0.8780

H_{2}O_{2}  OO  0.4795  1.2454  4.7984  0.9790  0.7797  0.6020  0.8130 
 OH  0.6401  2.1690  6.4239  0.9457  0.9081  0.5111  0.7093 
N_{2}H_{4}  NN  0.6277  2.0719  7.9519  0.9922  0.9006  0.5624  0.7173 
 NH  0.6897  2.6304  6.6531  0.9659  0.9347  0.5378  0.6777 
 NH  0.6877  2.6092  6.7421  0.9705  0.9338  0.5431  0.6789 
Table 0: \protectm parameters and the estimates of the density matrix
elements taking into account the bond asymmetry
Molecule  Bond  m  m^{¢}  (P^{rr}P^{ll})(SLG)  (P^{rr}P^{ll})(m)  P^{rl}(SLG)  P^{rl}(m)  G^{rrll}(SLG)  G^{rrll}(m) 
H_{2}  HH  0.0000  0.0000  0.0000  0.0000  0.9747  0.9747  0.6118  0.6118 
CH_{4}  CH  0.1540  0.0059  0.2281  0.2287  0.9822  0.9822  0.5669
 0.5667 
NH_{3}  NH  0.2827  0.0198  0.4051  0.4075  0.9701  0.9702  0.5432
 0.5422 
H_{2}O  OH  0.3952  0.0277  0.5559  0.5622  0.9523  0.9528  0.5187
 0.5151 
HF  FH  0.5149  0.0404  0.6998  0.7100  0.9293  0.9312  0.4896  0.4805 
LiH  HLi  0.3450  0.1653  0.2605  0.2463  0.9727  0.9739  0.5854  0.5865 
CH_{3}F  CH  0.1784  0.0149  0.2532  0.2536  0.9815  0.9815  0.5626
 0.5625 
 FC  0.5394  0.0593  0.7206  0.7247  0.9265  0.9301  0.4814  0.4723

CH_{3}OH  OC  0.3519  0.0464  0.4896  0.4902  0.9634  0.9642  0.5209
 0.5191 
 OH  0.3966  0.0265  0.5554  0.5617  0.9519  0.9524  0.5203  0.5167

 CH  0.1832  0.0109  0.2645  0.2652  0.9802  0.9802  0.5635  0.5633

 CH  0.1661  0.0141  0.2344  0.2347  0.9822  0.9822  0.5654  0.5653

C_{2}H_{6}  CC  0.0000  0.0000  0.0000  0.0000  0.9941  0.9941  0.5545  0.5545 
 CH  0.1548  0.0069  0.2269  0.2275  0.9822  0.9822  0.5674  0.5673

F_{2}  FF  0.0000  0.0000  0.0000  0.0000  0.9484  0.9484  0.6586  0.6586 
H_{2}O_{2}  OO  0.0000  0.0000  0.0000  0.0000  0.9790  0.9790  0.6020  0.6020 
 OH  0.4275  0.0271  0.5999  0.6087  0.9457  0.9462  0.5111  0.5061

N_{2}H_{4}  NN  0.0000  0.0000  0.0000  0.0000  0.9922  0.9922  0.5624  0.5624 
 NH  0.3114  0.0209  0.4446  0.4479  0.9659  0.9661  0.5378  0.5364

 NH  0.2845  0.0234  0.4024  0.4042  0.9705  0.9707  0.5431  0.5421

Table 0: \protectq and \protecth parameters and the
corresponding estimates of the density matrix elements
Molecule  Bond  q  h  (P^{rr}P^{ll})(SLG)  (P^{rr}P^{ll})(qh)  P^{rl}(SLG)  P^{rl}(qh)  G^{rrll}(SLG)  G^{rrll}(qh) 
H_{2}  HH  0.1147  0.0000  0.0000  0.0000  0.9747  0.9738  0.6118  0.6137 
CH_{4}  CH  0.0666  0.1378  0.2281  0.2359  0.9822  0.9817  0.5669
 0.5664 
NH_{3}  NH  0.0373  0.2594  0.4051  0.4589  0.9701  0.9639  0.5432
 0.5373 
H_{2}O  OH  0.0024  0.3940  0.5559  0.7122  0.9523  0.9234  0.5187
 0.4976 
HF  FH  0.0709  0.5831  0.6998  1.0492  0.9293  0.8264  0.4896  0.4293 
LiH  HLi  0.0811  0.1745  0.2605  0.2867  0.9727  0.9719  0.5854  0.5808 
CH_{3}F  CH  0.0618  0.1534  0.2532  0.2647  0.9815  0.9807  0.5626
 0.5617 
 FC  0.0891  0.6171  0.7206  1.1112  0.9265  0.8026  0.4814  0.4114 
CH_{3}OH  OC  0.0087  0.3234  0.4896  0.5921  0.9634  0.9480  0.5209
 0.5087 
 OH  0.0012  0.3956  0.5554  0.7129  0.9519  0.9227  0.5203  0.4989 
 CH  0.0626  0.1614  0.2645  0.2776  0.9802  0.9793  0.5635  0.5625

 CH  0.0649  0.1416  0.2344  0.2433  0.9822  0.9817  0.5654  0.5647

C_{2}H_{6}  CC  0.0548  0.0000  0.0000  0.0000  0.9941  0.9940  0.5545  0.5547 
 CH  0.0672  0.1371  0.2269  0.2346  0.9822  0.9817  0.5674  0.5670

F_{2}  FF  0.1672  0.0000  0.0000  0.0000  0.9484  0.9446  0.6586  0.6641 
H_{2}O_{2}  OO  0.1042  0.0000  0.0000  0.0000  0.9790  0.9784  0.6020  0.6034 
 OH  0.0180  0.4429  0.5999  0.8023  0.9457  0.9029  0.5111  0.4820 
N_{2}H_{4}  NN  0.0629  0.0000  0.0000  0.0000  0.9922  0.9921  0.5624  0.5627 
 NH  0.0293  0.2913  0.4446  0.5180  0.9659  0.9562  0.5378  0.5293

 NH  0.0371  0.2574  0.4024  0.4557  0.9705  0.9644  0.5431  0.5371

Footnotes:
^{1}This contribution is dedicated to Prof. I.V. Abarenkov on occasion of his
75th birthday
File translated from
T_{E}X
by
T_{T}H,
version 2.67.
On 3 May 2007, 12:51.