IPP astice: OK

src/ao_many_one_e_ints/grad2_jmu_modif.irp.f

@@ -164,7 +164,7 @@ END_PROVIDER
 
 ! ---
 
-BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final_grid)]
+BEGIN_PROVIDER [double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final_grid)]
 
   BEGIN_DOC
   !

src/non_h_ints_mu/debug_fit.irp.f

@@ -18,13 +18,13 @@ program debug_fit
   PROVIDE mu_erf j1b_pen
 
   !call test_j1b_nucl()
-  call test_grad_j1b_nucl()
+  !call test_grad_j1b_nucl()
   !call test_lapl_j1b_nucl()
 
   !call test_list_b2()
-  !call test_list_b3()
+  call test_list_b3()
 
-  call test_fit_u()
+  !call test_fit_u()
   !call test_fit_u2()
   !call test_fit_ugradu()
 
@@ -236,16 +236,25 @@ subroutine test_list_b3()
   integer                    :: ipoint
   double precision           :: acc_ij, acc_tot, eps_ij, i_exc, i_tmp, i_num, normalz
   double precision           :: r(3)
-  double precision, external :: j1b_nucl
+  double precision           :: grad_num(3), eps_der, eps_lap, tmp_der, tmp_lap, i0, ip, im
+  double precision, external :: j1b_nucl_square
 
   print*, ' test_list_b3 ...'
 
+  eps_ij  = 1d-7
+
+  eps_der = 1d-5
+  tmp_der = 0.5d0 / eps_der
+
+  eps_lap = 1d-4
+  tmp_lap = 1.d0 / (eps_lap*eps_lap)
+
+  ! ---
+
   PROVIDE v_1b_list_b3
 
-  eps_ij  = 1d-7
   acc_tot = 0.d0
   normalz = 0.d0
-
   do ipoint = 1, n_points_final_grid
 
     r(1) = final_grid_points(1,ipoint)
@@ -253,8 +262,7 @@ subroutine test_list_b3()
     r(3) = final_grid_points(3,ipoint)
 
     i_exc  = v_1b_list_b3(ipoint) 
-    i_tmp  = j1b_nucl(r)
+    i_num  = j1b_nucl_square(r)
-    i_num  = i_tmp * i_tmp
     acc_ij = dabs(i_exc - i_num)
     if(acc_ij .gt. eps_ij) then
       print *, ' problem in list_b3 on', ipoint
@@ -267,8 +275,136 @@ subroutine test_list_b3()
     normalz += dabs(i_num)
   enddo
 
-  print*, ' acc_tot = ', acc_tot
+  print*, ' acc_tot on val = ', acc_tot
-  print*, ' normalz = ', normalz
+  print*, ' normalz on val = ', normalz
+
+  ! ---
+
+  PROVIDE v_1b_square_grad
+
+  acc_tot = 0.d0
+  normalz = 0.d0
+  do ipoint = 1, n_points_final_grid
+
+    r(1) = final_grid_points(1,ipoint)
+    r(2) = final_grid_points(2,ipoint)
+    r(3) = final_grid_points(3,ipoint)
+
+    i_exc  = v_1b_square_grad(ipoint,1)
+    r(1)   = r(1) + eps_der
+    ip     = j1b_nucl_square(r)
+    r(1)   = r(1) - 2.d0 * eps_der
+    im     = j1b_nucl_square(r)
+    r(1)   = r(1) + eps_der
+    i_num  = tmp_der * (ip - im)
+    acc_ij = dabs(i_exc - i_num)
+    if(acc_ij .gt. eps_ij) then
+      print *, ' problem in grad_x list_b3 on', ipoint
+      print *, ' r      = ', r
+      print *, ' r2     = ', r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
+      print *, ' analyt = ', i_exc
+      print *, ' numeri = ', i_num
+      print *, ' diff   = ', acc_ij
+    endif
+    acc_tot += acc_ij
+    normalz += dabs(i_num)
+
+    i_exc  = v_1b_square_grad(ipoint,2)
+    r(2)   = r(2) + eps_der
+    ip     = j1b_nucl_square(r)
+    r(2)   = r(2) - 2.d0 * eps_der
+    im     = j1b_nucl_square(r)
+    r(2)   = r(2) + eps_der
+    i_num  = tmp_der * (ip - im)
+    acc_ij = dabs(i_exc - i_num)
+    if(acc_ij .gt. eps_ij) then
+      print *, ' problem in grad_y list_b3 on', ipoint
+      print *, ' r      = ', r
+      print *, ' r2     = ', r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
+      print *, ' analyt = ', i_exc
+      print *, ' numeri = ', i_num
+      print *, ' diff   = ', acc_ij
+    endif
+    acc_tot += acc_ij
+    normalz += dabs(i_num)
+
+    i_exc  = v_1b_square_grad(ipoint,3)
+    r(3)   = r(3) + eps_der
+    ip     = j1b_nucl_square(r)
+    r(3)   = r(3) - 2.d0 * eps_der
+    im     = j1b_nucl_square(r)
+    r(3)   = r(3) + eps_der
+    i_num  = tmp_der * (ip - im)
+    acc_ij = dabs(i_exc - i_num)
+    if(acc_ij .gt. eps_ij) then
+      print *, ' problem in grad_z list_b3 on', ipoint
+      print *, ' r      = ', r
+      print *, ' r2     = ', r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
+      print *, ' analyt = ', i_exc
+      print *, ' numeri = ', i_num
+      print *, ' diff   = ', acc_ij
+    endif
+    acc_tot += acc_ij
+    normalz += dabs(i_num)
+  enddo
+
+  print*, ' acc_tot on grad = ', acc_tot
+  print*, ' normalz on grad = ', normalz
+
+  ! ---
+
+  PROVIDE v_1b_square_lapl
+
+  acc_tot = 0.d0
+  normalz = 0.d0
+  do ipoint = 1, n_points_final_grid
+
+    r(1) = final_grid_points(1,ipoint)
+    r(2) = final_grid_points(2,ipoint)
+    r(3) = final_grid_points(3,ipoint)
+    i0   = j1b_nucl_square(r)
+
+    i_exc = v_1b_square_lapl(ipoint)
+
+    r(1)  = r(1) + eps_lap
+    ip    = j1b_nucl_square(r)
+    r(1)  = r(1) - 2.d0 * eps_lap
+    im    = j1b_nucl_square(r)
+    r(1)  = r(1) + eps_lap
+    i_num = tmp_lap * (ip - 2.d0 * i0 + im)
+
+    r(2)  = r(2) + eps_lap
+    ip    = j1b_nucl_square(r)
+    r(2)  = r(2) - 2.d0 * eps_lap
+    im    = j1b_nucl_square(r)
+    r(2)  = r(2) + eps_lap
+    i_num = i_num + tmp_lap * (ip - 2.d0 * i0 + im)
+
+    r(3)  = r(3) + eps_lap
+    ip    = j1b_nucl_square(r)
+    r(3)  = r(3) - 2.d0 * eps_lap
+    im    = j1b_nucl_square(r)
+    r(3)  = r(3) + eps_lap
+    i_num = i_num + tmp_lap * (ip - 2.d0 * i0 + im)
+
+    acc_ij = dabs(i_exc - i_num)
+    if(acc_ij .gt. eps_ij) then
+      print *, ' problem in lapl list_b3 on', ipoint
+      print *, ' r      = ', r
+      print *, ' r2     = ', r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
+      print *, ' analyt = ', i_exc
+      print *, ' numeri = ', i_num
+      print *, ' diff   = ', acc_ij
+    endif
+
+    acc_tot += acc_ij
+    normalz += dabs(i_num)
+  enddo
+
+  print*, ' acc_tot on lapl = ', acc_tot
+  print*, ' normalz on lapl = ', normalz
+
+  ! ---
 
   return
 end subroutine test_list_b3

