Return to Main-Page


An Analytical Procedure to Evaluate Electronic Integrals for

Molecular Quantum Mechanical Calculations


Kleber Carlos Mundim

Instituto de Química, Universidade de Brasília, 70919-970 Brasília, DF, Brazil.



To be published on Physica-A (2005) 

Abstract 2

1.  Introduction. 2

2. Some Considerations about the Generalized Gaussian Function. 3

3.  Technical Limitations of the usual Hartree-Fock Approach. 4

4.  The Alternative Approach. 6

5.   Applications :   and BH Molecular Systems. 7

6.   Conclusions and Remarks. 12



The main contribution of this paper is to provide an alternative strategy to reduce the number of bi-electronic integrals in ab-initio quantum mechanics calculations. An analytical procedure to evaluate the energy of a molecule as well as two-electron integrals is proposed. This approach is based on the generalized exponential function and is particularly advantageous because it reduces substantially the CPU time in quantum mechanical calculations. Some examples of the effectiveness of this methodology are presented. It is important to point out that the new methodology is applicable to any kind of molecular system including relatively large molecular systems in the context of the Hartree-Fock (HF) and Density Functional (DFT) theories.



1.  Introduction


Theoretical and computational molecular sciences are necessary tools in developing molecular-level descriptions of chemical, physical and biological processes in natural and contaminated systems. A primary difficulty is that as the size of the systems becomes larger and more complex, the empirical methods become less reliable. One promising approach to avoid empiricism is to employ first principles (ab-initio) methods of quantum mechanics (QM). Solving the Schrödinger equation in QM is hence a major challenge. Part of the computational problem is in the construction of the electronic Hamiltonian. The Hamiltonian contains two terms that act on one electron at a time, the kinetic energy and the electron-nucleus attraction, and a term that describes the pair-wise repulsion of electrons. The latter depends on the coordinates of two electrons at the same time and has turned out to be a practical computational bottleneck accessible only for very small molecular systems. Since the time and memory/storage requirements vary with problem size as O(n4), a significant bottleneck in such calculations and a good candidate for algorithmic optimization is the 2-electron integral transformation, where n é the number of contracted Gaussian functions. Although the QM equations governing chemistry are well known, their mathematical solutions are heavily computer time demanding and hence impractical for most biochemical problems. Therefore, any methodology reducing computer time is welcome. In past years, different quantum mechanical methodologies have been proposed [1,2 ] to reduce the computational time and have met with partial success.

Here an alternative methodology for the calculation of multi-centric electronic integrals is being proposed using an analytical and universal function based on the generalized Gaussian function named here as qGaussian (to be defined below). The integrals become analytical functions of the inter-atomic distances. Therefore, when estimating certain quantities such as molecular energy, qGaussians avoid new calculations of the integrals: they are simply another value of the corresponding function. In the usual approach CPU time increases with n4 whereas in the proposed method the CPU time scales linearly with n. In the usual procedure the catastrophic dependence of the rank the Hamiltonian or Fock matrix with n4 two-electron integrals is a severe bottleneck for petaFLOPS computing time. Thus, by reducing the number of two-electron integrals, QM turns out more feasible for the description of macromolecular systems.

In section 2 we present some considerations about the generalized Gaussian function. In section 3, some technical limitations of the usual HF procedure which are the motivation for introducing an alternative methodology for QM calculations are being presented, and in the last sections a few examples of the effectiveness of the present methodology are presented. It is important to emphasize that in spite of the fact that the chosen examples are small systems, the method is equally applicable to systems of any size, including bio-molecules, solid materials and solutions, within the CI [3], MCSCF [4], HF [5,6], post-HF [7]  and DFT [8] theories. A computational code for quantum mechanics calculation using the new HF-procedure is being developed and with this code studies of relatively larger molecular systems should be possible.




2. Some Considerations about the Generalized Gaussian Function


