added charge-harmonizer one-body Jastrow

plugins/local/non_h_ints_mu/jast_1e.irp.f

@@ -7,6 +7,12 @@ BEGIN_PROVIDER [double precision, j1e_val, (n_points_final_grid)]
   integer          :: ipoint, i, j, p
   double precision :: x, y, z, dx, dy, dz, d2
   double precision :: a, c, tmp
+  double precision :: time0, time1
+
+  PROVIDE j1e_type
+
+  call wall_time(time0)
+  print*, ' providing j1e_val ...'
 
   if(j1e_type .eq. "none") then
 
@@ -46,29 +52,40 @@ BEGIN_PROVIDER [double precision, j1e_val, (n_points_final_grid)]
 
   else
 
-    print *, ' Error: Unknown j1e_type = ', j1e_type
+    print *, ' Error in j1e_val: Unknown j1e_type = ', j1e_type
     stop
 
   endif
 
+  call wall_time(time1)
+  print*, ' Wall time for j1e_val (min) = ', (time1 - time0) / 60.d0
+  call print_memory_usage()
+
 END_PROVIDER
 
 ! ---
 
- BEGIN_PROVIDER [double precision, j1e_dx, (n_points_final_grid)]
+ BEGIN_PROVIDER [double precision, j1e_gradx, (n_points_final_grid)]
-&BEGIN_PROVIDER [double precision, j1e_dy, (n_points_final_grid)]
+&BEGIN_PROVIDER [double precision, j1e_grady, (n_points_final_grid)]
-&BEGIN_PROVIDER [double precision, j1e_dz, (n_points_final_grid)]
+&BEGIN_PROVIDER [double precision, j1e_gradz, (n_points_final_grid)]
 
   implicit none
   integer                       :: ipoint, i, j, p
   double precision              :: x, y, z, dx, dy, dz, d2
   double precision              :: a, c, g, tmp_x, tmp_y, tmp_z
+  double precision              :: time0, time1
+  double precision, allocatable :: Pa(:,:), Pb(:,:), Pt(:,:)
+
+  PROVIDE j1e_type
+
+  call wall_time(time0)
+  print*, ' providing j1e_grad ...'
 
   if(j1e_type .eq. "none") then
 
-    j1e_dx = 0.d0
+    j1e_gradx = 0.d0
-    j1e_dy = 0.d0
+    j1e_grady = 0.d0
-    j1e_dz = 0.d0
+    j1e_gradz = 0.d0
 
   elseif(j1e_type .eq. "gauss") then
 
@@ -104,14 +121,105 @@ END_PROVIDER
         enddo
       enddo
 