src/non_h_ints_mu/grad_squared.irp.f

@@ -17,7 +17,7 @@ BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num, n_poi
   !
   ! if J(r1,r2) = u12 x v1 x v2
   !
-  ! gradu_squared_u_ij_mu = -0.50 x \int r2 \phi_i(2) \phi_j(2) [ v1^2  v2^2 ((grad_1 u12)^2 + (grad_2 u12^2)]) + u12^2  v2^2 (grad_1 v1)^2 + 2 u12 v1 v2^2 (grad_1 u12) . (grad_1 v1) ]
+  ! gradu_squared_u_ij_mu = -0.50 x \int r2 \phi_i(2) \phi_j(2) [ v1^2  v2^2 ((grad_1 u12)^2 + (grad_2 u12^2)) + u12^2  v2^2 (grad_1 v1)^2 + 2 u12 v1 v2^2 (grad_1 u12) . (grad_1 v1) ]
   !                       = -0.25 x        v1^2      \int r2 \phi_i(2) \phi_j(2) [1 - erf(mu r12)]^2 v2^2
   !                       + -0.50 x    (grad_1 v1)^2 \int r2 \phi_i(2) \phi_j(2)      u12^2          v2^2
   !                       + -1.00 x v1 (grad_1 v1)   \int r2 \phi_i(2) \phi_j(2)   (grad_1 u12)      v2^2
@@ -358,7 +358,8 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
 
   implicit none
   integer                       :: ipoint, i, j, k, l
-  double precision              :: weight1, ao_ik_r, ao_i_r
+  double precision              :: weight1, ao_k_r, ao_i_r
+  double precision              :: der_envsq_x, der_envsq_y, der_envsq_z, lap_envsq
   double precision              :: time0, time1
   double precision, allocatable :: b_mat(:,:,:), tmp(:,:,:)
 
@@ -373,16 +374,18 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
 
   else
 
+    ! ---
+
     PROVIDE int2_grad1_u12_square_ao
 
     allocate(b_mat(n_points_final_grid,ao_num,ao_num))
 
     b_mat = 0.d0
-   !$OMP PARALLEL               &
+    !$OMP PARALLEL               &
-   !$OMP DEFAULT (NONE)         &
+    !$OMP DEFAULT (NONE)         &
-   !$OMP PRIVATE (i, k, ipoint) &
+    !$OMP PRIVATE (i, k, ipoint) &
-   !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector)
+    !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector)
-   !$OMP DO SCHEDULE (static)
+    !$OMP DO SCHEDULE (static)
     do i = 1, ao_num
       do k = 1, ao_num
         do ipoint = 1, n_points_final_grid
@@ -390,13 +393,57 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
         enddo
       enddo
     enddo
-   !$OMP END DO
+    !$OMP END DO
-   !$OMP END PARALLEL
+    !$OMP END PARALLEL
 
     tc_grad_square_ao = 0.d0
     call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0                 &
               , int2_grad1_u12_square_ao(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid &
               , 0.d0, tc_grad_square_ao, ao_num*ao_num)
