An Analytical Procedure to Evaluate
Electronic Integrals for
Molecular Quantum Mechanical Calculations
To be published on Physica-A
(2005)
2. Some Considerations about the Generalized Gaussian
Function
3. Technical Limitations of the usual Hartree-Fock Approach
5. Applications : and
BH Molecular Systems
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.
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.
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,
(1)
The q parameter comes from the Tsallis
entropy [10] for the micro canonical ensemble that is given by
(2)
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,
(3)
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,
(4)
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,
(5)
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.
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.
(6)
where, .
In the Hartree-Fock calculations it is important to
evaluate the Fock Matrix,
defined as :
(7)
here (mn|ls) is a two-electron
integral, the one-electron integral defined as,
(8)
In case of two centers, the
two-electron integral, as a function of the inter-atomic distance, is given by,
.
(9)
here,
.
(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.
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,
(12)
and more properly,
(13)
here amnls, bmnls , 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,
(14)
and
(15)
By this methodology, the total electronic energy is given as a function of the inter-atomic
distance through the two-electron integrals, as follows
(16)
more precisely,
or
(17)
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}.
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.,
(18)
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 .
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].
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.
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 .
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,
(20)
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.
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.
Acknowledgment
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.
References
[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. Ostlund, Modern
Quantum Chemistry,
[5] D.R. Hartree,
Proc.
[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,
[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