added jast_type 0

src/ao_many_one_e_ints/grad2_jmu_modif.irp.f

@@ -1,4 +1,72 @@
 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2, (ao_num, ao_num, n_points_final_grid)]
+
+  BEGIN_DOC
+  !
+  ! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) [1 - erf(mu r12)]^2
+  !
+  END_DOC
+
+  implicit none
+  integer                       :: i, j, ipoint, i_fit
+  double precision              :: r(3), expo_fit, coef_fit
+  double precision              :: tmp
+  double precision              :: wall0, wall1
+
+  double precision, external    :: overlap_gauss_r12_ao
+
+  print*, ' providing int2_grad1u2_grad2u2 ...'
+  call wall_time(wall0)
+
+  provide mu_erf final_grid_points j1b_pen
+
+  int2_grad1u2_grad2u2 = 0.d0
+
+ !$OMP PARALLEL DEFAULT (NONE)                                               &
+ !$OMP PRIVATE (ipoint, i, j, i_fit, r, coef_fit, expo_fit, tmp)             & 
+ !$OMP SHARED  (n_points_final_grid, ao_num, final_grid_points, ng_fit_jast, &
+ !$OMP          expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2,int2_grad1u2_grad2u2)
+ !$OMP DO
+  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)
+
+    do i = 1, ao_num
+      do j = i, ao_num
+
+        tmp = 0.d0
+        do i_fit = 1, ng_fit_jast
+
+          expo_fit = expo_gauss_1_erf_x_2(i_fit)
+          coef_fit = coef_gauss_1_erf_x_2(i_fit)
+
+          tmp += -0.25d0 * coef_fit * overlap_gauss_r12_ao(r, expo_fit, i, j)
+        enddo
+
+        int2_grad1u2_grad2u2(j,i,ipoint) = tmp
+      enddo
+    enddo
+  enddo
+ !$OMP END DO
+ !$OMP END PARALLEL
+
+  do ipoint = 1, n_points_final_grid
+    do i = 2, ao_num
+      do j = 1, i-1
+        int2_grad1u2_grad2u2(j,i,ipoint) = int2_grad1u2_grad2u2(i,j,ipoint)
+      enddo
+    enddo
+  enddo
+
+  call wall_time(wall1)
+  print*, ' wall time for int2_grad1u2_grad2u2 =', wall1 - wall0
+
+END_PROVIDER 
+
 ! ---
 
 BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n_points_final_grid)]

src/ao_tc_eff_map/compute_ints_eff_pot.irp.f

@@ -53,13 +53,13 @@ subroutine compute_ao_tc_sym_two_e_pot_jl(j, l, n_integrals, buffer_i, buffer_va
       integral_erf = ao_two_e_integral_erf(i, k, j, l)
       integral     = integral_erf + integral_pot
 
-      if( j1b_type .eq. 1 ) then
+      !if( j1b_type .eq. 1 ) then
-        !print *, ' j1b type 1 is added'
+      !  !print *, ' j1b type 1 is added'
-        integral = integral + j1b_gauss_2e_j1(i, k, j, l)
+      !  integral = integral + j1b_gauss_2e_j1(i, k, j, l)
-      elseif( j1b_type .eq. 2 ) then
+      !elseif( j1b_type .eq. 2 ) then
-        !print *, ' j1b type 2 is added'
+      !  !print *, ' j1b type 2 is added'
-        integral = integral + j1b_gauss_2e_j2(i, k, j, l)
+      !  integral = integral + j1b_gauss_2e_j2(i, k, j, l)
-      endif
+      !endif
 
       if(abs(integral) < thr) then
         cycle

src/bi_ort_ints/one_e_bi_ort.irp.f

@@ -8,22 +8,22 @@ BEGIN_PROVIDER [double precision, ao_one_e_integrals_tc_tot, (ao_num,ao_num)]
 
   ao_one_e_integrals_tc_tot = ao_one_e_integrals      
 
-  provide j1b_type
+  !provide j1b_type
 
-  if( (j1b_type .eq. 1) .or. (j1b_type .eq. 2) ) then
+  !if( (j1b_type .eq. 1) .or. (j1b_type .eq. 2) ) then
+  !
+  !  print *, ' do things properly !'
+  !  stop
 
-    print *, ' do things properly !'
+  !  !do i = 1, ao_num
-    stop
+  !  !  do j = 1, ao_num
+  !  !    ao_one_e_integrals_tc_tot(j,i) += ( j1b_gauss_hermI  (j,i) &
+  !  !                                      + j1b_gauss_hermII (j,i) &
+  !  !                                      + j1b_gauss_nonherm(j,i) )
+  !  !  enddo
+  !  !enddo
 
-    !do i = 1, ao_num
+  !endif
-    !  do j = 1, ao_num
-    !    ao_one_e_integrals_tc_tot(j,i) += ( j1b_gauss_hermI  (j,i) &
-    !                                      + j1b_gauss_hermII (j,i) &
-    !                                      + j1b_gauss_nonherm(j,i) )
-    !  enddo
-    !enddo
-
-  endif
 
 END_PROVIDER
 

src/non_h_ints_mu/jast_deriv.irp.f

@@ -19,8 +19,12 @@
   END_DOC
 
   implicit none
-  integer          :: ipoint, jpoint
+  integer                    :: ipoint, jpoint
-  double precision :: r1(3), r2(3)
+  double precision           :: r1(3), r2(3)
+  double precision           :: v1b_r1, v1b_r2, u2b_r12
+  double precision           :: grad1_v1b(3), grad1_u2b(3)
+  double precision           :: dx, dy, dz
+  double precision, external :: j12_mu, j1b_nucl
 
   PROVIDE j1b_type
   PROVIDE final_grid_points_extra
@@ -28,12 +32,43 @@
   grad1_u12_num         = 0.d0
   grad1_u12_squared_num = 0.d0
 
-  if((j1b_type .ge. 100) .and. (j1b_type .lt. 200)) then
+  if(j1b_type .eq. 100) then
 
-    double precision           :: v1b_r1, v1b_r2, u2b_r12
+    !$OMP PARALLEL                                                                                    &
-    double precision           :: grad1_v1b(3), grad1_u2b(3)
+    !$OMP DEFAULT (NONE)                                                                              &
-    double precision           :: dx, dy, dz
+    !$OMP PRIVATE (ipoint, jpoint, r1, r2, v1b_r1, v1b_r2, u2b_r12, grad1_v1b, grad1_u2b, dx, dy, dz) &
-    double precision, external :: j12_mu, j1b_nucl
+    !$OMP SHARED (n_points_final_grid, n_points_extra_final_grid, final_grid_points,                  &
+    !$OMP         final_grid_points_extra, grad1_u12_num, grad1_u12_squared_num)
+    !$OMP DO SCHEDULE (static)
+    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)
+
+      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)
+
+        call grad1_j12_mu(r1, r2, grad1_u2b)
+
+        dx = grad1_u2b(1)
+        dy = grad1_u2b(2)
+        dz = grad1_u2b(3)
+
+        grad1_u12_num(jpoint,ipoint,1) = dx 
+        grad1_u12_num(jpoint,ipoint,2) = dy 
+        grad1_u12_num(jpoint,ipoint,3) = dz 
+
+        grad1_u12_squared_num(jpoint,ipoint) = dx*dx + dy*dy + dz*dz
+      enddo
+    enddo
+    !$OMP END DO
+    !$OMP END PARALLEL
+
+  elseif((j1b_type .gt. 100) .and. (j1b_type .lt. 200)) then
 
     !$OMP PARALLEL                                                                                    &
     !$OMP DEFAULT (NONE)                                                                              &

