Generic Molecular Mechanics as Based on Local Quantum Description of Molecular Electronic Structure

Content-type: text/HTML

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

Abstract

Molecular mechanics (MM) is nowadays widely used for modeling potential energy surfaces of large organic molecules. All the MM schemes use ad hoc forms for the potential functions and parameters specially adjusted to fit experimental or quantum chemical data. In this work we attempt to deduce a generic MM scheme starting from a local quantum mechanical description of molecular electronic structure. The basis for this derivation is the trial electronic wave function in the form of antisymmetrized product of strictly localized geminals. The MM scheme obtained does not require adjusting any parameters. The quality of numerical estimates obtained by this scheme is analyzed.

Introduction

The molecular mechanics (MM) [,] numerical schemes are widely used for modeling potential energy surfaces (PES) of large molecular systems. The practical usefulness of the MM schemes is obvious - they provide the equilibrium geometries, dynamic matrices and other features of the systems of large size when the quantum chemical procedures becomes inapplicable due to enormous computational costs. The MM results are sufficiently accurate for many purposes and for large objects MM schemes are the best considering the cost-efficiency ratio.

Despite the intense employment of the MM schemes they remain purely empirical and are only supported by the agreement of calculated and experimentally observed quantities. Such level of substantiation means that the success of the MM schemes in general must be considered as an experimental fact. This fact in its turn requires theoretical explanation which as any sequential description of molecular structure and properties must be based on quantum mechanics (QM). The question of the interrelations between the MM and QM approaches is well-known in the literature (see, for example, []). It is, however, not only of academic interest. Additional attention to the foundations of the MM schemes is caused by increasing interest to the hybrid QM/MM schemes [] where a chemically active part of the whole system is described by precise quantum chemical methods while the rest (more or less inert environment) is represented by an appropriate MM method. The hybrid QM/MM schemes use a combination of sometimes high quality QM and of an empirical MM procedure. The status of the combined approach thus remains unclear as is the form of the junction between the quantum and classical subsystems in the hybrid QM/MM schemes which is commonly chosen ad hoc since the relations between the QM and MM approaches are not clearly defined. To derive successively the explicit form of the junction it is necessary to assume also that of the wave function underlying the MM description. In this context the very concept of the MM must be clarified. The common feature of all MM schemes is that the contributions to the energy from the bonded and non-bonded atoms are treated differently, i.e., the bonding is somewhat externally assumed. Different MM schemes employ different specific forms of potentials and numerous systems of their parameters fitted to reproduce a variety of characteristics [,,,] for more or less wide class of molecules. Also during the evolution of the MM approach itself some increasingly sophisticated contributions are added to the original simple picture which allow to extend the MM approach to increasingly complex and not transparently tractable classes of the molecules where the metal containing complexes and other molecules with significant electron redistribution must be mentioned MMGK,CombaHambley,Kozelka. For the latter cases the modern MM schemes do not represent the energy as an explicit function of the geometry parameters but include some kind of optimization scheme (for example, in Ref. Kozelka the scheme adjusting the charges on atoms on the ground of some simple criterion is employed). We see that the definition of MM becomes even more diffuse. Summarizing, we can define the MM scheme as one representing the energy as a combination of bonding and non-bonding contributions which are either explicit functions of molecular geometry structure or can be obtained from the latter by simple non-iterative optimization procedure. In the present paper we suggest a derivation and a numerical test of a scheme of this type - a local parametric expression for the total energy of organic molecules on the basis of a quantum description which uses a special form of the trial wave function of electrons, namely that of the antisymmetrized product of strictly local geminals (APSLG).

The reason to choose the APSLG form of the underlying wave function is some intimate similarity in the energy expressions of the APSLG and MM schemes. In our previous works [,] we have derived the required form of the junction between the QM and MM subsystems in an assumption that the inert part of the combined system can be represented by the wave function in the APSLG form [,,]. At the same time the APSLG based QM method is rather complicated self-consistent iteration procedure with numerous diagonalizations of effective Hamiltonian matrices though of small dimensionality. Nevertheless, the APSLG approach seems to be appropriate for obtaining the MM-type schemes for different reasons: the APSLG approach provides the local description of the molecular electronic structure; it represents the energy as the sum of intrabond and interbond contributions and thus can be considered as a natural starting point for construction of different additive schemes. Since the MM scheme is itself a parametric empirical procedure, it is unnecessary to use sophisticated nonempirical description for deriving the latter. It is more logical to use a semiempirical method as a starting point. A semiempirical implementation of the APSLG trial wave function has been recently developed with use of the MINDO/3 form of the Hamiltonian and resulted in an O(N)-scaling method suitable for describing organic molecules. This method and its parameterization are described in details in Refs. [,,].

The present paper is organized as follows. In the next Section we first briefly summarize the main features of the APSLG approach. Next the one-step procedures for optimization of geminal amplitudes and hybridization matrices are discussed. In the Section 3 we propose a series of approaches for evaluation of the electronic energy of molecules in the APSLG approximation and evaluate their accuracy. Next we discuss the possibility of treating these approaches as generic MM-type descriptions. The last Section contains the summary of results.