+
+    ! ---
+
+    if((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) then
+
+      ! an additional term is added here directly instead of 
+      ! being added in int2_grad1_u12_square_ao for performance
+      ! note that the factor
+
+      PROVIDE int2_u2_j1b2
+
+      b_mat = 0.d0
+      !$OMP PARALLEL                                                                                     &
+      !$OMP DEFAULT (NONE)                                                                               &
+      !$OMP PRIVATE (i, k, ipoint, weight1, ao_i_r, ao_k_r)                                              &
+      !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector, &
+      !$OMP         v_1b_square_grad, v_1b_square_lapl, aos_grad_in_r_array_transp_bis)
+      !$OMP DO SCHEDULE (static)
+      do i = 1, ao_num
+        do k = 1, ao_num
+          do ipoint = 1, n_points_final_grid
+
+            weight1 = 0.25d0 * final_weight_at_r_vector(ipoint)
+
+            ao_i_r = aos_in_r_array_transp(ipoint,i)
+            ao_k_r = aos_in_r_array_transp(ipoint,k)
+
+            b_mat(ipoint,k,i) = weight1 * ( ao_k_r * ao_i_r * v_1b_square_lapl(ipoint)                                                                                   &
+                              + (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,1) + ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1)) * v_1b_square_grad(ipoint,1) &
+                              + (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,2) + ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2)) * v_1b_square_grad(ipoint,2) &
+                              + (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,3) + ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3)) * v_1b_square_grad(ipoint,3) )
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+      !$OMP END PARALLEL
+
+      call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0     &
+                , int2_u2_j1b2(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid &
+                , 1.d0, tc_grad_square_ao, ao_num*ao_num)
+    endif
+
+    ! ---
+
     deallocate(b_mat)
 
     call sum_A_At(tc_grad_square_ao(1,1,1,1), ao_num*ao_num)

src/non_h_ints_mu/j12_nucl_utils.irp.f

@@ -59,7 +59,7 @@ BEGIN_PROVIDER [ double precision, v_1b, (n_points_final_grid)]
 
   else
 
-    print*, 'j1b_type = ', j1b_pen, 'is not implemented'
+    print*, 'j1b_type = ', j1b_pen, 'is not implemented for v_1b'
     stop
 
   endif
@@ -172,7 +172,7 @@ BEGIN_PROVIDER [ double precision, v_1b_lapl, (n_points_final_grid)]
   implicit none
   integer          :: ipoint, i, j, phase
   double precision :: x, y, z, dx, dy, dz
-  double precision :: a, d, e, b
+  double precision :: a, e, b
   double precision :: fact_r
   double precision :: ax_der, ay_der, az_der, a_expo
 
@@ -283,6 +283,56 @@ BEGIN_PROVIDER [ double precision, v_1b_list_b3, (n_points_final_grid)]
 
 END_PROVIDER
 
+! ---
+
+ BEGIN_PROVIDER [double precision, v_1b_square_grad, (n_points_final_grid,3)]
+&BEGIN_PROVIDER [double precision, v_1b_square_lapl, (n_points_final_grid)  ]
+
+  implicit none
+  integer          :: ipoint, i
+  double precision :: x, y, z, dx, dy, dz, r2
+  double precision :: coef, expo, a_expo, tmp
+  double precision :: fact_x, fact_y, fact_z, fact_r
+
+  PROVIDE List_all_comb_b3_coef List_all_comb_b3_expo List_all_comb_b3_cent
+
+  do ipoint = 1, n_points_final_grid
+
+    x = final_grid_points(1,ipoint)
+    y = final_grid_points(2,ipoint)
+    z = final_grid_points(3,ipoint)
+
+    fact_x = 0.d0
+    fact_y = 0.d0
+    fact_z = 0.d0
+    fact_r = 0.d0
+    do i = 1, List_all_comb_b3_size
+
+      coef = List_all_comb_b3_coef(i)
+      expo = List_all_comb_b3_expo(i)
+
+      dx = x - List_all_comb_b3_cent(1,i)
+      dy = y - List_all_comb_b3_cent(2,i)
+      dz = z - List_all_comb_b3_cent(3,i)
+      r2 = dx * dx + dy * dy + dz * dz
+
+      a_expo = expo * r2
+      tmp    = coef * expo * dexp(-a_expo)
+
+      fact_x += tmp * dx
+      fact_y += tmp * dy
+      fact_z += tmp * dz
+      fact_r += tmp * (3.d0 - 2.d0 * a_expo)
+    enddo
+
+    v_1b_square_grad(ipoint,1) = -2.d0 * fact_x
+    v_1b_square_grad(ipoint,2) = -2.d0 * fact_y
+    v_1b_square_grad(ipoint,3) = -2.d0 * fact_z
+    v_1b_square_lapl(ipoint)   = -2.d0 * fact_r
+  enddo
+
+END_PROVIDER
+
 ! ---
 
 double precision function j12_mu_r12(r12)

src/non_h_ints_mu/jast_deriv.irp.f

@@ -164,7 +164,7 @@ double precision function j12_mu(r1, r2)
   double precision, intent(in) :: r1(3), r2(3)
   double precision             :: mu_tmp, r12
 
-  if((j1b_type .ge. 100) .and. (j1b_type .lt. 200)) then
+  if((j1b_type .ge. 0) .and. (j1b_type .lt. 200)) then
 
     r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) &
                + (r1(2) - r2(2)) * (r1(2) - r2(2)) &
@@ -175,7 +175,7 @@ double precision function j12_mu(r1, r2)
 
   else
 
-    print *, ' j1b_type = ', j1b_type, 'not implemented yet'
+    print *, ' j1b_type = ', j1b_type, 'not implemented for j12_mu'
     stop
 
   endif
@@ -196,7 +196,7 @@ subroutine grad1_j12_mu(r1, r2, grad)
 
   grad = 0.d0
 
-  if((j1b_type .ge. 100) .and. (j1b_type .lt. 200)) then
+  if((j1b_type .ge. 0) .and. (j1b_type .lt. 200)) then
 
     dx = r1(1) - r2(1)
     dy = r1(2) - r2(2)
@@ -252,7 +252,7 @@ double precision function j1b_nucl(r)
   integer                      :: i
   double precision             :: a, d, e, x, y, z
 
-  if(j1b_type .eq. 102) then
+  if((j1b_type .eq. 2) .or. (j1b_type .eq. 102)) then
 
     j1b_nucl = 1.d0
     do i = 1, nucl_num
@@ -263,7 +263,7 @@ double precision function j1b_nucl(r)
       j1b_nucl = j1b_nucl - dexp(-a*dsqrt(d))
     enddo
 
