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

Fast NDDO Method for Molecular Structure Calculations Based on Strictly Localized Geminals

Dedicated to Prof. R. Hoffmann on the occasion of his 65th birthday

Abstract

Introduction

Almost fourty years ago the Series of seminal papers by R. Hoffmann Hoff1,Hoff2,Hoff3 has launched the enterprise of all valence calculations of organic molecules with use of the single determinant approximation for the many-electron trial wave function. The progress of this enterprise was enormous and success in calculation of small molecules achieved by modern ab initio methods is impressive. At the same time a vast number of real molecular systems of chemists' interest remain practically unaccessible by quantum chemistry due to enormous computational costs. The reason of this is an unpleasant growth of computational resources required by ab initio techniques when the system size grows (N⁵¸N⁷, where N is the dimension of the one-electron states basis involved in the calculation). In the case of semiempirical self consistent field (SCF) methods the computational resources also grow as N³ due to diagonalization involvement in the procedure. Therefore, even the application of semiempirical methods to construction of potential energy surfaces (PESs) for large systems (especially, of biological significance) may as well become problematic.

Two principal types of solution of the above problems are proposed in the literature. The first one is to construct the methods with a weaker dependence of the required computational resources on the size of the system. It was shown that for large molecules calculation times should increase as slowly as N^1.3 []. There is a number of attempts to construct the methods with optimal N-scalability properties [,,,,]. Such schemes were proposed for different types of quantum chemical approaches - DFT, ab initio, and semiempirical ones. The most popular approach is based on ignoring some matrix elements of the Fock matrix between the basis functions centered on distant atoms [,,]. Another way to reduce computational costs can be based on the direct determination of localized Hartree-Fock orbitals []. Making use of local one-electron states allows to construct effective schemes for taking into account the electronic correlation. Acceleration of computation can be also achieved by pseudodiagonalization [] or by special tricks with wave function [,]. The methods with approximately linear N-scalability were applied to very large molecular systems [].

The second type of solution is based on construction of the so called hybrid QM/MM schemes where different parts of the system are treated with different level of quality and, therefore, with different computational costs. The reason of these methods is that the chemical transformation usually occurs only in a small region (reactive center) while the environment only slightly modifies the PES. After the pioneering work by Warshell and Levitt Warshel hybrid techniques become very popular Bakowies,Antes,Assfeld,Froese,Gao. At the same time the important question about junction between different subsystems is solved in these methods ad hoc and not by the successive separation of variables.

Recently we proposed a special procedure of deriving the junction between subsystems described by quantum mechanical and molecular mechanical methods []. It is based on the trial wave function having the form of the antisymmetrized product of strictly localized geminals (APSLG). This form is taken as an underlying one for the molecular mechanical part of the molecule. At the same time the APSLG wave function can be used to construct the O(N)-scalable methods since the strictly local structure of the trial wave function allows to eliminate the diagonalization of the matrices of the size N×N. In this paper we exploit this possibility.

The use of local orbitals is a direct way to good scalability properties []. These orbitals can be obtained, for example, by orthogonal transformations of canonical MOs [,]. At the same time these orbitals do not centered on the pair of atoms and have ''tails'' on other atoms. The localized orbitals without ''tails'' are called strictly localized molecular orbitals and can be obtained in special variational SCF determination of electronic structure []. The well known PCILO scheme [] also uses local orbitals but treats them by a sort of perturbation technique. In the framework of the present work we use special type of local orbitals - hybrid orbitals (HOs) obtained by transformations of minimal basis sets for each atom [].

The APSLG approximation is similar to other pair theories like the generalized method of valence bonds (GVB) []. The difference is in the way the one-electron states to be used for the bond function are chosen, and, therefore, in the degree of wave function localization. The ab initio version of the APSLG approach [] uses non-variational Pauling's HOs [] for constructing the geminals. This approach was applied only to a small number of very simple molecules. The results do not allow to make a conclusion about general applicability of the scheme to large molecules, since, even in the case of CH₃F molecule the APSLG electronic energy is significantly higher than that of the SCF approach Poirier. The calculations of the C-H bond dissociation show that the APSLG energy is lower than the GVB one for equilibrium bond distance but it is higher in the totally dissociated state [].