Construction of generic MM schemes

Wave function and energy of underlying APSLG based QM method

The APSLG wave function of electrons in the molecule has the form:

| F ñ =
Õ
m
g_m⁺|0 ñ
(0)
and is an example of the general ''group function method'' [] wave function. Each geminal corresponds to a chemical bond or to a lone electron pair. In the case of a chemical bond the geminal (singlet two-electron function) is constructed of two orbitals | r ñ and | l ñ assigned to the right and left ends of the bond, respectively. In the second quantization formalism it can be written as:

g_m⁺ = u_mr_ma⁺r_mb⁺+v_ml_ma⁺l_mb⁺+w_m( r_ma⁺l_mb⁺+l_ma⁺r_mb⁺) .
(0)
First two terms are ionic configurations (two electrons on the right or on the left end of the bond, respectively), while the last term is the covalent (Heitler-London) configuration. The normalization condition for the geminal Eq. (2) reads:

u_m²+v_m²+2w_m² = 1.
(0)
In the case of lone pair only the first (ionic) term makes sense and the geminal can be written as:

g_m⁺ = r_ma⁺r_mb⁺.
(0)
The orbitals | r ñ and | l ñ are chosen as strictly localized hybrid orbitals (HO) [], i.e. they are obtained by applying 4×4 orthogonal hybridization matrices h^A Î SO( 4) to the set of atomic orbitals (AOs) attached to each heavy (non-hydrogen) atom A. The amplitudes of the configurations (u_m, v_m, w_m) for each bond and the matrices h^A for each heavy atom are to be determined [,,] with use of the variational principle.

The above function is applied to find an estimate for the system energy with the Hamiltonian of the MINDO/3 approximation. In the HO basis the latter can be rewritten as a sum of one- and two-center contributions:

H =
å
A
H_A+ 1
2

å
A ¹ B
H_AB
(0)
where

H_A = å
\Sb t_m Î A

s\endSb æ
è U_m^t-
å
B ¹ A
g_ABZ_B ö
ø t_ms⁺t_ms-

- å
\Sb t_m₁,t_m₂^¢ Î A

m₁ < m₂\endSb
å
s
b_{t_m₁t_m₂^¢}^AA( t_m₁s⁺t_m₂s^¢+h.c.) +

+ 1
2
å
\Sb t_m₁,t_m₂^¢ Î A

t_m₃^¢¢,t_m₄^¢¢¢ Î A\endSb
å
st
( \QATOPtm₁\QATOPt^¢m₂ | \QATOPt^¢¢m₃\QATOPt^¢¢¢m₄)^At_m₁s⁺t_m₃t^¢¢+t_m₄t^¢¢¢t_m₂s^¢,

(0)
and

H_AB = - å
\Sb t_m₁ Î A

t_m₂^¢ Î B\endSb b_{t_m₁t_m₂^¢}^AB
å
s
( t_m₁s⁺t_m₂s^¢+h.c.) +

+g_AB å
\Sb t_m₁ Î A

t_m₂^¢ Î B\endSb
å
st
t_m₁s⁺t_m₂t^¢+t_m₂t^¢t_m₁s,

(0)
where h.c. stands for the hermitean conjugation and t refers tp the r or l HO. Molecular integrals entering this Hamiltonian depend on molecular geometry and on the hybridization matrices h^A. The parameter describing attraction of electrons to a core reads:

U_m^t =
å
i Î A
(h_mi^A)²U_i(A).
(0)
The electron attraction to other cores and repulsion of electrons on the orbitals of different atoms g_AB are invariant under the basis transformations allowed. The one-center two-electron matrix elements of the Coulomb repulsion are:

( t_m₁t_m₂^¢ | t_m₃^¢¢t_m₄^¢¢¢) ^A =
å
i Î A
h_m₁i^Ah_m₂i^Ah_m₃i^Ah_m₄i^A(ii | ii)^A+

+
å
i < j Î A
(ii | jj)^A(h_m₁i^Ah_m₂i^Ah_m₃j^Ah_m₄j^A+

+h_m₁j^Ah_m₂j^Ah_m₃i^Ah_m₄i^A)+
å
i < j Î A
(ij | ij)^A×

×(h_m₁i^Ah_m₂j^Ah_m₃i^Ah_m₄j^A+h_m₁i^Ah_m₂j^Ah_m₃j^Ah_m₄i^A+

+h_m₁j^Ah_m₂i^Ah_m₃i^Ah_m₄j^A+h_m₁j^Ah_m₂i^Ah_m₃j^Ah_m₄i^A),

(0)
where i and j are the AOs. The resonance integral describing one-electron transfers between the HO t_m₁ (centered on the atom A) and the HO t_m₂^¢ (centered on the atom B) is expressed through the resonance integrals in the AO basis:

b_{t_m₁t_m₂^¢}^AB =
å
i Î A

