# Gradient of integrals with Obara Saika recurrence scheme

Update I made a mistake in the first formula, now everything is correct as far as I can tell.

I am calculating integrals with the Obara–Saika recurrence scheme[1]. Now I want to calculate also the gradients of each integral. I suppose it is calculated like in the equation below which I derived with the product rule.

$$\hat\partial_i ( a || b ) = - 2\zeta_a ( a + {\bf 1}_i || b ) - 2\zeta_b (a ||b + {\bf 1}_i) + \sum^{l_a}_{j=1} \delta_{ij} ( a - {\bf 1}_j || b ) + \sum^{l_b}_{k=1} \delta_{ik} ( a || b - {\bf 1}_k )$$

I use the same nomenclature as Obara and Saika did, since the paper is somewhat loaded with math I will explain those here again. $\hat\partial_i$ is the derviative in $i$-direction, $a$ and $b$ are cartesian gaussian functions, $\zeta_a$ is the orbital exponent of the gaussian function $a$, $\zeta$ is a shortcut for $\zeta_a + \zeta_b$, $\delta_{ij}$ is the kronecker delta and ${\bf 1}_i$ is an increment or decrement of the polynom of the gaussian function in $i$ direction—that basicly means you get an $d_{ij}$ from $p_i+{\bf 1}_j$. Oh, I almost forgot that one, $(a||b)$ denotes an integral over two gaussian functions.

Update or What I have done so far:

I've done some calculations with the Obara–Saika recurrsion so far, that suggest that the expression I derived for the derivative is correct. If you calculate the derivative explicitly as Erik Kjellgren had done in his answer you get the expression I have shown above. It is a bit cryptic and but if you take a closer look you will see it. Please note that we both have tackled a different problem, I have calculated the derivative in x-direction, while he has calculated the derivative for the movement of the nucleus A in x direction. So mabye I am dealing with the wrong problem.

Let's solve both starting with my problem. The gauss integral as used in the paper from Obara and Saika is given as

$$\varphi({\bf r}; \zeta, {\bf n}, {\bf R}) = (x-R_x)^{n_x}(y-R_y)^{n_y}(z-R_z)^{n_z}\times\exp(-\zeta({\bf r-R})^2)$$

this is nearly the same as used by Erik, but I want to stay in this whole post in the same nomenclature. The next step is to calculate the derivative

$$\hat\partial_i ( a || b ) = \hat\partial_i \int{\rm d}{\bf r} N(\zeta_a,{\bf a})\varphi({\bf r}; \zeta_a, {\bf a}, {\bf A}) N(\zeta_b,{\bf b})\varphi({\bf r}; \zeta_b, {\bf b}, {\bf B}) \\ = N(\zeta_a,{\bf a})N(\zeta_b,{\bf b}) \hat\partial_i \int{\rm d}{\bf r} (x-A_x)^{a_x}(y-A_y)^{a_y}(z-A_z)^{a_z}\times\exp(-\zeta_a({\bf r-A})^2) \\\times (x-B_x)^{b_x}(y-B_y)^{b_y}(z-B_z)^{b_z}\times\exp(-\zeta_b({\bf r-B})^2) \\ = N(\zeta_a,{\bf a})N(\zeta_b,{\bf b}) \int{\rm d}{\bf r}\biggl( \hat\partial_i (x-A_x)^{a_x}(y-A_y)^{a_y}(z-A_z)^{a_z}\times\exp(-\zeta_a({\bf r-A})^2) \\+ (x-A_x)^{a_x}(y-A_y)^{a_y}(z-A_z)^{a_z}\times\hat\partial_i \exp(-\zeta_a({\bf r-A})^2) \\+ \hat\partial_i (x-B_x)^{b_x}(y-B_y)^{b_y}(z-B_z)^{b_z}\times\exp(-\zeta_b({\bf r-B})^2) \\+ (x-B_x)^{b_x}(y-B_y)^{b_y}(z-B_z)^{b_z}\times\hat\partial_i \exp(-\zeta_b({\bf r-B})^2)\biggr) \\ = - 2\zeta_a ( a + {\bf 1}_i || b ) - 2\zeta_b (a ||b + {\bf 1}_i) + \sum^{l_a}_{j=1} \delta_{ij} ( a - {\bf 1}_j || b ) + \sum^{l_b}_{k=1} \delta_{ik} ( a || b - {\bf 1}_k )$$

