From abc0e6b5ab2a0c0bc5112426c44320cfee2e31c8 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue, 27 Oct 2020 17:35:49 +0100
Subject: [PATCH] Optimized polynomials

---
 src/qmckl_ao.org | 43 ++++++++++++++++++++++++++++++++++---------
 1 file changed, 34 insertions(+), 9 deletions(-)

diff --git a/src/qmckl_ao.org b/src/qmckl_ao.org
index ee2c413..9f2399b 100644
--- a/src/qmckl_ao.org
+++ b/src/qmckl_ao.org
@@ -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
@@ -443,7 +468,7 @@ munit_assert_int(0, ==, test_qmckl_ao_polynomial_vgl(context));
     #+END_SRC
     #+END_SRC
    
-
+    
 * TODO Gaussian basis functions
 
 * TODO Slater basis functions