From 2820978da637e080dd244b561fb851dec1bd67ee Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Mon, 4 Nov 2024 16:44:28 +0100
Subject: [PATCH] Fixed bug in Jastrow

---
 org/qmckl_ao.org            | 11 +++++------
 org/qmckl_jastrow_champ.org |  4 ++--
 2 files changed, 7 insertions(+), 8 deletions(-)

diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org
index 284daa1..dd57585 100644
--- a/org/qmckl_ao.org
+++ b/org/qmckl_ao.org
@@ -3943,7 +3943,6 @@ function qmckl_compute_ao_basis_shell_gaussian_vgl( &
 
      do inucl=1,nucl_num
 
-        ! The shift below avoids having an exact zero on a node of the orbital
         x = coord(ipoint,1) - nucl_coord(inucl,1)
         y = coord(ipoint,2) - nucl_coord(inucl,2)
         z = coord(ipoint,3) - nucl_coord(inucl,3)
@@ -4572,7 +4571,7 @@ function qmckl_ao_polynomial_vgl_doc (context,    &
 
   ! The shift below is such that polynomials will not make the AO equal to zero at the nodes of the orbitals
   do i=1,3
-     Y(i) = X(i) - R(i) + 1.d-20
+     Y(i) = (X(i) - R(i)) + 1.d-20
   end do
 
   if (lmax == 0) then
@@ -4967,7 +4966,7 @@ qmckl_ao_polynomial_transp_vgl_hpc_inline (const qmckl_context context,
 
       if (ldv == 4) {
 
-        const double tmp[20] = {1.0, X[0]-R[0]+shift, X[1]-R[1]+shift, X[2]-R[2]+shift,
+        const double tmp[20] = {1.0, (X[0]-R[0])+shift, (X[1]-R[1])+shift, (X[2]-R[2])+shift,
                                 0.0, 1.0, 0.0, 0.0,
                                 0.0, 0.0, 1.0, 0.0,
                                 0.0, 0.0, 0.0, 1.0,
@@ -4991,9 +4990,9 @@ qmckl_ao_polynomial_transp_vgl_hpc_inline (const qmckl_context context,
           vgl5[k] = 0.0;
         }
         vgl1[0] = 1.0;
-        vgl1[1] = X[0]-R[0]+shift;
-        vgl1[2] = X[1]-R[1]+shift;
-        vgl1[3] = X[2]-R[2]+shift;
+        vgl1[1] = (X[0]-R[0])+shift;
+        vgl1[2] = (X[1]-R[1])+shift;
+        vgl1[3] = (X[2]-R[2])+shift;
         vgl2[1] = 1.0;
         vgl3[2] = 1.0;
         vgl4[3] = 1.0;
diff --git a/org/qmckl_jastrow_champ.org b/org/qmckl_jastrow_champ.org
index 508f7ff..76c5396 100644
--- a/org/qmckl_jastrow_champ.org
+++ b/org/qmckl_jastrow_champ.org
@@ -2999,7 +2999,6 @@ qmckl_compute_jastrow_champ_factor_ee_gl_hpc(const qmckl_context context,
           factor_ee_gl_2[i] = factor_ee_gl_2[i] + f*dx[2];
           factor_ee_gl_3[i] = factor_ee_gl_3[i] + f*dx[3];
         } else {
-          touched = true;
           factor_ee_gl_0[i] = f*dx[0];
           factor_ee_gl_1[i] = f*dx[1];
           factor_ee_gl_2[i] = f*dx[2];
@@ -3022,6 +3021,7 @@ qmckl_compute_jastrow_champ_factor_ee_gl_hpc(const qmckl_context context,
           factor_ee_gl_3[i] = factor_ee_gl_3[i] + f1*kf[k-1]*grad_c2;
         }
       }
+      touched = true;
     }
     if (!touched) {
       memset(&(factor_ee_gl[nw*4*elec_num]), 0, elec_num*4*sizeof(double));
@@ -4902,7 +4902,6 @@ qmckl_compute_jastrow_champ_factor_en_gl_hpc (const qmckl_context context,
           factor_en_gl_2[i] = factor_en_gl_2[i] + f*dx[2];
           factor_en_gl_3[i] = factor_en_gl_3[i] + f*dx[3];
         } else {
-          touched = true;
           factor_en_gl_0[i] = f*dx[0];
           factor_en_gl_1[i] = f*dx[1];
           factor_en_gl_2[i] = f*dx[2];
@@ -4927,6 +4926,7 @@ qmckl_compute_jastrow_champ_factor_en_gl_hpc (const qmckl_context context,
           factor_en_gl_3[i] = factor_en_gl_3[i] + f1*kf[k-1]*grad_c2;
         }
       }
+      touched = true;
     }
     if (!touched) {
       memset(&(factor_en_gl[nw*4*elec_num]), 0, elec_num*4*sizeof(double));