Minhhuy Hô, Julio Manuel Hernández-Pérez

III. Nuclear-Electron Attraction Integrals

This article carries out the evaluation of nuclear-electron attraction energy integrals using Gaussian-type functions with arbitrary Cartesian angular values. As an example, we calculate the corresponding matrix for the water molecule in the STO-3G basis set.


Evaluating molecular integrals has been an active field since the middle of the last decade. Efficient algorithms have been developed and implemented in various programs. Detailed accounts of molecular integrals can be found in the references of [1]. In this article, the third in a series describing algorithms for evaluating molecular integrals, we detail the evaluation of the nuclear-electron attraction energy integrals from a more didactic point of view, following the approach of Rys, Dupuis, and King [2] as implemented in the OpenMol program [3].

The energy caused by the attraction between an electron in the region described by the overlap of the orbitals , and a nuclear of charge located at is expressed by the nuclear-electron attraction integral


in which is an unnormalized Cartesian Gaussian primitive.

Using the Gaussian product (see, for example [1]) and defining the angular part as :


The pole problem can be solved by the Laplace transform


which turns the integral into


where for now, we have ignored the factor . In the following steps, we will make certain modifications, knowing in advance that they will help simplify the expressions later on. We first reduce the upper limit of to unity by making the changes of variable (recall from the Gaussian product that ):




Replace , , and in , to get


We now multiply by the factor :


Again, by inserting , we get


Having arrived at the desired form, we reinsert the value of the angular part into the expression and separate the term enclosed by the curly brackets into three components , , and :


Defining as the function of the component inside the bracket,


and similarly for and , we rewrite the integral as


We will show that the integrand in the expression for is in fact an overlap between two one-dimensional Gaussians, and we may use the results that have been developed in [1]. First, we expand the exponential parts of the integrand


regrouping in terms of and , we have


which becomes


where and . These definitions let us compare this equation with the result of the Appendix, in which we see that equation (15) is simply




Substituting and ,


into equation (13), we have


Substitute this result into the definition of to get


The integral has the same form as a one-dimensional overlap integral where the integrand is a Gaussian function centered at with an exponential coefficient .

Recurrence Relations

From the observation above, we make use of the results developed for overlap integrals in [1]. For example, for ,


In particular, we have the transfer equations


The and functions take similar forms. The product is a polynomial in , and if we replace , then the integral in equation (12) is


where is the said polynomial. The integral is a combination of the Boys function (see, for example, Reference 4 of [1])


a strictly positive and decreasing function.


Aside from the obvious choice of using Mathematica to evaluate the Boys function, there are several ways of evaluating the integral. In practice, most programs store pretabulated values of the function at different intervals and interpolation is done as needed (e.g. by Chebyshev polynomials). Here we use the Gauss-Chebyshev quadrature numerical integration [4]. For simplicity, we have adopted almost verbatim the F77 code in [4, p. 46].

The function Nea evaluates the nuclear-electron attraction integral of two Gaussian primitives; here alpha, beta, RA, RB, LA, and LB are , , , , , and as defined earlier; RR is the nuclear position.

As in our two earlier articles [1, 5], we use the same data for the water molecule (, , the geometry optimized at the HF/STO-3G level). The molecule lies in the plane with Cartesian coordinates in atomic units.

In the STO-3G basis set, each atomic orbital is approximated by a sum of three Gaussians; here are their unnormalized primitive contraction coefficients and orbital exponents.

Here are the basis function origins and Cartesian angular values of the orbitals, listed in the order , , , , , , and .

Specifically, for the nuclear-electron attraction energy integral between the first primitive of the orbital of hydrogen atom 1, , the first primitive of the orbital of the oxygen atom, , and atom 1 () is


We have

From the Gauss-Chebyshev quadrature, the integral in equation (23) yields . The nuclear-electron integral (25) is . This is calculated as follows.

We would first need the normalization factor before evaluating the nuclear-electron energy matrix.


We have provided a didactic introduction to the evaluation of nuclear-electron attraction-energy integrals involving Gaussian-type basis functions by use of recurrence relations and a numerical quadrature scheme. The results are sufficiently general so that no modification of the algorithm is needed when larger basis sets with more Gaussian primitives or primitives with larger angular momenta are employed.


Consider the Gaussian product: . Combine and expand the coefficients
to get

Let and substitute in the exponent to get

The first three terms inside the second bracket factor to , and the last two can be reduced to . The original Gaussian product is thus

Here is a verification.

If one opts to use Mathematica‘s own incomplete gamma function , equation (24) can be defined in closed form.

Here is an example.

This is the result.

Here is a comparison with the quadrature result.

Understandably, direct use of the Gamma function is much faster and yields more accurate results. We thank Paul Abbott for providing the example.


[1] M. Hô and J. M. Hernández-Pérez, “Evaluation of Gaussian Molecular Integrals I,” The Mathematica Journal, 14(3), 2012. doi:10.3888/tmj.14-3.
[2] J. Rys, M. Dupuis, and H. F. King, “Computation of Electron Repulsion Integrals Using the Rys Quadrature Method,” Journal of Computational Chemistry, 4(2), 1983 pp. 154-157. doi:10.1002/jcc.540040206.
[3] G. H. F. Diercksen and G. G. Hall, “Intelligent Software: The OpenMol Program,” Computers in Physics, 8(2), 1994 pp. 215-222. doi:10.1063/1.168520.
[4] J. Pérez-Jorda and E. San-Fabián, “A Simple, Efficient and More Reliable Scheme for Automatic Numerical Integration,” Computer Physics Communications, 77(1), 1993 pp. 46-56. doi:10.1016/0010-4655(93)90035-B.
[5] M. Hô and J. M. Hernández-Pérez, “Evaluation of Gaussian Molecular Integrals II,” The Mathematica Journal, 15(1), 2013. doi:10.3888/tmj.15-1.
M. Hô and J. M. Hernández-Pérez, “Evaluation of Gaussian Molecular Integrals,” The Mathematica Journal, 2014. dx.doi.org/doi:10.3888/tmj.16-9.

About the Authors

Minhhuy Hô received his Ph.D. in theoretical chemistry at Queen’s University, Kingston, Ontario, Canada, in 1998. He is currently a professor at the Centro de Investigaciones Químicas at the Universidad Autónoma del Estado de Morelos in Cuernavaca, Morelos, México.

Julio-Manuel Hernández-Pérez obtained his Ph.D. at the Universidad Autónoma del Estado de Morelos in 2008. He has been a professor of chemistry at the Facultad de Ciencias Químicas at the Benemérita Universidad Autónoma de Puebla since 2010.

Minhhuy Hô
Universidad Autónoma del Estado de Morelos
Centro de Investigaciones Químicas
Ave. Universidad, No. 1001, Col. Chamilpa
Cuernavaca, Morelos, Mexico CP 92010

Julio-Manuel Hernández-Pérez

Benemérita Universidad Autónoma de Puebla
Facultad de Ciencias Químicas
Ciudad Universitaria, Col. San Manuel
Puebla, Puebla, Mexico CP 72570