å
j Î B
h_m₁i^Ah_m₂j^Bb_ij^AB.
(0)

The electronic energy is then the sum of four terms:

E^el = E₁+E_Coul⁽¹⁾+E_res+E_Coul⁽²⁾,
(0)
originating from (i) the electron-cores attraction (E₁), (ii) one-center electron-electron repulsion (E_Coul⁽¹⁾), (iii) resonance intrabond interaction (E_res), and (iv) two-center electron-electron repulsion (E_Coul⁽²⁾). These contributions to the energy have the form:

E₁ = 2
å
A

å
t_m Î A
(U_m^t-
å
B ¹ A
g_ABZ_B)P_m^tt,

E_coul⁽¹⁾ =
å
A

å
t_m Î A
(t_mt_m | t_mt_m)^AG_m^tt+

+2
å
A
å
\Sb t_m₁,t_m₂^¢ Î A

m₁ < m₂\endSb [ 2(t_m₁t_m₁ | t_m₂^¢t_m₂^¢)^A-(t_m₁t_m₂^¢ | t_m₂^¢t_m₁)^A]P_m₁^ttP_m₂^{t^¢t^¢},

E_res = -4
å
m
b_{r_ml_m}^R_mL_mP_m^rl,

E_coul⁽²⁾ = 4
å
A < B
g_AB å
\Sb t_m₁ Î A

t_m₂^¢ Î B

m₁ ¹ m₂\endSb P_m₁^ttP_m₂^{t^¢t^¢}+
å
m
g_{R_mL_m}G_m^rl,

(0)
where we use the notation R_m for the right- and L_m for the left-end atoms respectively of the m-th bond and introduce the intrabond matrix elements of one- and two-electron density matrices:

P_m^{tt^¢} = á 0| g_mt_ms⁺t_ms^¢g_m⁺| 0 ñ ,

P_m^rr = u_m²+w_m², P_m^ll = v_m²+w_m²,

P_m^rl = P_m^lr = (u_m+v_m)w_m,

G_m^{tt^¢} = á 0| g_mt_mb⁺t_ma^¢+t_ma^¢t_mbg_m⁺| 0 ñ ,

G_m^rr = u_m², G_m^ll = v_m², G_m^rl = G_m^lr = w_m².

(0)

The nuclear-nuclear (core-core) repulsion in the MINDO/3 approximation differs from the pure Coulomb repulsion. It has the form:

Eⁿⁿ = 1
2

å
a ¹ b
C_ab,
(0)
where a and b are the indices of the respective nuclei and the quantities C_ab are defined by the expression:

C_ab = Z_a Z_b (g_ab+D_ab)
(0)
with

D_ab
g_ab
® 0, when R_ab® ¥,

D_ab® ¥, when R_ab® 0.

(0)
The standard MM potentials [,,] usually implicitly include the Coulomb interaction of some charges residing on the atoms. The Mulliken atomic charges in the framework of the APSLG approximation are simply given by:

Q_A = 2
å
t_m Î A
P_m^tt-Z_A.
(0)
The total energy can be written as a sum of contributions with obvious classification:

E =
å
m
E_m+
å
A

å
k < m Î A
E_km^A+E_nonbond+E_Coul.
(0)
This expression contains the intrabond contributions E_m, the contributions from the interaction between the vicinal chemical bonds at one atom A (E_km^A), the core interaction between nonbonded atoms, and the Coulomb interaction of atomic charges. Eq. (18) thus resembles the MM-type form of the potential energy. The explicit form of the contributions is the following:

E_m =
å
t Î { r,l}
[2U_m^tP_m^tt+(t_mt_m | t_mt_m)^T_mG_m^tt]+

+2g_{R_mL_m}[G_m^rl-2P_m^rrP_m^ll]-4b_{r_ml_m}^R_mL_mP_m^rl;

E_km^A = d_{AT_k}d_{T_kT_m^¢}[ 4(t_kt_k | t_m^¢t_m^¢)^T_k-2(t_kt_m^¢ | t_kt_m^¢)^T_k] P_k^ttP_m^{t^¢t^¢};

E_nonbond = 1
2

å
a ¹ b
Z_a Z_bD_ab;

E_Coul = 1
2

å
a ¹ b
Q_a Q_b g_ab.

(0)
At the same time this resemblance is relative because the contributions Eq. (19) to Eq. (18) are not explicit functions on geometry parameters but depend on the bond geminal amplitudes u_m, v_m, w_m through the matrix elements of the density matrices and on matrices h^A through the molecular integrals. In the framework of the APSLG approach [] these quantities are determined by iterative optimization procedure with alternating diagonalizations and multi-dimensional direct minimizations. At the same time our previous calculations [] have shown the bond amplitudes to be very close for the identical bonds in the similar environment. The hybridization is also close for the atoms of definite type in the definite environment (for example, carbon atoms in ethane, butane and isobutane molecules have similar parameters of hybridization). The quantities u_m, v_m, and w_m and h^A determined variationally completely characterize the molecular electronic structure in the APSLG approximation. These quantities can be termed as electronic structure parameters (ESP). An approximate transferability of these quantities can be confirmed by numerical estimates performed in the present and in our previous studies on the semiempirical version of the APSLG approximation. This allows to construct parametric expressions for the energy. Two classes of approaches can be constructed. In the first one the zero approximation transferable ESPs can be taken and inserted into Eq. (19). In the second class of approaches the zero approximation ESPs are tuned by some one-step (non-iterative) procedure in order to accomodate the variations of the ESPs under geometry variations. In the next Section we consider some procedures to be used for non-iterative ESP tuning.