The semiempirical implementation of the APSLG approach [] is performed with simple molecular Hamiltonian of the MINDO/3 type MINDOAppr,MINDO. The drawbacks of this approach are inherent from the MINDO/3 approximation: heats of formation for unsaturated organic compounds are too negative and for branched molecules are very positive; chemical bonds between atoms with lone electron pairs are too short, bond angles are not well reproduced. There are different possible ways to cure these drawbacks: to further adjust parameterization; to take into account perturbation corrections to the wave function/energy, or to use more elaborated Hamiltonians having new interactions. We choose the last possibility since in the case of the SCF approach it recommended itself as quite successful. Thus, we use in this paper the Hamiltonians of NDDO family taking in details the two-center Coulomb interactions as implemented in the well known MNDO [,], AM1 [], and PM3 [] schemes. The principal difference between the MINDO and NDDO schemes and theoretical justification of the NDDO Hamiltonian are discussed in details in Ref. []. To summarize, we try to obtain a quantum chemical method with weak dependence of computational costs on the size of the system by replacing the SCF wave function by the APSLG one. We also try to reach reliable (not worse than in the SCF method) description of molecular properties like heats of formation and molecular geometries. On the way the change of wave function should allow to cure such an unpleasant property of the SCF wave function as the incorrect asymptotic behaviour under cleavage of chemical bonds. The paper is organized as follows: in the next Section we consider general theoretical principles underlying the APSLG approach and its NDDO implementation; then we describe parameterization procedure and calculations of molecular properties using the APSLG-NDDO schemes; these results are discussed and, finally, the conclusions about general applicability and advantages of the procedure proposed are given.

Theory

The electronic wave function in the APSLG approximation has the form:

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

(0)

where the mth geminal is presented by a linear combination of singlet two-electron configurations given by products of two operators creating electrons on HOs corresponding to ''right'' (r) and ''left'' (l) atoms of chemical bond with the spin projections s(=a,b):

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

(0)

These geminals are mutually orthogonal and satisfy a normalization condition:

á 0| gmgm+| 0 ñ = um2+vm2+2wm2 = 1.

(0)

The amplitudes u_m, v_m and w_m correspond to two ionic configurations and one covalent (of Heitler-London type). This form of wave function was originally proposed by Weinbaum []. In the case of electron lone pair only one configuration survives (for the sake of definiteness we assume it to be the right-end ionic contribution) and the geminal has the form:

gm+ = rma+rmb+

(0)

with normalization condition automatically fulfilled.

The important question is about particular construction of HOs r_m and l_m . These one-electron functions form the carrier space for geminals (the so-called Arai subspaces []). Strictly local character of geminals (and the wave function itself) assumes that the HOs have not ''tails'', i.e. they are expressed through the basis functions centered on one atom only. In the case of minimal basis set used the orthogonality of geminals immediately leads to the mathematical structure of HOs as produced by orthogonal (SO(4)) transformations of the initial set of atomic orbitals (AOs) for each ''heavy'' (non-hydrogen) atom. These transformations h^A act in the four-dimensional spaces spanned by one s- and three p-AOs:

tms+ = å i Î A  hmiAais+,

(0)

where t denotes a HO (right r or left l) located on the atom A.

The transformation of the basis set produces transformation of molecular integrals entering the NDDO Hamiltonian. These integrals in the HO basis are linear combinations of the same type integrals in the AO basis with coefficients taken as products of the elements of transformation matrices. Here we produce only integrals which are actually necessary for estimation of the electronic energy. The attraction of an electron on the HO t_m to its own core is:

UtmtmA = å i Î A  (hmiA)2UiiA

(0)

or, using properties of SO(4) matrix h^A, it can be expressed as function of the weight of the s-AO only:

UtmtmA = UppA+(UssA-UppA)(hmsA)2.

(0)

Two-electron one-center molecular integrals for the sp-shell can be expressed through five Slater-Condon parameters []. The integrals actually entering the expressions for the energy (see below) can be expressed through the expansion coefficients for the s-function only:

(tmtm | tmtm)A = C1A+C2A(hmsA)2+C3A(hmsA)4, gtktm¢A = 2(tktk | tm¢tm¢)A-(tktm¢ | tm¢tk)A = = C4A+C5A[(hmsA)2+(hksA)2]+C3A(hmsAhksA)2,