-      j1e_dx(ipoint) = tmp_x
+      j1e_gradx(ipoint) = 2.d0 * tmp_x
-      j1e_dy(ipoint) = tmp_y
+      j1e_grady(ipoint) = 2.d0 * tmp_y
-      j1e_dz(ipoint) = tmp_z
+      j1e_gradz(ipoint) = 2.d0 * tmp_z
+    enddo
+
+  elseif(j1e_type .eq. "charge-harmonizer") then
+
+    ! The - sign is in the integral over r2
+    ! [(N-1)/2N] x \sum_{\mu,\nu} P_{\mu,\nu} \int dr2 [-1 * \grad_r1 J(r1,r2)] \phi_\mu(r2) \phi_nu(r2) 
+
+    PROVIDE elec_alpha_num elec_beta_num elec_num
+    PROVIDE mo_coef
+    PROVIDE int2_grad1_u2b_ao
+
+    allocate(Pa(ao_num,ao_num), Pb(ao_num,ao_num), Pt(ao_num,ao_num))
+
+    call dgemm( 'N', 'T', ao_num, ao_num, elec_alpha_num, 1.d0       &
+              , mo_coef, size(mo_coef, 1), mo_coef, size(mo_coef, 1) &
+              , 0.d0, Pa, size(Pa, 1))
+
+    if(elec_alpha_num .eq. elec_beta_num) then
+      Pb = Pa
+    else
+      call dgemm( 'N', 'T', ao_num, ao_num, elec_beta_num, 1.d0        &
+                , mo_coef, size(mo_coef, 1), mo_coef, size(mo_coef, 1) &
+                , 0.d0, Pb, size(Pb, 1))
+    endif
+    Pt = Pa + Pb
+
+    g = 0.5d0 * (dble(elec_num) - 1.d0) / dble(elec_num)
+
+    call dgemv("T", ao_num*ao_num, n_points_final_grid, g, int2_grad1_u2b_ao(1,1,1,1), ao_num*ao_num, Pt, 1, 0.d0, j1e_gradx, 1)
+    call dgemv("T", ao_num*ao_num, n_points_final_grid, g, int2_grad1_u2b_ao(1,1,1,2), ao_num*ao_num, Pt, 1, 0.d0, j1e_grady, 1)
+    call dgemv("T", ao_num*ao_num, n_points_final_grid, g, int2_grad1_u2b_ao(1,1,1,3), ao_num*ao_num, Pt, 1, 0.d0, j1e_gradz, 1)
+
+    deallocate(Pa, Pb, Pt)
+
+  else
+
+    print *, ' Error in j1e_grad: Unknown j1e_type = ', j1e_type
+    stop
+
+  endif
+
+  call wall_time(time1)
+  print*, ' Wall time for j1e_grad (min) = ', (time1 - time0) / 60.d0
+  call print_memory_usage()
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [double precision, j1e_lapl, (n_points_final_grid)]
+
+  implicit none
+  integer          :: ipoint, i, j, p
+  double precision :: x, y, z, dx, dy, dz, d2
+  double precision :: a, c, g, tmp
+
+  if(j1e_type .eq. "none") then
+
+    j1e_lapl = 0.d0
+
+  elseif(j1e_type .eq. "gauss") then
+
+    ! - \sum_{A} (r - R_A) \sum_p c_{p_A} \exp(-\alpha_{p_A} (r - R_A)^2)
+
+    PROVIDE j1e_size j1e_coef j1e_expo
+
+    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)
+
+      tmp = 0.d0
+      do j = 1, nucl_num
+
+        dx = x - nucl_coord(j,1)
+        dy = y - nucl_coord(j,2)
+        dz = z - nucl_coord(j,3)
+        d2 = dx*dx + dy*dy + dz*dz
+
+        do p = 1, j1e_size
+
+          c = j1e_coef(p,j)
+          a = j1e_expo(p,j)
+          g = c * a * dexp(-a*d2)
+
+          tmp = tmp + (2.d0 * a * d2 - 3.d0) * g
+        enddo
+      enddo
+
+      j1e_lapl(ipoint) = tmp
     enddo
 
   else
 
-    print *, ' Error: Unknown j1e_type = ', j1e_type
+    print *, ' Error in j1e_lapl: Unknown j1e_type = ', j1e_type
     stop
 
   endif
@@ -120,4 +228,3 @@ END_PROVIDER
 
 ! ---
 
-

plugins/local/non_h_ints_mu/jast_1e_utils.irp.f