src/non_h_ints_mu/tc_integ.irp.f

@@ -41,7 +41,34 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
 
   else
 
-    if(j1b_type .eq. 3) then
+    if(j1b_type .eq. 0) then
+
+      PROVIDE v_ij_erf_rk_cst_mu x_v_ij_erf_rk_cst_mu
+
+      int2_grad1_u12_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_u12_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_u12_ao(i,j,ipoint,1) = 0.5d0 * (tmp1 * x - x_v_ij_erf_rk_cst_mu(i,j,ipoint,1))
+            int2_grad1_u12_ao(i,j,ipoint,2) = 0.5d0 * (tmp1 * y - x_v_ij_erf_rk_cst_mu(i,j,ipoint,2))
+            int2_grad1_u12_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(j1b_type .eq. 3) then
 
       PROVIDE v_1b_grad v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b x_v_ij_erf_rk_cst_mu_j1b
 
@@ -172,7 +199,27 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p
 
   PROVIDE j1b_type
 
-  if(j1b_type .eq. 3) then
+  if(j1b_type .eq. 0) then
+
+    PROVIDE int2_grad1u2_grad2u2
+
+    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, int2_grad1u2_grad2u2)
+    !$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) = int2_grad1u2_grad2u2(i,j,ipoint)
+        enddo
+      enddo
+    enddo
+    !$OMP END DO
+    !$OMP END PARALLEL
+
+  elseif(j1b_type .eq. 3) then
 
     PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12
 

src/non_h_ints_mu/total_tc_int.irp.f

@@ -11,6 +11,13 @@ BEGIN_PROVIDER [double precision, ao_vartc_int_chemist, (ao_num, ao_num, ao_num,
   call wall_time(wall0)
   
   if(test_cycle_tc) then
+
+    PROVIDE j1b_type
+    if(j1b_type .ne. 3) then
+      print*, ' TC integrals with cycle can not be used for j1b_type =', j1b_type
+      stop
+    endif
+
     do j = 1, ao_num
       do l = 1, ao_num
         do i = 1, ao_num
@@ -20,7 +27,9 @@ BEGIN_PROVIDER [double precision, ao_vartc_int_chemist, (ao_num, ao_num, ao_num,
         enddo
       enddo
     enddo
+
   else
+
     do j = 1, ao_num
       do l = 1, ao_num
         do i = 1, ao_num
@@ -30,6 +39,7 @@ BEGIN_PROVIDER [double precision, ao_vartc_int_chemist, (ao_num, ao_num, ao_num,
         enddo
       enddo
     enddo
+
   endif
 
   call wall_time(wall1)
@@ -50,6 +60,12 @@ BEGIN_PROVIDER [double precision, ao_tc_int_chemist, (ao_num, ao_num, ao_num, ao
   
   if(test_cycle_tc) then
 
+    PROVIDE j1b_type
+    if(j1b_type .ne. 3) then
+      print*, ' TC integrals with cycle can not be used for j1b_type =', j1b_type
+      stop
+    endif
+
     ao_tc_int_chemist = ao_tc_int_chemist_test
 
   else

src/tc_keywords/EZFIO.cfg

@@ -176,7 +176,7 @@ default: 1.e-7
 type: logical
 doc: If |true|, the integrals of the three-body jastrow are computed with cycles
 interface: ezfio,provider,ocaml
-default: Flase
+default: False
 
 [thresh_biorthog_diag]
 type: Threshold