Working on TC ints. Not well tested

plugins/local/tc_int/jast_grad_full.irp.f

@@ -59,7 +59,7 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz)
   double precision              :: r2(3)
   double precision              :: dx, dy, dz, r12, tmp
   double precision              :: rn(3), f1A, grad1_f1A(3), f2A, grad2_f2A(3), g12, grad1_g12(3)
-  double precision              :: tmp1, tmp2
+  double precision              :: tmp1, tmp2, dist
   integer                       :: powmax1, powmax, powmax2
   double precision, allocatable :: f1A_power(:), f2A_power(:), double_p(:), g12_power(:)
 
@@ -90,30 +90,105 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz)
     gradx(jpoint) = 0.d0
     grady(jpoint) = 0.d0
     gradz(jpoint) = 0.d0
+
+    call jBH_elem_fct_grad_alpha1(r1, r2, g12, grad1_g12)
+
+!    dist =   (r1(1) - r2(1)) * (r1(1) - r2(1)) &
+!           + (r1(2) - r2(2)) * (r1(2) - r2(2)) &
+!           + (r1(3) - r2(3)) * (r1(3) - r2(3))
+!
+!    if(dist .ge. 1d-15) then
+!      dist = dsqrt( dist )
+!
+!      tmp1 = 1.d0 / (1.d0 + dist)
+!
+!      g12 = dist * tmp1
+!      tmp2 = tmp1 * tmp1 / dist
+!      grad1_g12(1) = tmp2 * (r1(1) - r2(1))
+!      grad1_g12(2) = tmp2 * (r1(2) - r2(2))
+!      grad1_g12(3) = tmp2 * (r1(3) - r2(3))
+!
+!    else
+!
+!      grad1_g12(1) = 0.d0
+!      grad1_g12(2) = 0.d0
+!      grad1_g12(3) = 0.d0
+!      g12 = 0.d0
+!
+!    endif
+!
+    do p = 1, powmax2
+      g12_power(p) = g12_power(p-1) * g12
+    enddo
+
     do i_nucl = 1, nucl_num
 
       rn(1) = nucl_coord(i_nucl,1)
       rn(2) = nucl_coord(i_nucl,2)
       rn(3) = nucl_coord(i_nucl,3)
 
-      call jBH_elem_fct_grad(jBH_en(i_nucl), r1, rn, f1A, grad1_f1A)
+        call jBH_elem_fct_grad_alpha1(r1, rn, f1A, grad1_f1A)
-      call jBH_elem_fct_grad(jBH_en(i_nucl), r2, rn, f2A, grad2_f2A)
+!      dist =   (r1(1) - rn(1)) * (r1(1) - rn(1)) &
-      call jBH_elem_fct_grad(jBH_ee(i_nucl), r1, r2, g12, grad1_g12)
+!             + (r1(2) - rn(2)) * (r1(2) - rn(2)) &
+!             + (r1(3) - rn(3)) * (r1(3) - rn(3))
+!      if (dist > 1.d-15) then
+!        dist = dsqrt( dist )
+!
+!        tmp1 = 1.d0 / (1.d0 + dist)
+!
+!        f1A = dist * tmp1
+!        tmp2 = tmp1 * tmp1 / dist
+!        grad1_f1A(1) = tmp2 * (r1(1) - rn(1))
+!        grad1_f1A(2) = tmp2 * (r1(2) - rn(2))
+!        grad1_f1A(3) = tmp2 * (r1(3) - rn(3))
+!
+!      else
+!
+!        grad1_f1A(1) = 0.d0
+!        grad1_f1A(2) = 0.d0
+!        grad1_f1A(3) = 0.d0
+!        f1A = 0.d0
+!
+!      endif
+
+        call jBH_elem_fct_grad_alpha1(r2, rn, f2A, grad2_f2A)
+!      dist =   (r2(1) - rn(1)) * (r2(1) - rn(1)) &
+!             + (r2(2) - rn(2)) * (r2(2) - rn(2)) &
+!             + (r2(3) - rn(3)) * (r2(3) - rn(3))
+!
+!      if (dist > 1.d-15) then
+!        dist = dsqrt( dist )
+!
+!        tmp1 = 1.d0 / (1.d0 + dist)
+!
+!        f2A = dist * tmp1
+!        tmp2 = tmp1 * tmp1 / dist
+!        grad2_f2A(1) = tmp2 * (r2(1) - rn(1))
+!        grad2_f2A(2) = tmp2 * (r2(2) - rn(2))
+!        grad2_f2A(3) = tmp2 * (r2(3) - rn(3))
+!
+!      else
+!
+!        grad2_f2A(1) = 0.d0
+!        grad2_f2A(2) = 0.d0
+!        grad2_f2A(3) = 0.d0
+!        f2A = 0.d0
+!
+!      endif
 
       ! Compute powers of f1A and f2A
       do p = 1, powmax1
         f1A_power(p) = f1A_power(p-1) * f1A
         f2A_power(p) = f2A_power(p-1) * f2A
       enddo
