I wanna to evaluate molecular integrals (overlap, kinetic and etc) using Gaussian basis functions. To do so, there are two approaches: contracted and uncontracted basis functions in which the former means the primitive Gaussian functions ($G_1$ and $G_2$) are summed to form contracted ones and finally by summing the contracted functions ($\phi_1$ and $\phi_2$) using appropriate coefficients, the final molecular orbital is formed, for two contracted basis functions and three primitives we have:

$$ \psi=c_1\phi_1+c_2\phi_2 $$ $$ \phi_1=d_{11}G_{11}+d_{12}G_{12}+d_{13}G_{13} $$ $$ \phi_2=d_{21}G_{21}+d_{22}G_{22}+d_{23}G_{23} $$ $$ G_{11}=e^{-\alpha_{11}r^2} $$ where the exponents are given as a 2x3 matrix (generally exponents matrix in contraction basis functions approach is a [number of electrons]x[number of primitive Gaussians] matrix and for uncontracted functions approache is a [number of electrons]x[number of basis functions] matrix). In this way the final overlap matrix for example is

$$ \begin{bmatrix} c_1c_1\phi_1\phi_1 && c_1c_2\phi_1\phi_2\\ c_2c_1\phi_2\phi_1 && c_2c_2\phi_2\phi_2\\ \end{bmatrix} $$

For simulating total wave function using two uncontracted basis functions, we have

$$ \psi=c_1\phi_1+c_2\phi_2 $$ where it's not necessary to use any summation to obtain $\phi_1$ and $\phi_2$, instead we can construct them directly by knowing Gaussian exponents:

$$ \phi_1=e^{-\alpha_{11}r^2} $$

Again the overlap matrix has a form like above matrix.

Are the mentioned explanations true? How should we apply these differences in a Mathematica code?


1 Answer 1


A lot of stuff is going on in this question. So for clarity let's make this answer very explicit. We can try to consider the two cases of using an uncontracted basis and a contracted basis to clarify the differences. As a case let's look at STO-$2$G for $\mathrm{H_2}$. In the STO-$2$G basis set $\mathrm{H_2}$ would normally be described by two basis functions each consisting of two primitives.

Uncontracted case

If we use the STO-$2$G basis set in the case of using uncontracted Gaussian basis functions, each hydrogen will be described by two basis functions. Our overlap matrix will therefore be of the form:

$$ S= \begin{bmatrix} \left<\phi_1^\mathrm{H_1}\left|\phi_1^\mathrm{H_1}\right.\right> & \left<\phi_1^\mathrm{H_1}\left|\phi_2^\mathrm{H_1}\right.\right> & \left<\phi_1^\mathrm{H_1}\left|\phi_3^\mathrm{H_2}\right.\right> & \left<\phi_1^\mathrm{H_1}\left|\phi_4^\mathrm{H_2}\right.\right> \\ \left<\phi_2^\mathrm{H_1}\left|\phi_1^\mathrm{H_1}\right.\right> & \left<\phi_2^\mathrm{H_1}\left|\phi_2^\mathrm{H_1}\right.\right> & \left<\phi_2^\mathrm{H_1}\left|\phi_3^\mathrm{H_2}\right.\right> & \left<\phi_2^\mathrm{H_1}\left|\phi_4^\mathrm{H_2}\right.\right> \\ \left<\phi_3^\mathrm{H_2}\left|\phi_1^\mathrm{H_1}\right.\right> & \left<\phi_3^\mathrm{H_2}\left|\phi_2^\mathrm{H_1}\right.\right> & \left<\phi_3^\mathrm{H_2}\left|\phi_3^\mathrm{H_2}\right.\right> & \left<\phi_3^\mathrm{H_2}\left|\phi_4^\mathrm{H_2}\right.\right> \\ \left<\phi_4^\mathrm{H_2}\left|\phi_1^\mathrm{H_1}\right.\right> & \left<\phi_4^\mathrm{H_2}\left|\phi_2^\mathrm{H_1}\right.\right> & \left<\phi_4^\mathrm{H_2}\left|\phi_3^\mathrm{H_2}\right.\right> & \left<\phi_4^\mathrm{H_2}\left|\phi_4^\mathrm{H_2}\right.\right> \end{bmatrix} $$

In this notation we have that:

$$\left<\phi_i^\mathrm{H_k}\left|\phi_j^\mathrm{H_l}\right.\right> = N_iN_j\int_{-\infty}^\infty P_i\exp\left(-\alpha_ir_i^2\right)P_j\exp\left(-\alpha_jr_j^2\right) d\tau$$

With $P$ being a function of the angular momentum, in cartesian coordinates:


It can be seen that $P=1$ for s-functions. As we can see we have 4 basis functions in total (two on each hydrogen). You wrote something about, you would expect that in the case of uncontracted that we would have ""[number of electrons]x[number of basis functions] matrix"". As we can see this is not the case. Right now our matrix is in the atomic basis. When we solve the Hartree-Fock equations and find the coefficients that will take us to a molecular basis, we will see that the orbitals that are occupied in the atomic orbital basis, will have a contribution to occupied orbitals in the molecular orbital basis, and are therefore still important.

Contracted case

Lets now take a look at how the STO-$2$G basis set will be used in the case of contracted basis functions. Now our overlap matrix will take the following form:

$$ S= \begin{bmatrix} \left<\phi_1^\mathrm{H_1}\left|\phi_1^\mathrm{H_1}\right.\right> & \left<\phi_1^\mathrm{H_1}\left|\phi_2^\mathrm{H_2}\right.\right> \\ \left<\phi_2^\mathrm{H_2}\left|\phi_1^\mathrm{H_1}\right.\right> & \left<\phi_2^\mathrm{H_2}\left|\phi_2^\mathrm{H_2}\right.\right> \end{bmatrix} $$

At first glance the overlap matrix look just like the one in the uncontracted case, but smaller. What has changed is how the elements are defined:

$$\left<\phi_i^\mathrm{H_k}\left|\phi_j^\mathrm{H_l}\right.\right> = \sum_I\sum_J c_Ic_JN_IN_J\int_{-\infty}^\infty P_I\exp\left(-\alpha_Ir_I^2\right)P_J\exp\left(-\alpha_Jr_J^2\right) d\tau$$

We can see that we now sum over all the primitives included in the contraction. To make it more clear, let's write out how the second element in the overlap matrix looks like. Now remember that each of the basis functions consists of two primitives:

$$ \begin{eqnarray} \left<\phi_1^\mathrm{H_1}\left|\phi_2^\mathrm{H_2}\right.\right> &=& c_{11}c_{21}N_{11}N_{21}\int_{-\infty}^\infty \exp\left(-\alpha_{11}r_{11}^2\right)\exp\left(-\alpha_{21}r_{21}^2\right) d\tau \\ &+&c_{11}c_{22}N_{11}N_{22}\int_{-\infty}^\infty \exp\left(-\alpha_{11}r_{11}^2\right)\exp\left(-\alpha_{22}r_{22}^2\right) d\tau \\ &+&c_{12}c_{21}N_{12}N_{21}\int_{-\infty}^\infty \exp\left(-\alpha_{12}r_{12}^2\right)\exp\left(-\alpha_{21}r_{21}^2\right) d\tau \\ &+&c_{12}c_{22}N_{12}N_{22}\int_{-\infty}^\infty\exp\left(-\alpha_{12}r_{12}^2\right)\exp\left(-\alpha_{22}r_{22}^2\right) d\tau \end{eqnarray} $$

In this notation the first number refers to the basis function and the second refers to the primitive. As we can see each element of the overlap matrix in the contracted basis is made from four integrals. Even though our overlap matrix is four times as small compared to that of the uncontracted basis, we still have to evaluate the same number of primitive integrals. In this case we have investigated here, where we used STO-$2$G, we have no unoccupied orbitals. In general, this will not be the case, most basis sets (also when used in the contracted basis) will have unoccupied orbitals.

Other comments

In your question I can see that you mention molecular orbitals. Just to be clear. When we do the integrals we are in the atomic orbital basis. The molecular orbital basis is the one we find when we solve the Hartree-Fock equations. I.e. the molecular orbitals have nothing to do with uncontracted vs contrancted basis functions, though we will get more molecular orbitals if we use uncontracted basis functions, because we will get the same number of molecular orbitals as we had atomic orbitals.

Regarding how to implement this in Mathematica I have no idea. Though as my understanding of how Mathematica is ought to be used, I would recommend that you would implement your integral code in a language like Python. For inspiration you can take a look at this guide, which implements all the integrals with contracted basis functions.

  • $\begingroup$ Many Thanks for your detailed answer. I think somewhere you misunderstood my mean, what is the exponents matrix in both cases? is that other than a 2x4 matrix? $\endgroup$
    – Wisdom
    Jun 15, 2018 at 3:48
  • $\begingroup$ @NSR I see. How you store the exponents in your code is pretty much free choice. You can store them in a matrix if you want. In the first case (uncontracted), you need to exponents for every integral (one exponent for each Gaussian). In the second case with STO$2$G you need two exponents for the first basis function, and two exponents for the second basis functions, i.e. you will need two vectors of length four. $\endgroup$ Jun 17, 2018 at 13:01
  • $\begingroup$ Yes, however, we can combine two vectors of length four and build a 2x4 matrix again, so finally, we encounter with the same exponent matrices, right? $\endgroup$
    – Wisdom
    Jun 18, 2018 at 4:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.