@@ -0,0 +1,181 @@
+
+! ---
+
+BEGIN_PROVIDER [double precision, int2_grad1_u2b_ao, (ao_num, ao_num, n_points_final_grid, 3)]
+
+  BEGIN_DOC
+  !
+  ! int2_grad1_u2b_ao(i,j,ipoint,:) = \int dr2 [-1 * \grad_r1 J_2b(r1,r2)] \phi_i(r2) \phi_j(r2) 
+  !
+  ! where r1 = r(ipoint)
+  !
+  END_DOC
+
+  implicit none
+  integer          :: ipoint, i, j, m, jpoint
+  double precision :: time0, time1
+  double precision :: x, y, z, r2 
+  double precision :: dx, dy, dz
+  double precision :: tmp_ct
+  double precision :: tmp0, tmp1, tmp2
+  double precision :: tmp0_x, tmp0_y, tmp0_z
+  double precision :: tmp1_x, tmp1_y, tmp1_z
+
+  PROVIDE j2e_type
+
+  call wall_time(time0)
+
+  print*, ' providing int2_grad1_u2b_ao ...'
+
+  if(tc_integ_type .eq. "numeric") then
+  
+    ! TODO combine 1shot & int2_grad1_u12_ao_num
+
+    PROVIDE int2_grad1_u12_ao_num
+    int2_grad1_u2b_ao = int2_grad1_u12_ao_num
+
+    !PROVIDE int2_grad1_u12_ao_num_1shot
+    !int2_grad1_u2b_ao = int2_grad1_u12_ao_num_1shot
+
+  elseif(tc_integ_type .eq. "semi-analytic") then
+
+    ! ---
+
+    if((j2e_type .eq. "rs-dft") .and. (env_type .eq. "none")) then
+
+      PROVIDE v_ij_erf_rk_cst_mu x_v_ij_erf_rk_cst_mu
+
+      int2_grad1_u2b_ao = 0.d0
+      !$OMP PARALLEL                                                &
+      !$OMP DEFAULT (NONE)                                          &
+      !$OMP PRIVATE (ipoint, i, j, x, y, z, tmp1)                   &
+      !$OMP SHARED ( ao_num, n_points_final_grid, final_grid_points &
+      !$OMP        , v_ij_erf_rk_cst_mu, x_v_ij_erf_rk_cst_mu, int2_grad1_u2b_ao)
+      !$OMP DO SCHEDULE (static)
+      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)
+        do j = 1, ao_num
+          do i = 1, ao_num
+            tmp1 = v_ij_erf_rk_cst_mu(i,j,ipoint)
+            int2_grad1_u2b_ao(i,j,ipoint,1) = 0.5d0 * (tmp1 * x - x_v_ij_erf_rk_cst_mu(i,j,ipoint,1))
+            int2_grad1_u2b_ao(i,j,ipoint,2) = 0.5d0 * (tmp1 * y - x_v_ij_erf_rk_cst_mu(i,j,ipoint,2))
+            int2_grad1_u2b_ao(i,j,ipoint,3) = 0.5d0 * (tmp1 * z - x_v_ij_erf_rk_cst_mu(i,j,ipoint,3))
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+      !$OMP END PARALLEL
+
+    elseif((j2e_type .eq. "rs-dft") .and. (env_type .eq. "prod-gauss")) then
+
+      PROVIDE env_type env_val env_grad
+      PROVIDE v_ij_erf_rk_cst_mu_env v_ij_u_cst_mu_env_an x_v_ij_erf_rk_cst_mu_env
+
+      int2_grad1_u2b_ao = 0.d0
+      !$OMP PARALLEL                                                                   &
+      !$OMP DEFAULT (NONE)                                                             &
+      !$OMP PRIVATE (ipoint, i, j, x, y, z, tmp0, tmp1, tmp2, tmp0_x, tmp0_y, tmp0_z)  &
+      !$OMP SHARED (ao_num, n_points_final_grid, final_grid_points, env_val, env_grad, &
+      !$OMP        v_ij_erf_rk_cst_mu_env, v_ij_u_cst_mu_env_an, x_v_ij_erf_rk_cst_mu_env, int2_grad1_u2b_ao)
+      !$OMP DO SCHEDULE (static)
+      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)
+        tmp0   =     0.5d0 * env_val(ipoint)
+        tmp0_x =          env_grad(1,ipoint)
+        tmp0_y =          env_grad(2,ipoint)
+        tmp0_z =          env_grad(3,ipoint)
+        do j = 1, ao_num
+          do i = 1, ao_num
+            tmp1 = tmp0 * v_ij_erf_rk_cst_mu_env(i,j,ipoint)
+            tmp2 = v_ij_u_cst_mu_env_an(i,j,ipoint)
+            int2_grad1_u2b_ao(i,j,ipoint,1) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_env(i,j,ipoint,1) - tmp2 * tmp0_x
+            int2_grad1_u2b_ao(i,j,ipoint,2) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_env(i,j,ipoint,2) - tmp2 * tmp0_y
+            int2_grad1_u2b_ao(i,j,ipoint,3) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_env(i,j,ipoint,3) - tmp2 * tmp0_z
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+      !$OMP END PARALLEL
+
+    elseif((j2e_type .eq. "rs-dft") .and. (env_type .eq. "sum-gauss")) then
+
+      PROVIDE mu_erf
+      PROVIDE env_type env_val env_grad
+      PROVIDE Ir2_rsdft_long_Du_0 Ir2_rsdft_long_Du_x Ir2_rsdft_long_Du_y Ir2_rsdft_long_Du_z Ir2_rsdft_long_Du_2
+      PROVIDE Ir2_rsdft_gauss_Du
+
+      tmp_ct = 0.5d0 / (dsqrt(dacos(-1.d0)) * mu_erf)
+
+      int2_grad1_u2b_ao = 0.d0
+
+      !$OMP PARALLEL                                                    &
+      !$OMP DEFAULT (NONE)                                              &
+      !$OMP PRIVATE (ipoint, i, j, x, y, z, r2, dx, dy, dz, tmp1, tmp2, & 
+      !$OMP         tmp0_x, tmp0_y, tmp0_z, tmp1_x, tmp1_y, tmp1_z)     &
+      !$OMP SHARED (ao_num, n_points_final_grid, final_grid_points,     &
+      !$OMP         tmp_ct, env_val, env_grad, Ir2_rsdft_long_Du_0,     &
+      !$OMP         Ir2_rsdft_long_Du_x, Ir2_rsdft_long_Du_y,           &
+      !$OMP         Ir2_rsdft_long_Du_z, Ir2_rsdft_gauss_Du,            &
+      !$OMP         Ir2_rsdft_long_Du_2, int2_grad1_u2b_ao)
+      !$OMP DO SCHEDULE (static)
+      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)
+        r2 = x*x + y*y + z*z
+
+        dx = env_grad(1,ipoint)
+        dy = env_grad(2,ipoint)
+        dz = env_grad(3,ipoint)
+
+        tmp0_x = 0.5d0 * (env_val(ipoint) * x + r2 * dx)
+        tmp0_y = 0.5d0 * (env_val(ipoint) * y + r2 * dy)
+        tmp0_z = 0.5d0 * (env_val(ipoint) * z + r2 * dz)
+
+        tmp1 = 0.5d0 * env_val(ipoint)
+
+        tmp1_x = tmp_ct * dx
+        tmp1_y = tmp_ct * dy
+        tmp1_z = tmp_ct * dz
+
+        do j = 1, ao_num
+          do i = 1, ao_num
+ 
+            tmp2 = 0.5d0 * Ir2_rsdft_long_Du_2(i,j,ipoint) - x * Ir2_rsdft_long_Du_x(i,j,ipoint) - y * Ir2_rsdft_long_Du_y(i,j,ipoint) - z * Ir2_rsdft_long_Du_z(i,j,ipoint)
+
+            int2_grad1_u2b_ao(i,j,ipoint,1) = -Ir2_rsdft_long_Du_0(i,j,ipoint) * tmp0_x + tmp1 * Ir2_rsdft_long_Du_x(i,j,ipoint) - dx * tmp2 + tmp1_x * Ir2_rsdft_gauss_Du(i,j,ipoint)
+            int2_grad1_u2b_ao(i,j,ipoint,2) = -Ir2_rsdft_long_Du_0(i,j,ipoint) * tmp0_y + tmp1 * Ir2_rsdft_long_Du_y(i,j,ipoint) - dy * tmp2 + tmp1_y * Ir2_rsdft_gauss_Du(i,j,ipoint)
+            int2_grad1_u2b_ao(i,j,ipoint,3) = -Ir2_rsdft_long_Du_0(i,j,ipoint) * tmp0_z + tmp1 * Ir2_rsdft_long_Du_z(i,j,ipoint) - dz * tmp2 + tmp1_z * Ir2_rsdft_gauss_Du(i,j,ipoint)
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+      !$OMP END PARALLEL
+
+    else 
+
+      print *, ' Error in int2_grad1_u2b_ao: Unknown Jastrow'
+      stop
+
+    endif ! j2e_type
+
+  else
+
+    write(*, '(A, A, A)') ' Error: The integration type ', trim(tc_integ_type), ' has not been implemented yet'
+    stop
+
+  endif ! tc_integ_type
+
+  call wall_time(time1)
+  print*, ' wall time for int2_grad1_u2b_ao (min) =', (time1-time0)/60.d0
+  call print_memory_usage()
+
+END_PROVIDER
+
+! ---
+