Before starting the application sections, let us give a brief introduction to the generalized Gaussian function. The generalization of the exponential function was proposed by Borges [9] based on the Tsallis statistic [10], in the following way,


            The q parameter comes from the Tsallis entropy [10] for the micro canonical ensemble that is given by


were k is a positive constant and W is the total number of microscopic possibilities of the system. In the limit , the q-entropy goes to Boltzmann entropy, as follows,



In this sense, the generalized q-exponential function (Eq. 1) is equivalent to the usual exponential function in the limit . The Tsallis Statistic has been successfully applied in the description of a variety of problems, particularly in global optimization processes [8-13]. Following this idea we introduce the generalized Gaussian function proper to the Hartree-Fock problem as,


where,  and  are adjustable parameters. In some case it is necessary to introduce the additional parameter  to better fit the one- and two-electron integrals. In this case the Eq. 4 will be,


It is important to point out that the success of the fitting is due, in part, to the qGaussian function is significantly more delocalized than the standard representation, as shown in the Fig. 1.

Fig. 1- The qGaussian function   is plotted for different values of q-parameter. In case of  the qGaussian is equivalent to the usual Gaussian function, that is represented by black points.




3.  Technical Limitations of the usual Hartree-Fock Approach


Quantum chemistry, particularly ab-initio molecular electronic structure theory has well-defined scaling parameters for the computation of total molecular energies, wave functions, and other derived properties. Such calculations correspond to approximate solutions of the quantum mechanical equations of motion (under the Schrödinger or Dirac form) for a system consisting of atoms, electrons and nuclei, the latter being fixed in space (Born-Oppenheimer approximation).

The most common approximation is the Hartree-Fock where molecular orbitals are written as linear combinations of atomic orbital (LCAO-MO), the expansion coefficients being optimized in a self-consistent way. The primitive basis functions for polyatomic molecules are most often chosen to be Cartesian Gaussian functions , where g is a constant and x,y,z, and r measure a nucleus-electron distance. Contracted Gaussian basis functions (CGTFs) [14-18], with fixed linear combinations of primitive functions that are based on atomic wave functions, are commonly used to reduce the effective number of coefficients to be optimized in HF and post-HF procedures. The one electron integrals giving the overlap between different states,  and , i.e.


where, .

In the Hartree-Fock calculations it is important to evaluate the Fock Matrix, defined as :


here (mn|ls) is a two-electron integral, the one-electron integral defined as,


In case of two centers, the two-electron integral, as a function of the inter-atomic distance, is given by,

.              (9)


.                                           (10)

On the right side of the Eq. 9, note that the number of the two-electron integrals (Áijkl) scales with n4. More precisely, the total number of one- and two-electron integrals scales to some power of the number of primitive Gaussian functions as,

           and                      .                  (11)

The “bottleneck” of this type of calculation is usually the amount of time consumed to evaluate the two-electron integrals.




4.  The Alternative Approach

         The approach proposed here is used to reduce the number of two-electron integrals (Eq. 9) using qGaussian functions defined as follows, in case of two-centers,


and more properly,


here amnlsbmnls , qmnls ,  and Rmnls are parameters adjusted to fit the sum of all two-electron integrals in the contracted Gaussian form. It is important to point out that the two-electron integrals depend on up to four centers. In case of three and four centers, the two-electrons integrals are represented by  hyper-surface that depends on the inter-atomic distances between the centers. In particular, for three atomic centers A, B and C the two-electronic integral will be . See one example for this case in the last section.

In the usual procedure it is necessary to evaluate n4 terms as in Eqs. 9 and 10, but in the present approach we have now only one term, i.e. an analytical function for the two-electron integral. As Eq. 13 expresses the value of the two-electron integral as a function of the inter-nuclear distance, the value of any integral Imnls  can be obtained, at any given distance r, just by inserting this value in Eq. 13.

         Similarly, it is possible to evaluate the  and  element of the matrix in terms of qGaussian functions respectively as,




         By this methodology, the total electronic energy is given as a function of the inter-atomic distance through the two-electron integrals, as follows


more precisely,