(0)

where the combinations of the Slater-Condon parameters are introduced:

C1A = F0A(pp)+4F2A(pp), C2A = 2F0A(sp)+4G1A(sp)-2F0A(pp)-8F2A(pp), C3A = F0A(ss)-2F0A(sp)-4G1A(sp)+F0A(pp)+4F2A(pp), C4A = 2F0A(pp)-7F2A(pp), C5A = 2F0A(sp)-G1A(sp)-2F0A(pp)+7F2A(pp).

(0)

The formulae Eqs. (7), (8) show that one-center molecular integrals (and therefore one-center energy) are independent of the directions of HOs. The dependence on the whole structure of HOs is given by two-center molecular integrals. Diagonal element of attraction of an electron on the HO t_m to other cores is:

Vtmtm,BA = å i,j Î A  Vij,BAhmiAhmjA.

(0)

Other matrix elements depend on the form of HOs for pairs of atoms. The resonance (electron transfer) matrix elements between the ''right'' and ''left'' HOs of the mth bond have the form:

brmlmAB = å i Î A  å j Î B  hmiAhmjBbijAB.

(0)

The matrix elements of the Coulomb repulsion of electrons located on different atoms A and B are:

(tm1tm1 | tm2¢tm2¢)AB = å \Sb i,j Î A k,l Î B\endSb (ij | kl)ABhm1iAhm1jAhm2kBhm2lB.

(0)

In the case of multiple bonds additional two-center matrix elements become necessary:

(tm1tm2¢ | ~ t   ¢ m2  ~ t   m1  )AB = å \Sb i,j Î A k,l Î B\endSb (ij | kl)ABhm1iAhm2jAhm2kBhm1lB,

(0)

where [(t)\tilde] = l if t = r and [(t)\tilde] = r if t = l. These matrix elements correspond to electron transfers within a pair of single bonds between the same pair of atoms. It should be noted that the molecular integrals in the AO basis entering the above expressions are taken the same as they are in the corresponding SCF predecessor procedures MNDOAppr,MNDO,AM1,PM3.

The electronic Hamiltonian for a molecular system can be written in the transformed basis (basis of HOs). It is a sum of one-center and two-center contributions:

H = å A  HA+ 1 2 å A ¹ B  HAB.

(0)

In the second quantization notation they are:

HA = å tm Î A  (UtmtmA+ å B ¹ A  Vtmtm,BA) å s  tms+tms+ + å \Sb tk,tm¢ Î A k < m\endSb å B ¹ A  Vtktm¢,BA å s  (tks+tms¢+h.c.)+ + 1 2 å \Sb tm1,tm2¢ Î A tm3¢¢,tm4¢¢¢ Î A\endSb(tm1tm2¢ | tm3¢¢tm4¢¢¢)A å st  tm1s+tm3t¢¢+tm4t¢¢¢tm2s¢

(0)

and

HAB = - å \Sb tm1 Î A tm2¢ Î B\endSb btm1tm2¢AB å s  (tm1s+tm2s+h.c.)+ + å \Sb tm1,tm2¢ Î A tm3¢¢,tm4¢¢¢ Î B\endSb(tm1tm2¢ | tm3¢¢tm4¢¢¢)AB å st  tm1s+tm3t¢¢+tm4t¢¢¢tm2s¢,

(0)

where h.c. stands for hermitean conjugation.

The total energy of a molecule is a sum of electronic energy and that of core-core interaction. The specific forms of the last term are respectively taken without changes from Refs. [,,]. The electronic energy is obtained by averaging the Hamiltonian Eq. (14) over the APSLG trial wave function Eq. (1). To do so we evaluate first the intrabond elements of one- and two-electron density matrices:

Pmtt¢ = á 0| gmtms+tms¢gm+| 0 ñ , Gmtt¢ = á 0| gmtmb+tma¢+tma¢tmbgm+| 0 ñ , Pmrr = um2+wm2, Pmll = vm2+wm2, Pmrl = Pmlr = (um+vm)wm, Gmrr = um2, Gmll = vm2, Gmrl = Gmlr = wm2.