-  elseif(j1b_type .eq. 103) then
+  elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then
 
     j1b_nucl = 1.d0
     do i = 1, nucl_num
@@ -275,7 +275,7 @@ double precision function j1b_nucl(r)
       j1b_nucl = j1b_nucl * e
     enddo
 
-  elseif(j1b_type .eq. 104) then
+  elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then
 
     j1b_nucl = 1.d0
     do i = 1, nucl_num
@@ -286,7 +286,7 @@ double precision function j1b_nucl(r)
       j1b_nucl = j1b_nucl - dexp(-a*d)
     enddo
 
-  elseif(j1b_type .eq. 105) then
+  elseif((j1b_type .eq. 5) .or. (j1b_type .eq. 105)) then
 
     j1b_nucl = 1.d0
     do i = 1, nucl_num
@@ -300,7 +300,7 @@ double precision function j1b_nucl(r)
 
   else
 
-    print *, ' j1b_type = ', j1b_type, 'not implemented yet'
+    print *, ' j1b_type = ', j1b_type, 'not implemented for j1b_nucl'
     stop
 
   endif
@@ -310,6 +310,75 @@ end function j1b_nucl
 
 ! ---
 
+double precision function j1b_nucl_square(r)
+
+  implicit none
+  double precision, intent(in) :: r(3)
+  integer                      :: i
+  double precision             :: a, d, e, x, y, z
+
+  if((j1b_type .eq. 2) .or. (j1b_type .eq. 102)) then
+
+    j1b_nucl_square = 1.d0
+    do i = 1, nucl_num
+      a = j1b_pen(i)
+      d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) &
+          + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) &
+          + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) )
+      j1b_nucl_square = j1b_nucl_square - dexp(-a*dsqrt(d))
+    enddo
+    j1b_nucl_square = j1b_nucl_square * j1b_nucl_square
+
+  elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then
+
+    j1b_nucl_square = 1.d0
+    do i = 1, nucl_num
+      a = j1b_pen(i)
+      d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) &
+          + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) &
+          + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) )
+      e = 1.d0 - dexp(-a*d)
+      j1b_nucl_square = j1b_nucl_square * e
+    enddo
+    j1b_nucl_square = j1b_nucl_square * j1b_nucl_square
+
+  elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then
+
+    j1b_nucl_square = 1.d0
+    do i = 1, nucl_num
+      a = j1b_pen(i)
+      d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) &
+          + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) &
+          + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) )
+      j1b_nucl_square = j1b_nucl_square - dexp(-a*d)
+    enddo
+    j1b_nucl_square = j1b_nucl_square * j1b_nucl_square
+
+  elseif((j1b_type .eq. 5) .or. (j1b_type .eq. 105)) then
+
+    j1b_nucl_square = 1.d0
+    do i = 1, nucl_num
+      a = j1b_pen(i)
+      x = r(1) - nucl_coord(i,1)
+      y = r(2) - nucl_coord(i,2)
+      z = r(3) - nucl_coord(i,3)
+      d = x*x + y*y + z*z
+      j1b_nucl_square = j1b_nucl_square - dexp(-a*d*d)
+    enddo
+    j1b_nucl_square = j1b_nucl_square * j1b_nucl_square
+
+  else
+
+    print *, ' j1b_type = ', j1b_type, 'not implemented for j1b_nucl_square'
+    stop
+
+  endif
+
+  return
+end function j1b_nucl_square
+
+! ---
+
 subroutine grad1_j1b_nucl(r, grad)
 
   implicit none
@@ -321,7 +390,7 @@ subroutine grad1_j1b_nucl(r, grad)
   double precision              :: fact_x, fact_y, fact_z
   double precision              :: ax_der, ay_der, az_der, a_expo
 
-  if(j1b_type .eq. 102) then
+  if((j1b_type .eq. 2) .or. (j1b_type .eq. 102)) then
 
     fact_x = 0.d0
     fact_y = 0.d0
@@ -343,7 +412,7 @@ subroutine grad1_j1b_nucl(r, grad)
     grad(2) = fact_y
     grad(3) = fact_z
 
-  elseif(j1b_type .eq. 103) then
+  elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then
 
     x = r(1)
     y = r(2)
@@ -382,7 +451,7 @@ subroutine grad1_j1b_nucl(r, grad)
     grad(2) = fact_y
     grad(3) = fact_z
 
-  elseif(j1b_type .eq. 104) then
+  elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then
 
     fact_x = 0.d0
     fact_y = 0.d0
@@ -404,7 +473,7 @@ subroutine grad1_j1b_nucl(r, grad)
     grad(2) = 2.d0 * fact_y
     grad(3) = 2.d0 * fact_z
 
-  elseif(j1b_type .eq. 105) then
+  elseif((j1b_type .eq. 5) .or. (j1b_type .eq. 105)) then
 
     fact_x = 0.d0
     fact_y = 0.d0
@@ -428,7 +497,7 @@ subroutine grad1_j1b_nucl(r, grad)
 
   else
 
-    print *, ' j1b_type = ', j1b_type, 'not implemented yet'
+    print *, ' j1b_type = ', j1b_type, 'not implemented for grad1_j1b_nucl'
     stop
 
   endif
@@ -520,3 +589,98 @@ subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der)
 end subroutine mu_r_val_and_grad
 
 ! ---