here Hmn is the one-electron part of the Hamiltonian.

According to the proposed methodology, the two-electron integrals are represented as simple analytical functions of the inter-atomic distances r and the process of optimizing geometries and/or of PES construction is therefore greatly simplified. Hence, the total energy of the molecule becomes a function of the set of coefficients {} and of the integrals, which in turn depend parametrically on the inter-nuclear distances {r}.  



5.   Applications :   and BH Molecular Systems

In the first application we evaluate the total energy for the HeH+ system using the usual Hartree-Fock and our procedure. The usual calculation is done expanding the atomic orbital (1s) in a linear combinations of three Gaussian function (STO-3G), i.e.,


where the  and  parameters are given in reference [4]. To evaluate the total energy it is also necessary to compute the following integrals: the overlap matrix , the one-electron part   and the two-electron integrals that are expressed by

 .                    (19)

Since, in this case, each index m,n,l and s can assume two values (1 and 2), there will be 16 Imnls  integrals. The right side of the above equation contains 34 terms. In this example, to construct the 16 integrals it is necessary to evaluate  Gaussian integrals like Áijkl, for each inter-atomic distance . Normally, if one is interested in optimizing the molecular geometry or in generating the potential energy surface (PES) for the molecule, the computational time increases drastically because it would be necessary to re-evaluate all of the terms  in Eq. 14 for each inter-atomic distance .

Next , , and have been computed using the proposed procedure to evaluate the total energy for the HeH+ system. To fit these integrals, only the parameters amnls, bmnls , qmnls , Rmnls and gmnls are needed. In Fig. 2(a-c) the value of some of the  integrals as well as the total energy for the HeH+ system obtained using the well known  STO-3G and our qGaussian functions has been compared, with a  chi-square error equal to .


Fig. 2-  Integrals STO-3G and fitted qGaussian for (a) overlap matrix element S12,

(b)  two-electron integral matrix element I2121, and

(c) total energy for HeH+ using STO-3G and fitted qGaussian integrals.


In Fig. 3(a-c) a similar analysis is done for the BH molecule using a STO-6G basis set. For this basis it is necessary to evaluate  non-zero two-electron integrals (Áijkl) to build up the full matrix Imnls, by the usual procedure. On the other hand, by using the proposed methodology, only 83 qIntegrals are needed to build up the non-zero values of Imnls, as defined in Eq. 8. An example of the quality of this fitting is shown in Fig. 3(a-c) and Fig. 4. All the calculations for STO-6G basis set were carried out using GAMESS ab-initio program [19].

Fig. 3- Integrals evaluated through STO-6G basis set and its qGaussian fitting for BH molecule: (a) Overlap matrix element S16. (b) One-electron integral, matrix element H16, and (c) Two-electron integral (element matrix I6221).

Using the qIntegrals, at the equilibrium geometry (=2.229 Bohr), the molecular energy for BH is -24.9893657480 Hartrees, which is exactly the same value obtained by the usual Hartree-Fock calculation. The total energy is plotted in the Fig. 4. 

Fig. 4- Total energy for BH molecule using a STO-6G base.


The new procedure to evaluate the two-electron integrals for large base set is also presented. Fig. 5 shows one of the bi-electronic integrals for BH molecule using a double-zeta base set, which was fitted through a qGaussian function with a chi-square error equal to .

Fig. 5- Two-electronic integral for BH molecule using a Double-zeta base set.


In case of three and four centers, the bi-electronic integrals are represented by hyper-surface that depends on the inter-atomic distances. For example, in case of the water molecule (H-O-H) all three-centers integrals are expressed by , that is,


 The Fig. 6 shows one three-centre bi-electronic integral for the atomic orbital , ,  and . The calculation was done using a STO-6G basis set. The surface in Fig. 6 is a projection of the  in the plane formed by the distances  and .