(0)

These matrix elements are spin-independent. Taking into consideration different contributions to the Hamiltonian Eqs. (15) and ( 16) one can represent an electronic energy as a sum of five different terms:

E = Ec-attr+Eoc-rep+Eres+Etc-rep+Em-b.

(0)

The first contribution is attributed to electron repulsion to cores:

Ec-attr = 2 å A  å tm Î A  æ è UtmtmA+ å B ¹ A  Vtmtm,BA ö ø Pmtt.

(0)

The one-center electron-electron repulsion is a sum of contributions from repulsion of electrons on one or two different HOs:

Eoc-rep = å A  å tm Î A  (tmtm | tmtm)AGmtt+2 å \Sb tm1,tm2¢ Î A m1 < m2\endSb gtm1tm2¢APm1ttPm2t¢t¢.

(0)

The overall contribution to the energy from the resonance interaction is:

Eres = -4 å A < B  å m Î AB  brmlmABPmrl,

(0)

where notation m Î AB means that the mth bond is one between atoms A and B. The interatomic contribution from repulsion of electrons also depends on the type of interaction - between different bonds (or lone pairs) or inside one chemical bond:

Etc-rep = 2 å A < B  å tm1 Î A  å tm2¢ Î B  (tm1tm1 | tm2¢tm2¢)AB× ×[2(1-dm1m2)Pm1ttPm2t¢t¢+dm1m2Gm1rl].

(0)

All All the above contributions have their analogues in the APSLG construction [] based on the MINDO [,] implementation. The NDDO Hamiltonian is superior to MINDO one. It leads to arising a specific contribution corresponding to interaction of single bonds constituting one multiple bond:

Em-b = -4 å \Sb tm1tm2¢ Î A m1 < m2\endSb (tm1tm2¢ | ~ t   ¢ m2  ~ t   m1  )ABPm1rlPm2rl.

(0)

The electronic (and total) energy thus depends on the electronic structure parameters (ESPs) - amplitudes u_m, v_m and w_m of Eq. (2) through elements of density matrices Eq. (17) and elements of SO(4) matrices h^A defining hybridized one-electron states through molecular integrals. The total number of independent amplitudes is 2M (M is a number of chemical bonds) due to normalization condition Eq. ( 3) imposed on the geminal amplitudes. The total number of hybridization defining parameters is 6L (L is a number of heavy atoms) since the SO(4) group is a six-parametric one. We use parametric representation of the SO(4) group based on six subsequent Jacobi rotations in two-dimensional subspaces of a four-dimensional space of AOs. Therefore, six parameters are the corresponding angles of Jacobi rotations. The determination of the ESPs is performed by using a variational principle by a series of iterations. The first step is calculation of geminal amplitudes by diagonalizing of 3×3 effective bond Hamiltonians for each geminal representing a chemical bond. The next step is a series of energy minimizations with respect to sextuples of parameters defining SO(4) transformations for each heavy atom. These minimizations are performed with use of analytical gradients of the energy with respect to the Jacobi angles. The alternating diagonalizations and minimizations are performed until convergence. The number of iterations remains approximately constant with increasing of the size of the molecular system. The procedure of equilibrium geometry determination based on the analytical gradients of the total energy with respect to geometric parameters is also implemented.

Results and Discussion

In the previous Section we have constructed a scheme of determination of the ESPs for APSLG-NDDO method. Clearly, this scheme avoids diagonalization of N×N matrices. Number of elementary steps (construction and diagonalization of effective bond Hamiltonians and minimizations of energy with respect to sextuples of hybridization angles) is proportional to the size of the system. Each step however has a contribution requiring computational resources proportional to the size of the system with small coefficient. Therefore, the scaling of computational costs is almost linear. If the molecular integrals between the basis functions centered on distant atoms are ignored like in Ref. [] the scaling becomes precisely linear. We do not give here any benchmark calculations since they are platform dependent. At the same time we note that the comparisons of SCF and APSLG calculation times are given for MINDO/3 implementation in Ref. APSLGJCC. It was shown that for system with 122 basis functions the APSLG procedure is 30 times faster than the SCF one.