Non-iterative tuning of ESPs

Bond (geminal) amplitudes

The optimal values of u_m, v_m, and z_m are the solutions of the eigenvalue problem:

æ
ç
ç
ç
è

A_m
0
D_m

0
B_m
D_m

D_m
D_m
C_m
ö
÷
÷
÷
ø æ
ç
ç
ç
è

u_m

v_m

z_m
ö
÷
÷
÷
ø = E^m æ
ç
ç
ç
è

u_m

v_m

z_m
ö
÷
÷
÷
ø ,
(0)
where w_m is replaced by z_m = Ö2w_m, which simplifies the normalization condition for the geminals:

u_m²+v_m²+z_m² = 1,
(0)
corresponding to the lowest eigenvalue. The 3×3 matrix here is that of the effective two-electron Hamiltonian for the m-th bond. The quantities A_m, B_m, C_m, and D_m are:

A_m = 2(U_m^r-
å
B ¹ R_m
g_{R_mB}Z_B)+(r_mr_m | r_mr_m)^R_m+

+ å
\Sb t_m₁ Î R_m

m₁ ¹ m\endSb [ 4(r_mr_m | t_m₁t_m₁)^R_m-2(r_mt_m₁ | r_mt_m₁)^R_m] P_m₁^tt+

+4
å
B ¹ R_m
g_{R_mB} å
\Sb t_m₁ Î B

m₁ ¹ m\endSb P_m₁^tt;

B_m = 2(U_m^l-
å
B ¹ L_m
g_{L_mB}Z_B)+(l_ml_m | l_ml_m)^L_m+

+ å
\Sb t_m₁ Î L_m

m₁ ¹ m\endSb [ 4(l_ml_m | t_m₁t_m₁)^L_m-2(l_mt_m₁ | l_mt_m₁)^L_m] P_m₁^tt+

+4
å
B ¹ L_m
g_{L_mB} å
\Sb t_m₁ Î B

m₁ ¹ m\endSb P_m₁^tt;

C_m = 1
2
(A_m+B_m)+g_{R_mL_m}- 1
2

å
t_m = r_mlr_m
(t_mt_m | t_mt_m)^T_m;

D_m = -Ö2b_{r_ml_m}^R_mL_m.

(0)
As it has been mentioned above the geminal amplitudes are remarkably transferable in a series of similar compounds. The quality of this approximation can be further improved by adjusting the zero approximation values of the geminal amplitudes. In order to do so one can notice that all the matrix elements Eq. (22) of the effective Hamiltonian Eq. ( 20) depend only on the diagonal elements of one-electron density matrix P_m₁^rr and P_m₁^ll. Keeping in mind the normalization condition for the diagonal one-electron densities:

P_m^rr+P_m^ll = 1
(0)
we (following Ref. []) introduce for each single chemical bond the polarity parameter:

q_m = P_m^rr-P_m^ll = u_m²-v_m²
(0)
Then the diagonal densities can be expressed through these parameters:

P_m^rr = 1+q_m
2
, P_m^ll = 1-q_m
2
.
(0)
The calculation of A_m, B_m, and C_m for the m-th geminal requires the information on polarities of all other bonds (geminals). These polarities are defined by the transferable values of the zero approximation amplitudes u_m, v_m, and w_m for all the chemical bonds. Also intuitively transparent values of the effecitive atomic charges can be conveniently used. They also can be expressed through the bond polarities (see above):

Q_A =
å
t_m Î A
[ 1+(d_tr-d_tl)q_m] -Z_A.
(0)
Using this notation A_m and B_m can be rewritten:

A_m = 2U_m^r+(r_mr_m | r_mr_m)^R_m+2
å
B ¹ R_m
g_{R_mB}Q_B-2g_{R_mL_m}+

+ å
\Sb t_m₁ Î R_m

m₁ ¹ m\endSb [ 2(r_mr_m | t_m₁t_m₁)^R_m-(r_mt_m₁ | r_mt_m₁)^R_m] [ 1+(d_tr-d_tl)q_m₁];

B_m = 2U_m^l+(l_ml_m | l_ml_m)^L_m+2
å
B ¹ L_m
g_{L_mB}Q_B-2g_{R_mL_m}+

+ å
\Sb t_m₁ Î L_m

m₁ ¹ m\endSb [ 2(l_ml_m | t_m₁t_m₁)^L_m-(l_mt_m₁ | l_mt_m₁)^L_m] [ 1+(d_tr-d_tl)q_m₁].