Okay this was straigth forward. Because I don't know anything about the derivative I applied kronecker-deltas. This should be okay because you can also write something like that for the above expression:

$$\hat\partial_i ( a || b ) = - 2\zeta_a ( a + {\bf 1}_i || b ) - 2\zeta_b (a ||b + {\bf 1}_i) + \delta_{ix} a_x ( a - {\bf 1}_x || b ) + \delta_{iy} a_y ( a - {\bf 1}_y || b ) + \delta_{iz} a_z ( a - {\bf 1}_z || b ) + \delta_{ix} b_x ( a || b - {\bf 1}_x ) + \delta_{iy} b_y ( a || b - {\bf 1}_y ) + \delta_{iz} b_z ( a || b - {\bf 1}_z )$$

I am now at the same point as Erik Kjellgren, nothing new so far, but maybe a bit more general. Now I am going to get rid of the integrals with the Obara–Saika recurrsion, because that is what it is all about.

$$\hat\partial_i ( a || b ) = - 2\zeta_a ( a + {\bf 1}_i || b ) - 2\zeta_b (a ||b + {\bf 1}_i) + \sum^{l_a}_{j=1} \delta_{ij} ( a - {\bf 1}_j || b ) + \sum^{l_b}_{k=1} \delta_{ik} ( a || b - {\bf 1}_k ) \\ = -2\zeta_a \Bigl( (P_i-A_i)(a||b) + \sum^{l_a}_{j=1} \delta_{ij}\frac 1{2\zeta} ( a - {\bf 1}_j || b ) + \sum^{l_b}_{k=1} \delta_{ik}\frac 1{2\zeta} ( a || b - {\bf 1}_k ) \Bigr)\\ -2\zeta_b \Bigl( (P_i-B_i)(a||b) + \sum^{l_a}_{j=1} \delta_{ij}\frac 1{2\zeta} ( a - {\bf 1}_j || b ) + \sum^{l_b}_{k=1} \delta_{ik}\frac 1{2\zeta} ( a || b - {\bf 1}_k ) \Bigr)\\ + \sum^{l_a}_{j=1} \delta_{ij}( a - {\bf 1}_j || b ) + \sum^{l_b}_{k=1} \delta_{ik}( a || b - {\bf 1}_k )$$

After applying the basic recurrence scheme I got something that looks interesting and promising

$$\hat\partial_i ( a || b ) = -2\zeta_a (P_i-A_i)(a||b)-2\zeta_b (P_i-B_i)(a||b) \\ - \sum^{l_a}_{j=1} \delta_{ij} ( a - {\bf 1}_j || b ) - \sum^{l_b}_{k=1} \delta_{ik} ( a || b - {\bf 1}_k ) \\ + \sum^{l_a}_{j=1} \delta_{ij}( a - {\bf 1}_j || b ) + \sum^{l_b}_{k=1} \delta_{ik}( a || b - {\bf 1}_k )$$

The rest is simple math and so I derived the following expression