Although the proposed methodology was illustrated by means of very simple systems, it can clearly be extended to systems of any sizes. In fact, this methodology should be particularly useful for studying large systems and also for performing molecular quantum dynamics calculations at much reduced costs. A general quantum mechanical computational code, using the new procedure, is in development. With this computational code we will be able to use the alternative Hartree-Fock approach in the calculations of relative large molecular systems.

Fig. 6- Bi-electronic integral in case of three-centers was evaluated using a STO-6G basis set and fitted through a qIntegral. The inter-atomic distance RAB is given in a.u.


6.   Conclusions and Remarks


Based on the results of the present study one could draw the following general conclusions:

- That the new methodology is applicable to any kind of molecular system including relatively large molecular systems in the context of the Hartree-Fock (HF), a post Hartree-Fock, Multi-Configuration Self-Consistent Field (MCSCF), the configuration interaction (CI) as well as the Density Functional (DFT) theories.

- Comparing the usual HF calculations with those using the q Integrals, it is clear that the alternative procedure is computationally more advantageous for larger GTO basis set functions. Because the number of bi-electronic integrals scale with n4 in the usual case, whereas in the present methodology this number scale with n.

-  For the general case, by using the qGaussian procedure, the full Hartree-Fock problem, as well as all molecular properties like total energy and potential energy surface (PES), can be determined without calculating any one- or two-electron integrals. It is believed that this procedure will be suitable and useful for molecular dynamics studies in relatively large molecular systems.



The author acknowledges fruitful discussions with L.A.C. Malbouisson, M.A. Chaer Nascimento, R. Garg. Some financial support from CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico- Brazil) and CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior-Brazil) is thankfully acknowledge.



[1] H.J. Werner, F.R. Manby, P.J. Knowles, J. Chem. Phys. 118 (18) (2003) 8149.

[2] N.A. Anikin, V.M. Anisimov, V.L. Bugaenko, V.V. Bobrikov, A.M. Andreyev,  J. Chem. Phys. 121 (2004) 1266-1270.

[3]   J. Paldus, J. Chem. Phys. 61 (1974) 5321.

[4]  A. Szabo,  N.S. OstlundModern Quantum Chemistry, Dover Publications, NY, (1996).

[5]  D.R. Hartree, Proc. Cambridge Phil. Soc. 24 (1928) 111  and  25 (1929) 225.

[6]  V. Fock, Z. Physik, 61 (1930) 126.

[7]  J.A. Pople, M. Head-Gordon, J. Chem. Phys. 87 (1987) 5968.

[8]  L.J. Sham, W. Kohn,  Phys. Rev. 145 (1966) 561.

[9]  E.P. Borges,  J. Phys. A: Math. Gen. 31  (1998) 5281-5288.

[10]  C. Tsallis,  J. Stat. Phys.  52 (1988) 479-487.

[11]  K.C. Mundim, C. Tsallis,  Int. J. Quantum Chem. 58: (4) (1996) 373-381.

[12] M.A. Moret, P.G. Pascutti, K.C. Mundim, P.M. Bisch , E. Nogueira Jr., Phys. Rev. E  63  (2001) 020901 (R).

[13] L.E. Espinola,  J.J. Soares Neto, K.C. Mundim, M.S.P. Mundim, R. Gargano, Chem. Phys. Letters 361 (2002)  271-276.

[14] S.F. Boys,  Proc. R. Soc.  A 200 (1950) 542.

[15]  W.J. Hehre, R.F. Stewart, J.A. Pople, J. Chem. Phys. 51 (1969) 2657.

[16]  M.S. Gordon, J.S. Binkley, J.A. Pople, W.J. Pietro, W.J. Hehre, J. Am. Chem. Soc. 104 (1982) 2797.

[17]  T.H. Dunning,  J. Chem. Phys. 90 (1989) 1007

[18]  J.S. Binkley, J.A. Pople, W.J. Hehre,  J. Am. Chem. Soc. 102 (1980) 939.

[19]  GAMESS-2000 Iowa State University, M.W.Schmidt, et al., J.Comp. Chem. 14 (1993) 1347-1363.