(0)
Analytical formulae for the roots of the corresponding cubic equation are used in calculation. The lowest of the three roots E_m must be taken. Then the amplitudes of the geminal are given by the formulae:

u_m = z_m D_m
A_m-E_m
,

v_m = z_m D_m
B_m-E_m
,

z_m = 1

Ö

1+( [(D_m)/(A_m-E_m)]) ²+( [(D_m)/(B_m-E_m)]) ²

.

(0)
These expressions represent explicit formulae for the amplitudes u_m, v_m, and w_m.

Hybridization matrices

According to Ref. [] each hybridization matrix can be represented by the product of six Jacobi rotation matrices acting in the two-dimensional subspaces of the whole four-dimensional space. Therefore it depends on six angles w_sx^A, w_sy^A, w_sz^A, w_xy^A, w_xz^A, and w_yz^A. These six angles, however, play different roles: the first triple determines the hybridization itself and must be at least approximately transferable while three others rotate the set of the HOs attached to the atom A as a whole and can not possess any transferability. It is congenial to the modern MM schemes to consider the six angles determining hybridization matrices as quantities which are determined by minimizing a simple functional:

E^A =
å
t_m Î A
[2U_m^tP_m^tt+(t_mt_m | t_mt_m)^AG_m^tt-4b_{r_ml_m}^R_mL_mP_m^rl]+

+2 å
\Sb t_m₁,t_m₂^¢ Î A

m₁ < m₂\endSb [ 2(t_m₁t_m₁ | t_m₂^¢t_m₂^¢)^A-(t_m₁t_m₂^¢ | t_m₂^¢t_m₁)^A]P_m₁^ttP_m₂^{t^¢t^¢},

(0)
which is a sum of the intraatomic energy and of the resonance energy of bonds incident to the atom A. This approach is implemented in the present paper. In the next Section we consider a series of approximate schemes based on the energy expression Eq. (19) but employing different approximate procedures for the ESPs (geminal amplitudes and hybridization matrices).

Parametric forms for the energy

The important feature of the APSLG-MINDO/3 method [] is a remarkable transferability of its relevant ESPs - the geminal amplitudes and hybridization angles. It distinguishes the APSLG method from others (based, for example, on the SCF approximation where the corresponding ESPs - MO LCAO coefficients - are not transferable). The geminal amplitudes vary slightly when going from one molecule to another: for example, the quantities u_m, v_m, and w_m for the C-H bond in methane are 0.4897, 0.4178, and 0.5412, respectively, while in ethane the analogous quantities are 0.4805, 0.4272, and 0.5416. These amplitudes also remain approximately the same in the case of the C-H bonds in significantly more different environment (for example, in methylamin molecule these amplitudes are 0.4814, 0.4223, and 0.5431). That degree of transferability is even more accurate than one could guess.

At this point we are in a position to formulate a generic MM scheme with certain QM foundation. Indeed, if we assume the expression Eq. ( 19) for the energy to be an exact QM one, then schemes fitting into a diffuse concept of MM can be suggested. The key point with them is to accept specific procedures for defining the ESPs entering Eq. ( 19).

Within such a setting the parameters of the method can be subdivided in two groups: the first one comprises the parameters of the Hamiltonian of the underlying semiempirical QM approach (in our case those of the APSLG-MINDO/3 []). They do not require special adjusting; other parameters reflects the transferability of the electronic structure characteristics obtained by the APSLG approach and can be taken from calculations of simple molecules by the APSLG approach or specially adjusted.

Using such definition of the HO matrices we can define zero-order approximation to the energy Eq. (19) by the following procedure: the amplitudes u_m, v_m, w_m are chosen to be constant for each type of the bond and the angles w_sx^A, w_sy^A, w_sz^A, determining the hybridization are also constant for each type of the atom (type of the atom can be determined in the way analogous to that of MM: it depends on the closest environment). In other way, the values of the hybridization angles define the type of the atom in terms of its hybridization. The quantities w_xy^A, w_xz^A, w_yz^A are chosen on the basis of purely geometric correspondence between directions of HOs and those of the chemical bonds. The above assumption of perfect transferability of ESPs u_m, v_m, w_m and w_sx^A, w_sy^A, w_sz^A allows to represent the total energy as an explicit function of geometric structure. At the same time this approximation may turn out to be rather crude because the transferability of the zero approximation ESPs is good only for similar geometric structure and the large distortions of the molecule can impair the quality of this simple model. In fact the transferability is characteristic only for the interatomic separations close to the equilibrium bond lengths. If the bond becomes elongated the ionic contributions diminish while the amplitude (w_m) of the homeopolar configuration tends to 1/Ö2. In this case the transferability of the zero approximation amplitudes is broken and they must be readjusted. Analogously, if the valence angles are far from the equilibrium ones the degree of transferability of the hybridization angles is questionable as well. Incidentally, the quality of the approximation to the amplitudes u_m, v_m, w_m can be significantly improved without loss of the explicit character of the energy expression.