$$\hat\partial_i (a||b) = -\bigl((2\zeta_a({\rm P}_i - {\rm A}_i) + 2\zeta_b({\rm P}_i - {\rm B}_i)\bigr) (a||b)$$

where $\bf A$ and $\bf B$ are the centers of the gaussian functions $a$ and $b$ and $\bf P$ is the center of the product of both gaussian functions.

Orginal question

Have any of you experience with gradients of molecular integrals and can give me a hint on that issue? Is this the right way to calculate the gradient?

[1] S. Obara and A. Saika, J. Chem. Phys. 84 (7), 1986

• If you could put your derivation from (a||b) all the way to the final equation, into your question that would be awesome – Erik Kjellgren Aug 11 '17 at 19:17
• I believe eqs. 41 and 42 as well as chapter 14 of this review have the answer to your question. It seems like a lot of work to actually understand everything that's going on and translate notations so maybe someone braver than me will do it. – jheindel Aug 12 '17 at 5:25
• @ErikKjellgren I added my derivation for both problems, now my question is a wall of math. I tried to get away from the Obara–Saika notation, but it would be more formulas. Hope it is helpful. – awvwgk Aug 12 '17 at 22:35
• @awvwgk I naivly thought that derivative in the x direction implied with the derivaitve with respect to the geometry in the x direction, but can see I missinterpreted that. When you write "the basic recurrence scheme", what relations are that excactly? When I think about the OS relations for overlap integrals I think of these – Erik Kjellgren Aug 12 '17 at 23:06
• @awvwgk When I look at your final expression I got a feeling something is not right. $2\zeta (P_i - A_i)$ is "just" a constant. So it can be written as: d/dA(a|b)=C*(a|b). Since a and b are Gaussian type functions, this implies that d/dA(a) = C*(a), which we know is not true. – Erik Kjellgren Aug 12 '17 at 23:23

## 2 Answers

I see you have not gotten any answers for some time. I have not myself become familiar enough with the notation used in S. Obara and A. Saika, J. Chem. Phys. 84 (7), 1986. But I have a suggestion for how you can help your self proceed. If you use your derived formula to calculate some derivative, and you try to calculate the same derivatives with finite difference, then you can check and see if you get the same result. If the results are identical to some nummerical error, your equation are proberly correct.

UPDATE:

I have tried to make a derivation of the derivative myself. Here is what I have come up with: The overlap integral is given as:

$$S_{\mu,\nu}=\left(\mu|\nu\right)$$

Here $\mu$ and $\nu$ represent Gaussian functions of the type:

$$\phi_{l,m,n}=N_{l,m,n}x_{a}^{l}y_{a}^{m}z_{a}^{n}\exp\left(-\alpha\left|r-R_{A}\right|^{2}\right)$$

Now the derivative with repsect to the movement of nucleus $A$ can be calculated by using the chain rule:

$$\left.\frac{\partial S_{\mu,\nu}}{\partial X_{A}}\right|_{X_{\mu}\neq X_{\nu}}=\left.\left(\left.\frac{\partial\mu}{\partial X_{A}}\right|\nu\right)\right|_{X_{\mu}=X_{A}}+\left.\left(\mu\left|\frac{\partial\nu}{\partial X_{A}}\right.\right)\right|_{X_{\nu}=X_{A}}$$

Note that $X_{\mu}\neq X_{\nu}$, this is imposed by the invariance of moving all centers [1]:

$$\frac{\partial S_{\mu,\nu}}{\partial X_{A}}+\frac{\partial S_{\mu,\nu}}{\partial X_{B}}=0$$

Here $A=B$. Also a wavefunction must have a depends on $A$ for a derivative to give a value, thus for atom centered wavefunctions, the wavefunction must be placed on the atom for which the derivative is taken. Now the the derivative of the Gaussian functions is given as [2]:

$$\frac{\partial\phi_{l,m,n}}{\partial X_{A}}=\sqrt{(2l+1)\alpha}\phi_{l+1,m,n}-\left.2l\sqrt{\left(\frac{\alpha}{2l-1}\right)}\phi_{l-1,m,n}\right|_{l>0}$$

It should be noted that the factors in front of the Gaussians is because of normalization. Thus now giving:

$$\left.\frac{\partial S_{\mu,\nu}}{\partial X_{A}}\right|_{X_{A}=X_{\mu}\neq X_{\nu}}=\left(\left.\frac{\partial\mu_{l,m,n}}{\partial X_{A}}\right|\nu_{l,m,n}\right)$$

$$=\left(\left.\sqrt{(2l+1)\alpha}\mu_{l+1,m,n}-\left.2l\sqrt{\left(\frac{\alpha}{2l-1}\right)}\mu_{l-1,m,n}\right|_{l>0}\right|\nu_{l,m,n}\right)$$

$$=\sqrt{(2l+1)\alpha}\left(\left.\mu_{l+1,m,n}\right|\nu_{l,m,n}\right)-2l\sqrt{\left(\frac{\alpha}{2l-1}\right)}\left(\left.\left.\mu_{l-1,m,n}\right|_{l>0}\right|\nu_{l,m,n}\right)$$

Now $S_{i,,j}=\left(\mu_{l=i,m,n}|\nu_{l=j,m,n}\right)$and by insertion into the found expression for the derivative:

$$\left.\frac{\partial S_{\mu,\nu}}{\partial X_{A}}\right|_{X_{A}=X_{\mu}\neq X_{\nu}}=\sqrt{(2l+1)\alpha}S_{i+1,j}-\left.2l\sqrt{\left(\frac{\alpha}{2l-1}\right)}S_{i-1.j}\right|_{i>0}$$

Now from here on I was unable to find the equation you found. Also my notation is abit different. I thought you might want it as an alternative way to find the derivative (or maybe you will be able see from this if your equation is right or wrong). Again if you do the calculations with numbers you can see if the two different equations give the same results. Now the Obara-Saika recurrence relation can be considered [3]:

$$S_{i+1,j}=X_{PA}S_{i,j}+\frac{1}{2p}(iS_{i-1,j}+jS_{i,j-1})$$

This can be insterted into the above derivative to eliminate the $i+1$ term:

$$\left.\frac{\partial S_{\mu,\nu}}{\partial X_{A}}\right|_{X_{A}=X_{\mu}\neq X_{\nu}}=\sqrt{(2l+1)\alpha}\left( X_{PA}S_{i,j}+\frac{1}{2p}(iS_{i-1,j}+jS_{i,j-1}) \right)-\left.2l\sqrt{\left(\frac{\alpha}{2l-1}\right)}S_{i-1.j}\right|_{i>0}$$

[1] Trygve Helgaker, Poul Jorgensen, Jeppe Olsen; Molecular Electronic-Structure Theory. Equation: 9.3.4

[2] Attila Szabo, Neil S. Ostlund;Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory. page. 442

[3] Trygve Helgaker, Poul Jorgensen, Jeppe Olsen; Molecular Electronic-Structure Theory. Equation: 9.3.8

• That is clear, this is almost the first formula in my question, now in the correct form. I’m more interested in the application of the Obara–Saika-recursion on the $S_{i+1,j}$ part. – awvwgk Aug 12 '17 at 22:40
• Thanks for your input. I was looking at the wrong problem all the time. – awvwgk Aug 13 '17 at 0:06
• @awvwgk I have tried to insert the OS relations into the derivative. Now my derivation, was not very general at all. Also because of my ignorence (imfamilarity) with your notation, it is a bit hard to compare the results. To complicate even more, I have normalized the expression. But I think our two expressions look alike very much. – Erik Kjellgren Aug 13 '17 at 0:07
• They are the same. – awvwgk Aug 13 '17 at 0:09

I moved my extended explanation from my question to a real answer. The given problem was to calculate the gradient of overlap integrals. Thanks to Erik Kjellgren I noticed that I am looking at the false problem.

So far I have calculated derivatives for cartesian coordinates. This is possible but not very helpful in case I want to have gradient for geometry optimizations. I can use a bit from my previous work with

$$\frac{\partial}{\partial R_i}\varphi ({\bf r}; \zeta, {\bf n}, {\bf R}) = - \frac{\partial}{\partial r_i}\varphi ({\bf r}; \zeta, {\bf n}, {\bf R})$$

and get

$$\frac{\partial}{\partial A_i} (a || b) = 2\zeta_a ( a + {\bf 1}_i || b ) - \sum^{l_a}_{j=1} \delta_{ij}( a - {\bf 1}_j || b ).$$

This is the exactly same as in Erik Kjellgren’s answer. If I use the Obara–Saika recurrsion formula on the expression above I get the following formula.

$$\frac{\partial}{\partial A_i} (a || b) = 2\zeta_a \Bigl( (P_i-A_i)(a||b) + \sum^{l_a}_{j=1} \delta_{ij}\frac 1{2\zeta}( a - {\bf 1}_j || b ) + \sum^{l_b}_{k=1} \delta_{ik}\frac 1{2\zeta}( a || b - {\bf 1}_k ) \Bigr) - \sum^{l_a}_{j=1} \delta_{ij}( a - {\bf 1}_j || b )$$

From here on I cannot go any further. The equation I have derived now does not look easy anymore, but I think it is fine, all intermediates are available from the calculation of the overlap integral itself, so I don’t need any additional expensive calculations.