The change of the trial wave function leads to changes in the calculated properties. It should be noted that the difference in the total energy for the SCF and APSLG wave functions (with same form of the Hamiltonian and parameterization) can be understood as a sum of two effects: better account of static intrabond electron correlation and neglect of interbond delocalization (electron transfer) in the APSLG scheme. The former lowers the energy while the latter increases it. In the case of H₂ molecule the interbond delocalization is absent and we always obtain lowering of the energy. We should state that the method after change of trial wave function without change of parameterization remains to give sane results. It allowed not to perform the total re-parameterization but to restrict ourselves by only slight tuning of parameters. We have studied three mostly known NDDO schemes: MNDO, AM1, and PM3 and restricted ourselves to molecules containing H, C, N, O, and F atoms. The SCF results on these compounds are given in Refs. [,,,,,]. We assume that the change of the wave function mostly affects the two-center contributions to the energy. So, we attempted to reach the reliable results by tuning only very small subset of parameters related to the resonance. The new resonance parameters are given in Table 1 and are compared with the SCF values. For all atomic orbitals the corrections are small. It should be noted that using the MINDO/3 type scheme for the resonance integrals which is more simple for parameterization is not possible here since the relation between different types of resonance integrals (b_ss and b_sp) is fixed by the AO ionization potentials which led to absurd molecular geometries - plane ammonia, linear water molecule etc.

The calculated and experimentally observed heats of formation are given in Table 2 for all three parameterizations and both (SCF and APSLG) trial wave functions for a set of typical molecules taken from Refs. [,]. The numerical results show that the APSLG parameterization is internally consistent. The estimates of the heats of formation obtained in the framework of the SCF and APSLG methods are of the same quality. It is seen that the use of the NDDO Hamiltonian in combination with the APSLG trial wave function cures the main problems of the APSLG-MINDO/3 approach - bad description of branched and unsaturated molecules. At the same time we note that problems of SCF-NDDO schemes are mirrored in the APSLG implementation. It is well seen on the example of cubane molecule where both the SCF-MNDO and APSLG-MNDO methods give very large error, while both the SCF-AM1 and APSLG-AM1 methods predict experimental value of the heat of formation relatively well. It should be noted that in the case of molecules with triple bonds (especially, N₂) the APSLG method leads to heats of formation significantly smaller than the SCF method. It is due to small importance of interbond electron transfer as compared to the intrabond correlation due to symmetry of the problem.

Another important characteristic of the quality of the quantum chemical method is its ability to reproduce correctly the parameters of molecular structure. The calculated and experimental geometric parameters are given in Table 3 for typical, most characteristic, and most difficult cases. The numerical results show that the APSLG-based method suits better to reproduce the molecular geometries than the SCF one. Using the APSLG scheme allows to cure significant problems of the SCF approach - incorrect description of torsion angles in cyclobutane and hydrogen peroxide molecules. In the case of ab initio SCF approach the acceptable result for cyclobutane can be achieved only by using large basis sets with polarization functions Cremer. Moreover, the description of bond lengths in many cases is significantly improved by taking into account the intrabond correlation (for example, the N-N bond in hydrazine, the F-O bond in F₂O). At the same time we find that using of the correlated APSLG wave function typically leads to increasing of bond length for atoms with electron lone pairs and to diminishing of valence angles (see NH₃ and H₂O molecules) in comparison with the SCF scheme. In some cases such one-directional change of geometry parameters can lead to a worse agreement with experiment. These changes should be considered as characteristic effect of the trial wave function not depending on the parameterization.

The APSLG trial wave functional operates with quantities widely used in chemistry but not in the quantum chemistry - bond characteristics and the hybridization parameters. Here we demonstrate how these characteristics determined on the ground of variational principle correspond to chemical intuition. We can rewrite the geminal expression Eq. (2) in the form []:

gm+ = Im Ö2 é ë Ö   1+lm   rma+rmb++ Ö   1-lm   lma+lmb+ ù û + + Cm Ö2 [ rma+lmb++lma+rmb+] ,

(0)