The general idea of constructing local parametric expressions based on the APSLG description of the molecular electronic structure can be verified by numerical estimates of the validity of different schemes for obtaining the ESPs. The results of previous Sections allow to propose a range of schemes for estimating the total energy by Eq. (19) differing by approximations used for obtaining the ESPs:

the geminal amplitudes and hybridization angles w are taken as perfectly transferable parameters and are not additionally tuned. This scheme can be termed as FAFO (fixed amplitudes and fixed orbitals);
the geminal amplitudes are additionally adjusted by using formulae Eq. (28) while the HOs remain fixed at their transferable values - the TAFO scheme (tuned amplitudes and fixed orbitals);
only the HOs are adjusted by minimizing the functionals Eq. ( 29) for each heavy atom but the geminal amplitudes are taken equal to their transferable values - the FATO scheme;
the geminal amplitudes are adjusted first and the HOs are corrected for the adjusted values of the bond amplitudes - the TATO scheme;
analogous to the previous scheme but the order of two adjustment procedures is changed - the TOTA scheme.

Results and Discussion

In the present paper we performed test calculations on some organic molecules with use of the above five schemes. As the first test we applied the above schemes to constructing the energy profile for the process of variation of one of the C-H bonds length in the methane molecule. The transferable hybridization matrix is taken to correspond to the sp³-hybridization and the zero approximation C-H bond amplitudes are taken from a calculation on the methane molecule at its experimental equilibrium geometry. When such a choice is made all the five schemes give in the equilibrium point the electronic structure (and other properties, for example, heat of formation) coinciding with those obtained by the original APSLG approach. Figure 1 represents the differences between heats of formation calculated by each of the five above MM-type schemes and by the APSLG method as a function of the C-H bond length. This Figure shows that all the five schemes work well in the vicinity of the equilibrium geometry. It can be concluded that both the minimum position and the heat of formation in the minimum remain the same as in the original QM APSLG approach JCC. Since the APSLG approach reproduces the experimental equilibrium geometry of methane with perfect accuracy the MM-type schemes do the same as well. At the same time the results of the zero-order approximation (FAFO) deviate from the precise APSLG one rather strongly for interatomic distances far from the equilibrium. This level of approximation can, however, satisfactorily cover the PES within a ±0.1 Å range near the equilibrium position. It is important to find the quantitative measure of deviation between PES calculated by different computational methods. In our case the difference between an MM-type PES and the APSLG PES must be characterized. As a convenient and representative measure in the one-dimensional case we choose the area (integral) between two energy profiles in the definite interval of variation of the geometric parameter. Table 1 contains such data for three different ranges of variation of the C-H bond length in methane for all five schemes. The methods TAFO and FATO somewhat improve the results of the FAFO scheme but the range of their validity is not strongly extended as compared to that of the FAFO scheme. The methods combining the optimization of HOs with adjustment of the geminal amplitudes work satisfactorily for all interatomic distances studied. At the same time the TOTA method yields the numerical values which are closer to those obtained by the underlying APSLG approach than the TATO scheme. The maximal deviation of the TOTA scheme from the exact APSLG one does not exceed 0.2 kcal mol^-1 in the whole considered range of geometry parameters.

The above example is not representative enough because in the methane molecule the sp³ hybridization can be a good approximation even for very large geometry distortions. The hybridization can be significantly perturbed when less symmetric molecule is considered. To show how these schemes work under larger deviation of the hybridization from one calculated by the APSLG at the equilibrium geometry, we consider as the next example distortions of the water molecule taken with the sp³ hybridization for the oxygen atom as the zero approximation. The amplitudes u_m, v_m, and w_m of the O-H bond were taken from the APSLG-MINDO/3 calculation of the water molecule: 0.5946, 0.3317, and 0.5179, respectively. The distortions considered were (i) the variation of the length of one O-H bond, (ii) simultaneous variation of the lengths of both O-H bonds, (iii) variation of the valence angle H-O-H (these three cover all the possible distortions of the water molecule). Figures 2-4 represent for these processes the deviations of the heats of formation calculated by MM-type schemes and the underlying QM APSLG method for these processes. Analogously, Table 1 contains the integral characteristics of these deviations.

These data show the importance of the quality of the zero approximation hybridization employed. The FAFO and TAFO methods result in very large deviations of the heats of formation from the APSLG method even in the case of geometries close to the equilibrium. Also in the cases of variation of the lengths of O-H bonds the minimum position is significantly displaced, and in the case of the variation of the valence angle the shift of the minimum position is very large. Much better results are obtained with use of the FATO scheme. This scheme in this case works even better than the TATO scheme. At the same time the best of these schemes is the TOTA scheme. It gives the results which are in perfect agreement with those of the underlying QM APSLG approach itself. The data of Figures 1-4 and Table 1 unambiguously demonstrate that only the TOTA scheme can be successfully and surely applied for description of the large portions of molecular PESs. Only this scheme is applied in our further analysis.