plugins/local/non_h_ints_mu/tc_integ.irp.f

@@ -119,8 +119,6 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
         !$OMP END DO
         !$OMP END PARALLEL
 
-        FREE v_ij_erf_rk_cst_mu_env v_ij_u_cst_mu_env_an x_v_ij_erf_rk_cst_mu_env
-
       elseif((j2e_type .eq. "rs-dft") .and. (env_type .eq. "sum-gauss")) then
 
         PROVIDE mu_erf
@@ -190,7 +188,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
 
         PROVIDE elec_num
         PROVIDE ao_overlap
-        PROVIDE j1e_dx j1e_dy j1e_dz
+        PROVIDE j1e_gradx j1e_grady j1e_gradz
 
         tmp_ct = 1.d0 / (dble(elec_num) - 1.d0)
 
@@ -198,12 +196,12 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
         !$OMP DEFAULT (NONE)                                 &
         !$OMP PRIVATE (ipoint, i, j, tmp0_x, tmp0_y, tmp0_z) &
         !$OMP SHARED (ao_num, n_points_final_grid, tmp_ct,   &
-        !$OMP         j1e_dx, j1e_dy, j1e_dz, ao_overlap, int2_grad1_u12_ao)
+        !$OMP         j1e_gradx, j1e_grady, j1e_gradz, ao_overlap, int2_grad1_u12_ao)
         !$OMP DO SCHEDULE (static)
         do ipoint = 1, n_points_final_grid