where I_m² can be considered as bond ionicity, C_m²( = 1-I_m²) as bond covalency, and l_m as bond polarity. In Table 4 we show the bond order (2P_m^rl), bond covalency and bond polarity for some typical chemical bonds. It can be seen that the bond order is very close to unity for very large number of s-bonds. The results show that different parameterizations lead to quite different description of chemical bond structure. For example, the AM1 scheme is prone to polarize bonds significantly; the PM3 scheme leads to electronegativity of the carbon atom exceeding that of the nitrogen atom and also predicts the fluorine and oxygen atoms to have very close electronegativities. It is seen that p-bonds are significantly more covalent and more polarizable than s-bonds in accordance with usual chemical intuition.

Another important parameter of electronic structure of the APSLG method is hybridization matrices. Here they are determined variationally. This approach allows to consider a question about the structure of multiple chemical bonds. Our calculations show that the s/p separated chemical bonds are more preferable than those of the bent-type (''banana'') in accordance with predictions of Ref. [] based on the full GVB consideration. Analogously we can confirm an old observation CoulsonBook that in the case of cyclopropane molecule the HOs are not directed along the C-C bonds.

It is typical to represent the ratio of s- and p-AOs in the HO as sp^x. We use this form of representation in Table 5 for some typical HOs. It can be seen that this type of representation is not very suitable for large and small weights of s-function due to its significantly non-linear character. The data of Table 5 demonstrate that the most characteristic hybridization patterns (sp³, sp², and sp) are reproduced in the calculation with small deviations due to non-equivalence of the bonds involved.

The form of HOs can be used for semiquantitative prediction of different properties. As it was shown in Ref. [] the value of pK for hydrocarbon dissociation is approximately a linear function of the weight of s-function in the HO representing the C-H bond. Here we construct these linear relations:

MNDO: pK = 74.13-1.01hms2; AM1: pK = 73.91-0.98hms2; PM3: pK = 73.64-0.94hms2.

(0)

The calculated by these relations data given in Table 6 are in good correspondence with the experimental ones.

It should be noted that the proposed APSLG approach has significant limitations. It applies in its current form only to systems with well defined chemical bonds without large delocalization. At the same time generalization to the molecules with groups with high level of electron delocalization is straightforward. The work in this direction is in progress now. Another problem is calculation of properties corresponding to significantly delocalized states like ionization potentials. At the same time in this case some types of configuration interaction procedures can be applied as it was demonstrated in Ref. [].

Conclusions

A semiempirical method for molecular electronic structure calculations is developed in this work. It is based on the trial APSLG wave function with three NDDO-type Hamiltonians - MNDO, AM1 and PM3. It can be considered as different from the SCF point on the plane with axes Hamiltonian - wave function. Using local one-electron states allows to achieve approximately linear dependence of computational costs on the size of the molecule. The APSLG method by construction has correct asymptotic behaviour under cleavage of chemical bonds. It is shown that the quality of results on the heats of formation and molecular geometries obtained by the SCF and APSLG approaches is comparable. In the framework of the approach proposed many chemically sensible concepts like hybridization and bond characteristics found their theoretical substantiation by determination on the ground of variational principle.

Acknowledgement.

This work has been performed with partial financial support of RAS through the grant # 6-120 dispatched by its Young Researchers' Commission. Financial support for AMT on the part of the Haldor Tops oe A/S is acknowledged as well.

References

[]

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

[]

Hoffmann, R. J. Chem. Phys. 1964, 40, 2474-2489.

[]

Hoffmann, R. J. Chem. Phys. 1964, 40, 2745-2745.

[]

Steigerwald, M. L.; Goddard III, W. A. J. Am. Chem. Soc. 1984, 106, 308-311.

[]

Kohn, W. Phys. Rev. Lett. 1996, 76, 3168-3171.

[]

Ordejón, P.; Artacho, E.; Soler, J. M. Phys. Rev. B 1996, 53, R10441-R10444.

[]

Hierse, W.; Stechel, E. B. Phys. Rev. B 1994, 50, 17811-17819.

[]

Millam, J. M.; Scuseria, G. E. J. Chem. Phys. 1997, 107, 5569-5577.

[]

Lee, T. S.; York, D.; Yang, W. J. Chem. Phys. 1996, 105, 2744-2750.

[]

Dixon, S. L.; Merz, K. M., Jr. J. Chem. Phys. 1997, 107, 879-893.

[]