The molecules considered previously have only one type of bonds. Now we consider how the TOTA scheme works with different types of bonds present. Table 2 shows the results of calculations on the heats of formation for the ethane molecule as a function of the C-C interatomic separation. The zero approximation amplitudes u_m, v_m, and w_m for the C-H bond were taken from the APSLG calculation on the methane molecule in its equilibrium geometry and the analogous quantities for the C-C bond were taken from the calculation on the ethane molecule (0.4502, 0.4502, and 0.5453). The data in Table 2 show that the ethane molecule is well described by the TOTA scheme for all interatomic distances considered. Moreover, also for the large C-C distances the difference between the APSLG and the TOTA heats of formation becomes very small. The position of the minimum on the PES remains unchanged.

Table 3 contains the results of calculations on a series of organic molecules (for their experimental geometry parameters) by the QM APSLG-MINDO/3 method and by the TOTA scheme for the two sets of zero approximation ESPs (geminal amplitudes). The first type of ESPs contains no adjusted parameters. The geminal amplitudes for the C-H, C-C, O-H and C-O bonds are taken from the APSLG calculations on the methane, ethane, water, and methanol molecules at their experimental geometry parameters. The sp³ hybridization was taken as a zero approximation one for all heavy atoms. These results show that the TOTA scheme provides the heats of formation well corresponding to those of the APSLG approach. The deviation from the underlying APSLG scheme is significantly smaller than the precision of this scheme itself. It allows one to make a conclusion of a validity of the MM-type TOTA scheme for calculations of the heats of formation and equilibrium geometry structures of organic compounds with well defined chemical bonds and lone pairs. It can be noted that this parameterization takes zero approximation ESPs from calculations of the simplest representatives of organic molecules. At the same time these ESPs are somewhat different from those characteristic for more large molecules because the simplest molecules have somewhat specific structure in the class of homologues. Therefore, it can be reasonable to improve the results on the heats of formation by slight change of the parameters and thus, the second procedure is an attempt to improve the results of calculations by changing the initial amplitudes for C-H bond. They were taken as follows: 0.4739, 0.4306, and 0.5431. The change of the parameterization allows to improve the agreement between the heats of formation calculated by the APSLG and TOTA schemes. At the same type this improvement is not very significant and the initial geminal amplitudes can be left unadjusted.

Conclusion

The present work is aimed to constructing parametric expressions for the energy of molecules by making start from an underlying QM description of molecular electronic structure. Such a construct could serve as a generic form of MM. The APSLG trial wave function was proven to be the proper choice for underlying the desirable MM description. In such a way a range of different MM-type schemes has been constructed for description of the total molecular energy as a function of molecular geometry. These generic MM schemes can be characterized as the QM APSLG-MINDO/3 method [] accompanied by special procedures of selection of the relevant electronic structure parameters - geminal amplitudes and hybridization angles. Within such a setting the parameters of generic MM schemes are (i) those of the underlying QM Hamiltonian, and (ii) the electronic structure parameters (ESPs). The latter are calculated by the APSLG-MINDO/3 method for some simple molecules and are used further according to a procedure adopted. The deformations of methane and water molecules have been used to discriminate the quality of different schemes for selecting the ESPs. The TOTA scheme results in a reliable agreement with the underlying QM APSLG method. The performed calculations on the organic molecules of different classes confirm this conclusion. An attempt to change the initial electronic structure parameters leads to slight improvement of the results.

Acknowledgements

The authors gratefully acknowledge valuable discussions with Dr. I.V. Pletnev. This work is performed under partial financial support on the part of the Federal Target Program of Russia ''Integracia'' (grant No A0078) and on the part of the RAS grant program for younger researchers No 6 (Grant No 120). Financial support for AMT from the Haldor Tops oe A/S is acknowledged as well.

References