-      do p = 1, powmax2
-        g12_power(p) = g12_power(p-1) * g12
-      enddo
 
       do p = 1, jBH_size
         mpA = jBH_m(p,i_nucl)
         npA = jBH_n(p,i_nucl)
         opA = jBH_o(p,i_nucl)
         tmp = jBH_c(p,i_nucl)
+!        if (dabs(tmp) <= 1.d-10) cycle
+!
         if(mpA .eq. npA) then
           tmp = tmp * 0.5d0
         endif
@@ -132,3 +207,39 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz)
   return
 end
 
+subroutine jBH_elem_fct_grad_alpha1(r1, r2, fct, grad1_fct)
+
+  implicit none
+  double precision, intent(in)  :: r1(3), r2(3)
+  double precision, intent(out) :: fct, grad1_fct(3)
+  double precision              :: dist, tmp1, tmp2
+
+  dist =   (r1(1) - r2(1)) * (r1(1) - r2(1)) &
+         + (r1(2) - r2(2)) * (r1(2) - r2(2)) &
+         + (r1(3) - r2(3)) * (r1(3) - r2(3))
+
+
+  if(dist .ge. 1d-15) then
+    dist = dsqrt( dist )
+
+    tmp1 = 1.d0 / (1.d0 + dist)
+
+    fct = dist * tmp1
+    tmp2 = tmp1 * tmp1 / dist
+    grad1_fct(1) = tmp2 * (r1(1) - r2(1))
+    grad1_fct(2) = tmp2 * (r1(2) - r2(2))
+    grad1_fct(3) = tmp2 * (r1(3) - r2(3))
+
+  else
+
+    grad1_fct(1) = 0.d0
+    grad1_fct(2) = 0.d0
+    grad1_fct(3) = 0.d0
+    fct = 0.d0
+
+  endif
+
+  return
+end
+
+! ---

plugins/local/tc_int/jast_utils_bh.irp.f

@@ -1,31 +1,39 @@
 
 ! ---
 
+
+
 subroutine jBH_elem_fct_grad(alpha, r1, r2, fct, grad1_fct)
 
   implicit none
   double precision, intent(in)  :: alpha, r1(3), r2(3)
   double precision, intent(out) :: fct, grad1_fct(3)
-  double precision              :: dist, tmp1, tmp2
+  double precision              :: dist, tmp1, tmp2, dist_inv
 
-  dist = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) &
+  dist =   (r1(1) - r2(1)) * (r1(1) - r2(1)) &
          + (r1(2) - r2(2)) * (r1(2) - r2(2)) &
-              + (r1(3) - r2(3)) * (r1(3) - r2(3)) )
+         + (r1(3) - r2(3)) * (r1(3) - r2(3))
 
 
-  if(dist .ge. 1d-10) then
+  if(dist .ge. 1d-15) then
+    dist_inv = 1.d0/dsqrt( dist )
+    dist = dist_inv * dist
+
     tmp1 = 1.d0 / (1.d0 + alpha * dist)
 
     fct = alpha * dist * tmp1
-    tmp2 = alpha * tmp1 * tmp1 / dist
+    tmp2 = alpha * tmp1 * tmp1 * dist_inv
     grad1_fct(1) = tmp2 * (r1(1) - r2(1))
     grad1_fct(2) = tmp2 * (r1(2) - r2(2))
     grad1_fct(3) = tmp2 * (r1(3) - r2(3))
+
   else
+
     grad1_fct(1) = 0.d0
     grad1_fct(2) = 0.d0
     grad1_fct(3) = 0.d0
     fct = 0.d0
+
   endif
 
   return