+
+subroutine grad1_j1b_nucl_square_num(r1, grad)
+
+  implicit none
+  double precision, intent(in)  :: r1(3)
+  double precision, intent(out) :: grad(3)
+  double precision              :: r(3), eps, tmp_eps, vp, vm
+  double precision, external    :: j1b_nucl_square
+
+  eps     = 1d-5
+  tmp_eps = 0.5d0 / eps
+
+  r(1:3) = r1(1:3)
+
+  r(1) = r(1) + eps
+  vp   = j1b_nucl_square(r)
+  r(1) = r(1) - 2.d0 * eps
+  vm   = j1b_nucl_square(r)
+  r(1) = r(1) + eps
+  grad(1) = tmp_eps * (vp - vm)
+
+  r(2) = r(2) + eps
+  vp   = j1b_nucl_square(r)
+  r(2) = r(2) - 2.d0 * eps
+  vm   = j1b_nucl_square(r)
+  r(2) = r(2) + eps
+  grad(2) = tmp_eps * (vp - vm)
+
+  r(3) = r(3) + eps
+  vp   = j1b_nucl_square(r)
+  r(3) = r(3) - 2.d0 * eps
+  vm   = j1b_nucl_square(r)
+  r(3) = r(3) + eps
+  grad(3) = tmp_eps * (vp - vm)
+  
+  return
+end subroutine grad1_j1b_nucl_square_num
+
+! ---
+
+subroutine grad1_j12_mu_square_num(r1, r2, grad)
+
+  include 'constants.include.F'
+
+  implicit none
+  double precision, intent(in)  :: r1(3), r2(3)
+  double precision, intent(out) :: grad(3)
+  double precision              :: r(3)
+  double precision              :: eps, tmp_eps, vp, vm
+  double precision, external    :: j12_mu_square
+
+  eps     = 1d-5
+  tmp_eps = 0.5d0 / eps
+
+  r(1:3) = r1(1:3)
+
+  r(1)    = r(1) + eps
+  vp      = j12_mu_square(r, r2)
+  r(1)    = r(1) - 2.d0 * eps
+  vm      = j12_mu_square(r, r2)
+  r(1)    = r(1) + eps
+  grad(1) = tmp_eps * (vp - vm)
+
+  r(2)    = r(2) + eps
+  vp      = j12_mu_square(r, r2)
+  r(2)    = r(2) - 2.d0 * eps
+  vm      = j12_mu_square(r, r2)
+  r(2)    = r(2) + eps
+  grad(2) = tmp_eps * (vp - vm)
+
+  r(3)    = r(3) + eps
+  vp      = j12_mu_square(r, r2)
+  r(3)    = r(3) - 2.d0 * eps
+  vm      = j12_mu_square(r, r2)
+  r(3)    = r(3) + eps
+  grad(3) = tmp_eps * (vp - vm)
+
+  return
+end subroutine grad1_j12_mu_square_num
+
+! ---
+
+double precision function j12_mu_square(r1, r2)
+
+  implicit none
+  double precision, intent(in) :: r1(3), r2(3)
+  double precision, external   :: j12_mu
+
+  j12_mu_square = j12_mu(r1, r2) * j12_mu(r1, r2)
+
+  return
+end function j12_mu_square
+
+! ---
+

src/non_h_ints_mu/tc_integ.irp.f

@@ -221,18 +221,21 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p
 
   elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) then
 
-    PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12
+    ! the term u12_grad1_u12_j1b_grad1_j1b is added directly for performance
+    !PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12
+    PROVIDE u12sq_j1bsq grad12_j12
 
     int2_grad1_u12_square_ao = 0.d0
     !$OMP PARALLEL               &
     !$OMP DEFAULT (NONE)         &
     !$OMP PRIVATE (i, j, ipoint) &
-    !$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, u12sq_j1bsq, u12_grad1_u12_j1b_grad1_j1b, grad12_j12)
+    !$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, u12sq_j1bsq, grad12_j12)
     !$OMP DO SCHEDULE (static)
     do ipoint = 1, n_points_final_grid
       do j = 1, ao_num
         do i = 1, ao_num
-          int2_grad1_u12_square_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint)
+          !int2_grad1_u12_square_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint)
+          int2_grad1_u12_square_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint)
         enddo
       enddo
     enddo

src/non_h_ints_mu/test_non_h_ints.irp.f

@@ -1,15 +1,25 @@
 program test_non_h
- implicit none
+
+  implicit none
+
   my_grid_becke  = .True.
   my_n_pt_r_grid = 50
   my_n_pt_a_grid = 74
+  !my_n_pt_r_grid = 400
+  !my_n_pt_a_grid = 974
+
 !  my_n_pt_r_grid = 10 ! small grid for quick debug
 !  my_n_pt_a_grid = 26 ! small grid for quick debug
   touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
-!call routine_grad_squared
+
- call routine_fit
+  !call routine_grad_squared
+  !call routine_fit
+  
+  call test_ipp()
 end
 
+! ---
+
 subroutine routine_lapl_grad
  implicit none
  integer :: i,j,k,l
@@ -100,3 +110,445 @@ subroutine routine_fit
  enddo
 
 end
