Avoid FPE in Jastrow GL

org/qmckl_ao.org

@@ -5483,7 +5483,9 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
      if (L(3,j) > 1) then
         w = w + L(3,j) * (L(3,j)-1) * Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**(L(3,j)-2)
      end if
-     if (dabs(1.d0 - VGL(5,j) / w) > epsilon ) return
+     if (w /= 0.d0) then
+       if (dabs(1.d0 - VGL(5,j) / w) > epsilon ) return
+     endif
   end do
 
   test_qmckl_ao_polynomial_vgl = QMCKL_SUCCESS

org/qmckl_distance.org

@@ -86,7 +86,7 @@ int main() {
           const double* B,
           const int64_t ldb,
           double* const C,
-          const int64_t ldc ); 
+          const int64_t ldc );
     #+end_src
 
     #+begin_src f90 :tangle (eval f)
@@ -485,7 +485,7 @@ end function test_qmckl_distance_sq
           const double* B,
           const int64_t ldb,
           double* const C,
-          const int64_t ldc ); 
+          const int64_t ldc );
     #+end_src
 
 *** Source
@@ -909,7 +909,7 @@ end function test_qmckl_dist
           const int64_t ldb,
           double* const C,
           const int64_t ldc,
-          const double rescale_factor_kappa ); 
+          const double rescale_factor_kappa );
     #+end_src
 
 *** Source
@@ -1215,7 +1215,7 @@ integer(qmckl_exit_code) function test_qmckl_dist_rescaled(context) bind(C)
 
   integer*8, parameter :: elec_num = 10_8
   integer*8, parameter :: nucl_num = 2_8
-  
+
   double precision :: nucl_coord(nucl_num,3) = reshape( (/ &
         0.0d0,  0.0d0 ,  &
         0.0d0,  0.0d0 ,  &
@@ -1263,19 +1263,19 @@ integer(qmckl_exit_code) function test_qmckl_dist_rescaled(context) bind(C)
        1.20919677d0, 1.3111856d0, 1.26122875d0, 0.22122563d0, 0.55669168d0 /), shape(ref_en) )
 
   double precision, allocatable :: distance_ee(:,:), distance_en(:,:)
-  
+
   allocate( distance_ee(elec_num,elec_num), distance_en(elec_num,nucl_num) )
-  
+
   print *, 'ee'
   test_qmckl_dist_rescaled = &
        qmckl_distance_rescaled(context, 'N', 'N', elec_num, elec_num, elec_coord, &
        size(elec_coord,1)*1_8, elec_coord, size(elec_coord,1)*1_8, &
        distance_ee, size(distance_ee,1)*1_8, kappa)
-  
+
   if (test_qmckl_dist_rescaled /= QMCKL_SUCCESS) return
-  
+
   test_qmckl_dist_rescaled = QMCKL_FAILURE
-  
+
   do j=1,elec_num
      do i=1,elec_num
         print *, i,j,real(distance_ee(i,j)), real(ref_ee(i,j))
@@ -1294,7 +1294,7 @@ integer(qmckl_exit_code) function test_qmckl_dist_rescaled(context) bind(C)
   if (test_qmckl_dist_rescaled /= QMCKL_SUCCESS) return
 
   test_qmckl_dist_rescaled = QMCKL_FAILURE
-  
+
   do j=1,nucl_num
      do i=1,elec_num
         print *, i,j,real(distance_en(i,j)), real(ref_en(i,j))
@@ -1385,7 +1385,7 @@ end function test_qmckl_dist_rescaled
     - ~A~ is allocated with at least $3 \times m \times 8$ bytes
     - ~B~ is allocated with at least $3 \times n \times 8$ bytes
     - ~C~ is allocated with at least $4 \times m \times n \times 8$ bytes
-      
+
     #+CALL: generate_c_header(table=qmckl_distance_rescaled_gl_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
 
     #+RESULTS:
@@ -1402,7 +1402,7 @@ end function test_qmckl_dist_rescaled
           const int64_t ldb,
           double* const C,
           const int64_t ldc,