-          tmp0_x = tmp_ct * j1e_dx(ipoint)
+          tmp0_x = tmp_ct * j1e_gradx(ipoint)
-          tmp0_y = tmp_ct * j1e_dy(ipoint)
+          tmp0_y = tmp_ct * j1e_grady(ipoint)
-          tmp0_z = tmp_ct * j1e_dz(ipoint)
+          tmp0_z = tmp_ct * j1e_gradz(ipoint)
           do j = 1, ao_num
             do i = 1, ao_num
               int2_grad1_u12_ao(i,j,ipoint,1) = int2_grad1_u12_ao(i,j,ipoint,1) + tmp0_x * ao_overlap(i,j)
@@ -217,7 +215,13 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
 
       else
 
+        if((j2e_type .eq. "rs-dft") .and. (env_type .eq. "none")) then
+          FREE v_ij_erf_rk_cst_mu x_v_ij_erf_rk_cst_mu
+        elseif((j2e_type .eq. "rs-dft") .and. (env_type .eq. "prod-gauss")) then
+          FREE v_ij_erf_rk_cst_mu_env v_ij_u_cst_mu_env_an x_v_ij_erf_rk_cst_mu_env
+        elseif((j2e_type .eq. "rs-dft") .and. (env_type .eq. "sum-gauss")) then
           FREE Ir2_rsdft_long_Du_0 Ir2_rsdft_long_Du_x Ir2_rsdft_long_Du_y Ir2_rsdft_long_Du_z Ir2_rsdft_gauss_Du Ir2_rsdft_long_Du_2
+        endif
 
       endif ! j1e_type
 
@@ -519,7 +523,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p
 
       PROVIDE elec_num
       PROVIDE ao_overlap
-      PROVIDE j1e_dx j1e_dy j1e_dz
+      PROVIDE j1e_gradx j1e_grady j1e_gradz
 
       tmp_ct1 = 1.0d0 / (dsqrt(dacos(-1.d0)) * mu_erf)
       tmp_ct2 = 1.0d0 / (dble(elec_num) - 1.d0)