+
+
+subroutine test_ipp()
+
+  implicit none
+  integer                       :: i, j, k, l, ipoint
+  double precision              :: accu, norm, diff, old, new, eps, int_num
+  double precision              :: weight1, ao_i_r, ao_k_r
+  double precision, allocatable :: b_mat(:,:,:), I1(:,:,:,:), I2(:,:,:,:)
+
+  eps = 1d-7
+
+  allocate(b_mat(n_points_final_grid,ao_num,ao_num))
+  b_mat = 0.d0
+
+  ! ---
+
+  ! first way
+
+  allocate(I1(ao_num,ao_num,ao_num,ao_num))
+  I1 = 0.d0
+
+  PROVIDE u12_grad1_u12_j1b_grad1_j1b
+
+  !$OMP PARALLEL               &
+  !$OMP DEFAULT (NONE)         &
+  !$OMP PRIVATE (i, k, ipoint) &
+  !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector)
+  !$OMP DO SCHEDULE (static)
+  do i = 1, ao_num
+    do k = 1, ao_num
+      do ipoint = 1, n_points_final_grid
+        b_mat(ipoint,k,i) = final_weight_at_r_vector(ipoint) * aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,k)
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0                    &
+            , u12_grad1_u12_j1b_grad1_j1b(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid &
+            , 0.d0, I1, ao_num*ao_num)
+
+  ! ---
+
+  ! 2nd way
+
+  allocate(I2(ao_num,ao_num,ao_num,ao_num))
+  I2 = 0.d0
+
+  PROVIDE int2_u2_j1b2
+
+  b_mat = 0.d0
+  !$OMP PARALLEL                                                                                     &
+  !$OMP DEFAULT (NONE)                                                                               &
+  !$OMP PRIVATE (i, k, ipoint, weight1, ao_i_r, ao_k_r)                                              &
+  !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector, &
+  !$OMP         v_1b_square_grad, v_1b_square_lapl, aos_grad_in_r_array_transp_bis)
+  !$OMP DO SCHEDULE (static)
+  do i = 1, ao_num
+    do k = 1, ao_num
+      do ipoint = 1, n_points_final_grid
+
+        weight1 = 0.25d0 * final_weight_at_r_vector(ipoint)
+
+        ao_i_r = aos_in_r_array_transp(ipoint,i)
+        ao_k_r = aos_in_r_array_transp(ipoint,k)
+
+        b_mat(ipoint,k,i) = weight1 * ( ao_k_r * ao_i_r * v_1b_square_lapl(ipoint)                                                                                   &
+                          + (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,1) + ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1)) * v_1b_square_grad(ipoint,1) &
+                          + (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,2) + ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2)) * v_1b_square_grad(ipoint,2) &
+                          + (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,3) + ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3)) * v_1b_square_grad(ipoint,3) )
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0     &
+            , int2_u2_j1b2(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid &
+            , 0.d0, I2, ao_num*ao_num)
+ 
+  ! ---
+
+  deallocate(b_mat)
+
+  accu = 0.d0
+  norm = 0.d0
+  do i = 1, ao_num
+    do k = 1, ao_num
+      do l = 1, ao_num
+        do j = 1, ao_num
+
+          old = I1(j,l,k,i)
+          new = I2(j,l,k,i)
+
+          !print*, l, k, j, i
+          !print*, old, new
+          
+          diff = new - old
+          if(dabs(diff) .gt. eps) then
+            print*, ' problem on :', j, l, k, i
+            print*, ' diff      = ', diff
+            print*, ' old value = ', old
+            print*, ' new value = ', new
+            call I_grade_gradu_naive1(i, j, k, l, int_num)
+            print*, ' full num1 = ', int_num
+            call I_grade_gradu_naive2(i, j, k, l, int_num)
+            print*, ' full num2 = ', int_num
+            call I_grade_gradu_naive3(i, j, k, l, int_num)
+            print*, ' full num3 = ', int_num
+            call I_grade_gradu_naive4(i, j, k, l, int_num)
+            print*, ' full num4 = ', int_num
+            call I_grade_gradu_seminaive(i, j, k, l, int_num)
+            print*, ' semi num  = ', int_num
+          endif
+
+          accu += dabs(diff)
+          norm += dabs(old)
+        enddo
+      enddo
+    enddo
+  enddo
+
+  deallocate(I1, I2)
+
+  print*, ' accu = ', accu
+  print*, ' norm = ', norm
+
+  return
+end subroutine test_ipp
+
+! ---
+
+subroutine I_grade_gradu_naive1(i, j, k, l, int)
+
+  implicit none
+  integer,          intent(in)  :: i, j, k, l
+  double precision, intent(out) :: int
+  integer                       :: ipoint, jpoint
+  double precision              :: r1(3), r2(3)
+  double precision              :: weight1_x, weight1_y, weight1_z
+  double precision              :: weight2_x, weight2_y, weight2_z
+  double precision              :: aor_i, aor_j, aor_k, aor_l
+  double precision              :: e1_val, e2_val, e1_der(3), u12_val, u12_der(3)
+  double precision, external    :: j1b_nucl, j12_mu
+
+  int = 0.d0
+
+  do ipoint = 1, n_points_final_grid ! r1
+
+    r1(1) = final_grid_points(1,ipoint)
+    r1(2) = final_grid_points(2,ipoint)
+    r1(3) = final_grid_points(3,ipoint)
+
+    aor_i = aos_in_r_array_transp(ipoint,i)
+    aor_k = aos_in_r_array_transp(ipoint,k)
+
+    e1_val = j1b_nucl(r1)
+    call grad1_j1b_nucl(r1, e1_der)
+
+    weight1_x = aor_i * aor_k * e1_val * final_weight_at_r_vector(ipoint) * e1_der(1)
+    weight1_y = aor_i * aor_k * e1_val * final_weight_at_r_vector(ipoint) * e1_der(2)
+    weight1_z = aor_i * aor_k * e1_val * final_weight_at_r_vector(ipoint) * e1_der(3)
+
+    do jpoint = 1, n_points_extra_final_grid ! r2 
+
+      r2(1) = final_grid_points_extra(1,jpoint)
+      r2(2) = final_grid_points_extra(2,jpoint)
+      r2(3) = final_grid_points_extra(3,jpoint)
+
+      aor_j = aos_in_r_array_extra_transp(jpoint,j)
+      aor_l = aos_in_r_array_extra_transp(jpoint,l)
+
+      e2_val = j1b_nucl(r2)
+
+      u12_val = j12_mu(r1, r2)
+      call grad1_j12_mu(r1, r2, u12_der)
+
+      weight2_x = aor_j * aor_l * e2_val * e2_val * u12_val * final_weight_at_r_vector_extra(jpoint) * u12_der(1)
+      weight2_y = aor_j * aor_l * e2_val * e2_val * u12_val * final_weight_at_r_vector_extra(jpoint) * u12_der(2)
+      weight2_z = aor_j * aor_l * e2_val * e2_val * u12_val * final_weight_at_r_vector_extra(jpoint) * u12_der(3)
+
+      int = int - (weight1_x * weight2_x  + weight1_y * weight2_y + weight1_z * weight2_z)
+    enddo
+  enddo
+
+  return
+end subroutine I_grade_gradu_naive1
+
+! ---
+
+subroutine I_grade_gradu_naive2(i, j, k, l, int)
+
+  implicit none
+  integer,          intent(in)  :: i, j, k, l
+  double precision, intent(out) :: int
+  integer                       :: ipoint, jpoint
+  double precision              :: r1(3), r2(3)
+  double precision              :: weight1_x, weight1_y, weight1_z
+  double precision              :: weight2_x, weight2_y, weight2_z
+  double precision              :: aor_i, aor_j, aor_k, aor_l
+  double precision              :: e1_square_der(3), e2_val, u12_square_der(3)
+  double precision, external    :: j1b_nucl
+
+  int = 0.d0
+
+  do ipoint = 1, n_points_final_grid ! r1
+
+    r1(1) = final_grid_points(1,ipoint)
+    r1(2) = final_grid_points(2,ipoint)
+    r1(3) = final_grid_points(3,ipoint)
+
+    aor_i = aos_in_r_array_transp(ipoint,i)
+    aor_k = aos_in_r_array_transp(ipoint,k)
+
+    call grad1_j1b_nucl_square_num(r1, e1_square_der)
+
+    weight1_x = aor_i * aor_k * final_weight_at_r_vector(ipoint) * e1_square_der(1)
+    weight1_y = aor_i * aor_k * final_weight_at_r_vector(ipoint) * e1_square_der(2)
+    weight1_z = aor_i * aor_k * final_weight_at_r_vector(ipoint) * e1_square_der(3)
+
+    do jpoint = 1, n_points_extra_final_grid ! r2 
+
+      r2(1) = final_grid_points_extra(1,jpoint)
+      r2(2) = final_grid_points_extra(2,jpoint)
+      r2(3) = final_grid_points_extra(3,jpoint)
+
+      aor_j = aos_in_r_array_extra_transp(jpoint,j)
+      aor_l = aos_in_r_array_extra_transp(jpoint,l)
+
+      e2_val = j1b_nucl(r2)
+      call grad1_j12_mu_square_num(r1, r2, u12_square_der)
+
+      weight2_x = aor_j * aor_l * e2_val * e2_val * final_weight_at_r_vector_extra(jpoint) * u12_square_der(1)
+      weight2_y = aor_j * aor_l * e2_val * e2_val * final_weight_at_r_vector_extra(jpoint) * u12_square_der(2)
+      weight2_z = aor_j * aor_l * e2_val * e2_val * final_weight_at_r_vector_extra(jpoint) * u12_square_der(3)
+
+      int = int - 0.25d0 * (weight1_x * weight2_x  + weight1_y * weight2_y + weight1_z * weight2_z)
+    enddo
+  enddo
+
+  return
+end subroutine I_grade_gradu_naive2
+
+! ---
+
+subroutine I_grade_gradu_naive3(i, j, k, l, int)
+
+  implicit none
+  integer,          intent(in)  :: i, j, k, l
+  double precision, intent(out) :: int
+  integer                       :: ipoint, jpoint
+  double precision              :: r1(3), r2(3)
+  double precision              :: weight1, weight2
+  double precision              :: aor_j, aor_l
+  double precision              :: grad(3), e2_val, u12_val
+  double precision, external    :: j1b_nucl, j12_mu
+
+  int = 0.d0
+
+  do ipoint = 1, n_points_final_grid ! r1
+
+    r1(1) = final_grid_points(1,ipoint)
+    r1(2) = final_grid_points(2,ipoint)
+    r1(3) = final_grid_points(3,ipoint)
+
+    call grad1_aos_ik_grad1_esquare(i, k, r1, grad)
+
+    weight1 = final_weight_at_r_vector(ipoint) * (grad(1) + grad(2) + grad(3))
+
+    do jpoint = 1, n_points_extra_final_grid ! r2 
+
+      r2(1) = final_grid_points_extra(1,jpoint)
+      r2(2) = final_grid_points_extra(2,jpoint)
+      r2(3) = final_grid_points_extra(3,jpoint)
+
+      aor_j = aos_in_r_array_extra_transp(jpoint,j)
+      aor_l = aos_in_r_array_extra_transp(jpoint,l)
+
+      e2_val  = j1b_nucl(r2)
+      u12_val = j12_mu(r1, r2)
+
+      weight2 = aor_j * aor_l * e2_val * e2_val * u12_val * u12_val * final_weight_at_r_vector_extra(jpoint)
+
+      int = int + 0.25d0 * weight1 * weight2
+    enddo
+  enddo
+
+  return
+end subroutine I_grade_gradu_naive3
+
+! ---
+
+subroutine I_grade_gradu_naive4(i, j, k, l, int)
+
+  implicit none
+  integer,          intent(in)  :: i, j, k, l
+  double precision, intent(out) :: int
+  integer                       :: ipoint, jpoint
+  double precision              :: r1(3), r2(3)
+  double precision              :: weight1, weight2
+  double precision              :: aor_j, aor_l, aor_k, aor_i
+  double precision              :: grad(3), e2_val, u12_val
+  double precision, external    :: j1b_nucl, j12_mu
+
+  int = 0.d0
+
+  do ipoint = 1, n_points_final_grid ! r1
+
+    r1(1) = final_grid_points(1,ipoint)
+    r1(2) = final_grid_points(2,ipoint)
+    r1(3) = final_grid_points(3,ipoint)
+
+    aor_i = aos_in_r_array_transp(ipoint,i)
+    aor_k = aos_in_r_array_transp(ipoint,k)
+
+    weight1 = final_weight_at_r_vector(ipoint) * ( aor_k * aor_i * v_1b_square_lapl(ipoint)                                                          &
+            + (aor_k * aos_grad_in_r_array_transp_bis(ipoint,i,1) + aor_i * aos_grad_in_r_array_transp_bis(ipoint,k,1)) * v_1b_square_grad(ipoint,1) &
+            + (aor_k * aos_grad_in_r_array_transp_bis(ipoint,i,2) + aor_i * aos_grad_in_r_array_transp_bis(ipoint,k,2)) * v_1b_square_grad(ipoint,2) &
+            + (aor_k * aos_grad_in_r_array_transp_bis(ipoint,i,3) + aor_i * aos_grad_in_r_array_transp_bis(ipoint,k,3)) * v_1b_square_grad(ipoint,3) )
+
+    do jpoint = 1, n_points_extra_final_grid ! r2 
+
+      r2(1) = final_grid_points_extra(1,jpoint)
+      r2(2) = final_grid_points_extra(2,jpoint)
+      r2(3) = final_grid_points_extra(3,jpoint)
+
+      aor_j = aos_in_r_array_extra_transp(jpoint,j)
+      aor_l = aos_in_r_array_extra_transp(jpoint,l)
+
+      e2_val  = j1b_nucl(r2)
+      u12_val = j12_mu(r1, r2)
+
+      weight2 = aor_j * aor_l * e2_val * e2_val * u12_val * u12_val * final_weight_at_r_vector_extra(jpoint)
+
+      int = int + 0.25d0 * weight1 * weight2
+    enddo
+  enddo
+
+  return
+end subroutine I_grade_gradu_naive4
+
+! ---
+
+subroutine I_grade_gradu_seminaive(i, j, k, l, int)
+
+  implicit none
+  integer,          intent(in)  :: i, j, k, l
+  double precision, intent(out) :: int
+  integer                       :: ipoint
+  double precision              :: r1(3)
+  double precision              :: weight1
+  double precision              :: aor_i, aor_k
+
+  int = 0.d0
+
+  do ipoint = 1, n_points_final_grid ! r1
+
+    aor_i = aos_in_r_array_transp(ipoint,i)
+    aor_k = aos_in_r_array_transp(ipoint,k)
+
+    weight1 = 0.25d0 * final_weight_at_r_vector(ipoint) * ( aor_k * aor_i * v_1b_square_lapl(ipoint)                                                 &
+            + (aor_k * aos_grad_in_r_array_transp_bis(ipoint,i,1) + aor_i * aos_grad_in_r_array_transp_bis(ipoint,k,1)) * v_1b_square_grad(ipoint,1) &
+            + (aor_k * aos_grad_in_r_array_transp_bis(ipoint,i,2) + aor_i * aos_grad_in_r_array_transp_bis(ipoint,k,2)) * v_1b_square_grad(ipoint,2) &
+            + (aor_k * aos_grad_in_r_array_transp_bis(ipoint,i,3) + aor_i * aos_grad_in_r_array_transp_bis(ipoint,k,3)) * v_1b_square_grad(ipoint,3) )
+
+    int = int + weight1 * int2_u2_j1b2(j,l,ipoint)
+  enddo
+
+  return
+end subroutine I_grade_gradu_seminaive
+
+! ---
+
+subroutine aos_ik_grad1_esquare(i, k, r1, val)
+
+  implicit none
+  integer,          intent(in)  :: i, k
+  double precision, intent(in)  :: r1(3)
+  double precision, intent(out) :: val(3)
+  double precision              :: tmp
+  double precision              :: der(3), aos_array(ao_num), aos_grad_array(3,ao_num)
+
+  call give_all_aos_and_grad_at_r(r1, aos_array, aos_grad_array)
+  call grad1_j1b_nucl_square_num(r1, der)
+
+  tmp    = aos_array(i) * aos_array(k)
+  val(1) = tmp * der(1)
+  val(2) = tmp * der(2)
+  val(3) = tmp * der(3)
+
+  return
+end subroutine phi_ik_grad1_esquare
+
+! ---
+
+subroutine grad1_aos_ik_grad1_esquare(i, k, r1, grad)
+
+  implicit none
+  integer,          intent(in)  :: i, k
+  double precision, intent(in)  :: r1(3)
+  double precision, intent(out) :: grad(3)
+  double precision              :: r(3), eps, tmp_eps, val_p(3), val_m(3)
+
+  eps     = 1d-5
+  tmp_eps = 0.5d0 / eps
+
+  r(1:3) = r1(1:3)
+
+  r(1) = r(1) + eps
+  call aos_ik_grad1_esquare(i, k, r, val_p)
+  r(1) = r(1) - 2.d0 * eps
+  call aos_ik_grad1_esquare(i, k, r, val_m)
+  r(1) = r(1) + eps
+  grad(1) = tmp_eps * (val_p(1) - val_m(1))
+
+  r(2) = r(2) + eps
+  call aos_ik_grad1_esquare(i, k, r, val_p)
+  r(2) = r(2) - 2.d0 * eps
+  call aos_ik_grad1_esquare(i, k, r, val_m)
+  r(2) = r(2) + eps
+  grad(2) = tmp_eps * (val_p(2) - val_m(2))
+
+  r(3) = r(3) + eps
+  call aos_ik_grad1_esquare(i, k, r, val_p)
+  r(3) = r(3) - 2.d0 * eps
+  call aos_ik_grad1_esquare(i, k, r, val_m)
+  r(3) = r(3) + eps
+  grad(3) = tmp_eps * (val_p(3) - val_m(3))
+
+  return
+end subroutine grad1_aos_ik_grad1_esquare
+
+! ---
+
+
+
+
+
+
+

src/tc_scf/tc_scf_energy.irp.f

@@ -1,5 +1,5 @@
 
- BEGIN_PROVIDER [ double precision, TC_HF_energy]
+ BEGIN_PROVIDER [ double precision, TC_HF_energy      ]
 &BEGIN_PROVIDER [ double precision, TC_HF_one_e_energy]
 &BEGIN_PROVIDER [ double precision, TC_HF_two_e_energy]