-          const double rescale_factor_kappa ); 
+          const double rescale_factor_kappa );
     #+end_src
 
     #+begin_src f90 :tangle (eval f)
@@ -1504,6 +1504,10 @@ integer function qmckl_distance_rescaled_gl_f(context, transa, transb, m, n, &
            y = A(2,i) - B(2,j)
            z = A(3,i) - B(3,j)
            dist = dsqrt(x*x + y*y + z*z)
+           ! Avoid floating-point exception
+           if (dist == 0.d0) then
+             dist = 69.d0/rescale_factor_kappa
+           endif
            dist_inv = 1.0d0/dist
            rij = dexp(-rescale_factor_kappa * dist)
            C(1,i,j) = x * dist_inv * rij
@@ -1521,6 +1525,10 @@ integer function qmckl_distance_rescaled_gl_f(context, transa, transb, m, n, &
            y = A(i,2) - B(2,j)
            z = A(i,3) - B(3,j)
            dist = dsqrt(x*x + y*y + z*z)
+           ! Avoid floating-point exception
+           if (dist == 0.d0) then
+             dist = 69.d0/rescale_factor_kappa
+           endif
            dist_inv = 1.0d0/dist
            rij = dexp(-rescale_factor_kappa * dist)
            C(1,i,j) = x * dist_inv * rij
@@ -1538,6 +1546,10 @@ integer function qmckl_distance_rescaled_gl_f(context, transa, transb, m, n, &
            y = A(2,i) - B(j,2)
            z = A(3,i) - B(j,3)
            dist = dsqrt(x*x + y*y + z*z)
+           ! Avoid floating-point exception
+           if (dist == 0.d0) then
+             dist = 69.d0/rescale_factor_kappa
+           endif
            dist_inv = 1.0d0/dist
            rij = dexp(-rescale_factor_kappa * dist)
            C(1,i,j) = x * dist_inv * rij
@@ -1555,6 +1567,10 @@ integer function qmckl_distance_rescaled_gl_f(context, transa, transb, m, n, &
            y = A(i,2) - B(j,2)
            z = A(i,3) - B(j,3)
            dist = dsqrt(x*x + y*y + z*z)
+           ! Avoid floating-point exception
+           if (dist == 0.d0) then
+             dist = 69.d0/rescale_factor_kappa
+           endif
            dist_inv = 1.0d0/dist
            rij = dexp(-rescale_factor_kappa * dist)
            C(1,i,j) = x * dist_inv * rij
@@ -1570,7 +1586,7 @@ end function qmckl_distance_rescaled_gl_f
     #+end_src
 
     This function is more efficient when ~A~ and ~B~ are transposed.
-    
+
     #+CALL: generate_c_interface(table=qmckl_distance_rescaled_gl_args,fname=get_value("Name"))
 
    #+RESULTS:

org/qmckl_jastrow_champ.org

@@ -6095,14 +6095,23 @@ integer function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_f( &
   een_rescaled_e_gl     = 0.0d0
   do nw = 1, walk_num
     do j = 1, elec_num
-      do i = 1, elec_num
+      do i = 1, j-1
         rij_inv = 1.0d0 / ee_distance(i, j, nw)
         do ii = 1, 3
           elec_dist_gl(ii, i, j) = (coord_ee(i, ii, nw) - coord_ee(j, ii, nw)) * rij_inv
         end do
         elec_dist_gl(4, i, j) = 2.0d0 * rij_inv
       end do
+
       elec_dist_gl(:, j, j) = 0.0d0
+
+      do i = j+1, elec_num
+        rij_inv = 1.0d0 / ee_distance(i, j, nw)
+        do ii = 1, 3
+          elec_dist_gl(ii, i, j) = (coord_ee(i, ii, nw) - coord_ee(j, ii, nw)) * rij_inv
+        end do
+        elec_dist_gl(4, i, j) = 2.0d0 * rij_inv
+      end do
     end do
 
     ! prepare the actual een table