@@ -531,7 +535,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p
       !$OMP         tmp0_x, tmp0_y, tmp0_z)                            &
       !$OMP SHARED (ao_num, n_points_final_grid, final_grid_points,    &
       !$OMP         tmp_ct1, tmp_ct2, env_val, env_grad,               &
-      !$OMP         j1e_dx, j1e_dy, j1e_dz,                            &
+      !$OMP         j1e_gradx, j1e_grady, j1e_gradz,                   &
       !$OMP         Ir2_rsdft_long_Du_0, Ir2_rsdft_long_Du_2,          &
       !$OMP         Ir2_rsdft_long_Du_x, Ir2_rsdft_long_Du_y,          &
       !$OMP         Ir2_rsdft_long_Du_z, Ir2_rsdft_gauss_Du,           &
@@ -548,9 +552,9 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p
         dy1 = env_grad(2,ipoint)
         dz1 = env_grad(3,ipoint)
 
-        dx2 = j1e_dx(ipoint)
+        dx2 = j1e_gradx(ipoint)
-        dy2 = j1e_dy(ipoint)
+        dy2 = j1e_grady(ipoint)
-        dz2 = j1e_dz(ipoint)
+        dz2 = j1e_gradz(ipoint)
 
         dr12 = dx1*dx2 + dy1*dy2 + dz1*dz2
 

plugins/local/non_h_ints_mu/test_non_h_ints.irp.f

@@ -26,16 +26,20 @@ program test_non_h
   
   !call test_v_ij_u_cst_mu_env_an()
 
-  call test_int2_grad1_u12_square_ao()
+  !call test_int2_grad1_u12_square_ao()
-  call test_int2_grad1_u12_ao()
+  !call test_int2_grad1_u12_ao()
+
+  call test_j1e_grad()
 end
 
 ! ---
 
 subroutine routine_fit
+
   implicit none
   integer :: i,nx
   double precision :: dx,xmax,x,j_mu,j_mu_F_x_j,j_mu_fit_gauss
+ 
   nx = 500
   xmax = 5.d0
   dx = xmax/dble(nx)
@@ -48,6 +52,7 @@ subroutine routine_fit
 
 end
 
+! ---
 
 subroutine test_ipp()
 
@@ -561,7 +566,7 @@ subroutine test_int2_grad1_u12_square_ao()
   print*, ' accuracy(%) = ', 100.d0 * accu / norm
 
   return
-end subroutine test_int2_grad1_u12_square_ao
+end
 
 ! ---
 
@@ -605,7 +610,108 @@ subroutine test_int2_grad1_u12_ao()
   print*, ' accuracy(%) = ', 100.d0 * accu / norm
 
   return
-end subroutine test_int2_grad1_u12_ao
+end
+
+! ---
+
+subroutine test_j1e_grad()
+
+  implicit none
+  integer                       :: i, j, ipoint
+  double precision              :: g
+  double precision              :: x_loops, x_dgemm, diff, thr, accu, norm
+  double precision, allocatable :: pa(:,:), Pb(:,:), Pt(:,:)
+  double precision, allocatable :: x(:), y(:), z(:)
+  
+  PROVIDE int2_grad1_u2b_ao
+  PROVIDE mo_coef
+
+  allocate(Pa(ao_num,ao_num), Pb(ao_num,ao_num), Pt(ao_num,ao_num))
+
+  call dgemm( 'N', 'T', ao_num, ao_num, elec_alpha_num, 1.d0       &
+            , mo_coef, size(mo_coef, 1), mo_coef, size(mo_coef, 1) &
+            , 0.d0, Pa, size(Pa, 1))
+
+  if(elec_alpha_num .eq. elec_beta_num) then
+    Pb = Pa
+  else
+    call dgemm( 'N', 'T', ao_num, ao_num, elec_beta_num, 1.d0        &
+              , mo_coef, size(mo_coef, 1), mo_coef, size(mo_coef, 1) &
+              , 0.d0, Pb, size(Pb, 1))
+  endif
+  Pt = Pa + Pa
+
+
+  g = 0.5d0 * (dble(elec_num) - 1.d0) / dble(elec_num)
+
+  allocate(x(n_points_final_grid), y(n_points_final_grid), z(n_points_final_grid))
+
+  do ipoint = 1, n_points_final_grid
+    x(ipoint) = 0.d0
+    y(ipoint) = 0.d0
+    z(ipoint) = 0.d0
+    do i = 1, ao_num
+      do j = 1, ao_num
+        x(ipoint) = x(ipoint) + g * Pt(i,j) * int2_grad1_u2b_ao(i,j,ipoint,1)
+        y(ipoint) = y(ipoint) + g * Pt(i,j) * int2_grad1_u2b_ao(i,j,ipoint,2)
+        z(ipoint) = z(ipoint) + g * Pt(i,j) * int2_grad1_u2b_ao(i,j,ipoint,3)
+      enddo
+    enddo
+  enddo
+
+  deallocate(Pa, Pb, Pt)
+
+  ! ---
+
+  thr  = 1d-10
+  norm = 0.d0
+  accu = 0.d0
+  do ipoint = 1, n_points_final_grid
+
+    x_loops = x        (ipoint)
+    x_dgemm = j1e_gradx(ipoint)
+    diff    = dabs(x_loops - x_dgemm)
+    if(diff .gt. thr) then
+      print *, ' problem in j1e_gradx on:', ipoint
+      print *, ' loops :', x_loops
+      print *, ' dgemm :', x_dgemm
+      stop
+    endif
+    accu += diff
+    norm += dabs(x_loops)
+
+    x_loops = y        (ipoint)
+    x_dgemm = j1e_grady(ipoint)
+    diff    = dabs(x_loops - x_dgemm)
+    if(diff .gt. thr) then
+      print *, ' problem in j1e_grady on:', ipoint
+      print *, ' loops :', x_loops
+      print *, ' dgemm :', x_dgemm
+      stop
+    endif
+    accu += diff
+    norm += dabs(x_loops)
+
+    x_loops = z        (ipoint)
+    x_dgemm = j1e_gradz(ipoint)
+    diff    = dabs(x_loops - x_dgemm)
+    if(diff .gt. thr) then
+      print *, ' problem in j1e_gradz on:', ipoint
+      print *, ' loops :', x_loops
+      print *, ' dgemm :', x_dgemm
+      stop
+    endif
+    accu += diff
+    norm += dabs(x_loops)
+
+  enddo
+
+  deallocate(x, y, z)
+
+  print*, ' accuracy(%) = ', 100.d0 * accu / norm
+
+  return
+end
 
 ! ---