Stewart, J. J. P. Int. J. Quantum Chem. 1996, 58, 133-146.

[]

Yang, W.; Lee, T. S. J. Chem. Phys. 1995, 103, 5674-5678.

[]

Rubio, J.; Povill, A.; Malrieu, J. P.; Reinhardt, P. J. Chem. Phys. 1997, 107, 10044-10050.

[]

Strain, M. C.; Scuseria, G. E.; Frisch, M. J. Science 1996, 271, 51-53.

[]

Cullen, J. M. Int. J. Quantum Chem. 1995, 56, 97-113.

[]

Derecskei-Kovacs, A.; Woon, D. E.; Marynick, D. S. Int. J. Quantum Chem. 1997, 61, 67-76.

[]

York, D., Lee, T. S.; Yang, W. J. Am. Chem. Soc. 1996, 118, 10940-10941.

[]

Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227-249.

[]

Bakowies, D.; Thiel, W. J. Phys. Chem. 1996, 100, 10580-10594.

[]

Antes, I.; Thiel, W. J. Phys. Chem. A 1999, 103, 9290-9295.

[]

Assfeld, X.; Rivail, J. L. Chem. Phys. Lett. 1996, 263, 100-106.

[]

Froese, R. D. J.; Humbel, S.; Svensson, M.; Morokuma, K. J. Phys. Chem. A 1997, 101, 227-233.

[]

Gao, J.; Amara, P.; Alhambra, C.; Field, M. J. J. Phys. Chem. A 1998, 102, 4714-4721.

[]

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

[]

Boys, S. F. Rev. Mod. Phys. 1960, 32, 296-299.

[]

Edmiston, C.; Ruedenberg, K. Rev. Mod. Phys. 1963, 35, 457-465.

[]

Smits, G. F.; Altona, C. Theoret. Chim. Acta 1985, 67, 461-475.

[]

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

[]

Randi\'c, M.; Maksi\'c, Z. B. Chem. Rev. 1972, 72, 43-53.

[]

Bobrowicz, F. W.; Goddard, W. A., III In Methods of Electronic Structure Theory; Schaefer, H. F., III, Ed.; Plenum Press: New York, 1977.

[]

Surján, P. R. In Theoretical Models of Chemical Bonding; Maksi\'c, Z. B., Ed.; Springer: Heidelberg, 1989; part 2.

[]

Poirier, R. A.; Surján, P. R. J. Comput. Chem. 1987, 8, 436-441.

[]

Tokmachev, A. M.; Tchougréeff, A. L. J. Comput. Chem. 2001, 22, 752-764.

[]

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

[]

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

[]

Dewar, M. J. S.; Thiel, W. J. Am. Chem. Soc. 1977, 99, 4899-4907.

[]

Dewar, M. J. S.; Thiel, W. J. Am. Chem. Soc. 1977, 99, 4907-4917.

[]

Dewar, M. J. S.; Zoebisch, E. G.; Healy, E.; Stewart, J. J. P. J. Am. Chem. Soc. 1985, 107, 3902-3909.

[]

Stewart, J. J. P. J. Comput. Chem. 1989, 10, 209-220.

[]

Weinbaum, S. J. Chem. Phys. 1933, 1, 593-596.

[]

Arai, T. J. Chem. Phys. 1960, 33, 95-98.

[]

Condon, E. U.; Shortley, G. H. Theory of Atomic Spectra; Cambridge, 1951.

[]

Dewar, M. J. S.; Rzepa, H. S. J. Am. Chem. Soc. 1978, 100, 58-67.

[]

Dewar, M. J. S.; Zoebisch, E. G. Theochem 1988, 49, 1-21.

[]

Cremer, D. J. Am. Chem. Soc. 1977, 99, 1307-1309.

[]

Ogliaro, F.; Cooper, D. L.; Karadakov, P. B. Int. J. Quantum Chem. 1999, 74, 223.

[]

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

[]

Maksi\'c, Z. B.; Eckert-Maksi\'c, M. Tetrahedron 1969, 25, 5113-5114.

[]

Tokmachev, A. M.; Tchougréeff, A. L.; Misurkin, I. A. Int. J. Quantum Chem. 2001, 85, 109-117.