Gradient of integrals with Obara Saika recurrence scheme

Question

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

Answer 1

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

Answer 2

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.