[]: Dashevskii V.G. Conformations of Organic Molecules, Khimiya, Moscow, 1974 [in Russian].
[]: Burkert U.; Allinger N.L. Molecular mechanics, ACS, Washington, DC, 1982.
[]: Bernardi F.; Olivucci M.; Robb M.A. J Amer Chem Soc 1992, 114, 1606.
[]: CECAM-NSF Meeting on QC/MM methods, Int J Quantum Chem 1996, 60, No 6, Special Issue.
[]: Allinger N.L. J Amer Chem Soc 1977, 99, 8127.
[]: AllingerN.L.; Yuh Y.H.; Lii J.-H. J Amer Chem Soc 1989, 111, 8551.
[]: Brooks B.R.; Bruccoleri R.E.; Olafson B.D.; States D.J.; Swaminathan S.; Karplus M. J Comp Chem, 1983, 4, 187.
[]: Jorgensen W.L.; Tirado-Rives J. J Amer Chem Soc 1988, 110, 1657.
[]: Pletnev I.V.; Melnikov V.L. Russ Chem Bull 1997, 7, 1278.
[]: Comba P.; Hambley T.W. Molecular Modeling, VCH, Weinheim, New York, Basel, Cambridge, Tokyo, 1995.
[]: Kozelka J. in Metal Ions in Biological Systems, vol. 33, Sigel A.; Sigel H. Eds., Marcel Dekker Inc, New York, Basel, Hong Kong, 1996, p. 1.
[]: Tokmachev A.M.; Tchougréeff A.L.; Misurkin I.A. J Mol Struct (Theochem) 2000, 506, 17.
[]: Tokmachev A.M.; Tchougréeff A.L. Int J Quantum Chem 2001, Accepted.
[]: Surján P.R. in The Concept of the Chemical Bond, in Theoretical Models of Chemical Bonding, Part 2. Z.B. Maksi\'c Ed., Springer, Heidelberg, 1989, p. 205.
[]: Wu W.; McWeeny R. J Chem Phys 1994, 101, 4826.
[]: Tokmachev A.M.; Tchougréeff A.L. J Comp Chem 2001, 22, 752.
[]: Tokmachev A.M.; Tchougréeff A.L. Russ J Phys Chem 1999, 73, 259.
[]: Tokmachev A.M.; Tchougréeff A.L.; Misurkin I.A. Russ J Phys Chem 2000, 74, Suppl. 2, S205.
[]: McWeeny R. Methods of Molecular Quantum Mechanics, 2nd Edition, AP, London, 1992.
[]: Náray-Szabó G. Acta Chim Acad Sci Hung 1976, 40, 261.
[]: Smirnov V.I. Course of Higher Mathematics, vol. 3, Part. I, Moscow, 1958 [in Russian].
[]: Parks J.M.; Parr R.G. J Chem Phys 1958, 28, 335.

Table 0: Areas between energy profiles, kcal mol^-1 Å (or kcal mol^-1 deg.).

Variation range for a
geometry parameter	FAFO	TAFO	FATO	TATO	TOTA
(Å or deg.)
change of one C-H bond length in methane
0.994¸1.194	0.031	0.021	0.013	0.004	0.0006
0.794¸1.594	2.119	1.241	1.075	0.246	0.038
0.794¸3.094	73.973	17.618	56.527	2.121	0.207
change of one O-H bond length in water
0.880¸1.060	3.131	2.676	0.047	0.405	0.037
0.800¸1.200	7.904	6.732	0.229	1.112	0.113
change of two O-H bond lengths in water
0.880¸1.060	3.039	2.835	0.028	0.409	0.006
0.800¸1.200	9.325	7.978	0.273	1.230	0.014
change of angle H-O-H in water
101.0¸108.0	115.23	98.48	0.32	13.67	0.30
97.0¸112.0	246.64	211.00	0.81	29.29	0.66

Table 0: Heats of formation of ethane for different values of the C-C bond length calculated by the APSLG and TOTA methods, kcal mol^-1

r(C-C), Å	APSLG	TOTA	D
1.370	-8.456	-8.158	0.298
1.420	-15.585	-15.383	0.202
1.450	-17.902	-17.037	0.865
1.470	-18.739	-18.044	0.695
1.490	-19.065	-18.519	0.546
1.510	-18.923	-18.508	0.415
1.520	-18.689	-18.333	0.356
1.530	-18.354	-18.051	0.303
1.536	-18.106	-17.833	0.273
1.540	-17.922	-17.668	0.254
1.550	-17.397	-17.188	0.209
1.560	-16.785	-16.615	0.170
1.580	-15.312	-15.209	0.103
1.600	-13.538	-13.484	0.054
1.620	-11.493	-11.471	0.022
1.650	-7.977	-7.973	0.004
1.700	-1.138	-1.087	0.051

Table 0: Heats of formation of molecules calculated by the APSLG and TOTA methods with two sets of initial amplitudes for C-H bond, kcal mol^-1

Molecule	APSLG	TOTA(1)	D(1)	TOTA(2)	D(2)
CH₄	-8.414	-8.414	0.000	-8.397	0.017
H₂O	-55.648	-55.605	0.043	-55.605	0.043
C₂H₆	-18.106	-17.833	0.273	-17.844	0.262
C₃H₈	-23.438	-22.805	0.633	-22.869	0.569
C₄H₁₀	-26.951	-26.200	0.751	-26.295	0.656
C₅H₁₂	-32.093	-31.157	0.936	-31.283	0.855
cyclopropane	18.520	18.877	0.357	18.806	0.286
CH₃OH	-47.677	-47.014	0.653	-47.001	0.666
C₂H₅OH	-59.112	-57.905	1.207	-57.997	1.115
C₃H₇OH	-63.811	-62.386	1.425	-62.520	1.291

File translated from T_EX by T_TH, version 2.67.
On 19 Jul 2001, 17:56.

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

Generic Molecular Mechanics as Based on Local Quantum Description of Molecular Electronic Structure

Abstract

Introduction

Construction of generic MM schemes

Wave function and energy of underlying APSLG based QM method

Non-iterative tuning of ESPs

Bond (geminal) amplitudes

Hybridization matrices

Parametric forms for the energy

Results and Discussion

Conclusion

Acknowledgements

References

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