1
0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-08-17 02:41:43 +02:00

Optimized polynomials

This commit is contained in:
Anthony Scemama 2020-10-27 17:35:49 +01:00
parent 1db7d327e7
commit abc0e6b5ab

View File

@ -39,7 +39,23 @@ MunitResult test_qmckl_ao() {
* Polynomials
\[ P_l(\mathbf{r},\mathbf{R}_i) = (x-X_i)^a (y-Y_i)^b (z-Z_i)^c \]
\[
P_l(\mathbf{r},\mathbf{R}_i) = (x-X_i)^a (y-Y_i)^b (z-Z_i)^c
\]
\begin{eqnarray*}
\frac{\partial }{\partial x} P_l\left(\mathbf{r},\mathbf{R}_i \right) & = & a (x-X_i)^{a-1} (y-Y_i)^b (z-Z_i)^c \\
\frac{\partial }{\partial y} P_l\left(\mathbf{r},\mathbf{R}_i \right) & = & b (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c \\
\frac{\partial }{\partial z} P_l\left(\mathbf{r},\mathbf{R}_i \right) & = & c (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1} \\
\end{eqnarray*}
\begin{eqnarray*}
\left( \frac{\partial }{\partial x^2} +
\frac{\partial }{\partial y^2} +
\frac{\partial }{\partial z^2} \right) P_l
\left(\mathbf{r},\mathbf{R}_i \right) & = &
a(a-1) (x-X_i)^{a-2} (y-Y_i)^b (z-Z_i)^c + \\
&& b(b-1) (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c + \\
&& c(c-1) (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1}
\end{eqnarray*}
** =qmckl_ao_powers=
@ -244,6 +260,7 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL,
integer :: lmax_array(3)
real*8 :: pows(-2:lmax,3)
integer, external :: qmckl_ao_powers_f
double precision :: xy, yz, xz
info = 0
@ -293,16 +310,24 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL,
l(2,n) = b
l(3,n) = c
vgl(1,n) = pows(a,1) * pows(b,2) * pows(c,3)
xy = pows(a,1) * pows(b,2)
yz = pows(b,2) * pows(c,3)
xz = pows(a,1) * pows(c,3)
vgl(2,n) = dble(a) * pows(a-1,1) * pows(b ,2) * pows(c ,3)
vgl(3,n) = dble(b) * pows(a ,1) * pows(b-1,2) * pows(c ,3)
vgl(4,n) = dble(c) * pows(a ,1) * pows(b ,2) * pows(c-1,3)
vgl(1,n) = xy * pows(c,3)
xy = dble(c) * xy
xz = dble(b) * xz
yz = dble(a) * yz
vgl(2,n) = pows(a-1,1) * yz
vgl(3,n) = pows(b-1,2) * xz
vgl(4,n) = pows(c-1,3) * xy
vgl(5,n) = &
dble(a) * dble(a-1) * pows(a-2,1) * pows(b ,2) * pows(c ,3) + &
dble(b) * dble(b-1) * pows(a ,1) * pows(b-2,2) * pows(c ,3) + &
dble(c) * dble(c-1) * pows(a ,1) * pows(b ,2) * pows(c-2,3)
dble(a-1) * pows(a-2,1) * yz + &
dble(b-1) * pows(b-2,2) * xz + &
dble(c-1) * pows(c-2,3) * xy
exit
end if
end do