jast 4 added

src/ao_many_one_e_ints/NEED

@@ -3,3 +3,4 @@ ao_two_e_ints
 becke_numerical_grid
 mo_one_e_ints
 dft_utils_in_r
+tc_keywords

src/ao_many_one_e_ints/listj1b.irp.f

@@ -5,8 +5,25 @@ BEGIN_PROVIDER [ integer, List_all_comb_b2_size]
 
   implicit none
 
+  PROVIDE j1b_type
+
+  if(j1b_type .eq. 3) then
+
     List_all_comb_b2_size = 2**nucl_num
 
+  elseif(j1b_type .eq. 4) then
+
+    List_all_comb_b2_size = nucl_num + 1
+
+  else 
+
+    print *, 'j1b_type = ', j1b_pen, 'is not implemented'
+    stop
+
+  endif
+
+  print *, ' nb of linear terms in the envelope is ', List_all_comb_b2_size
+
 END_PROVIDER
 
 ! ---
@@ -50,6 +67,8 @@ END_PROVIDER
   List_all_comb_b2_expo = 0.d0
   List_all_comb_b2_cent = 0.d0
 
+  if(j1b_type .eq. 3) then
+
     do i = 1, List_all_comb_b2_size
 
       tmp_cent_x = 0.d0
@@ -102,6 +121,26 @@ END_PROVIDER
       List_all_comb_b2_coef(i) = (-1.d0)**dble(phase) * dexp(-List_all_comb_b2_coef(i))
     enddo
 
+  elseif(j1b_type .eq. 4) then
+
+    List_all_comb_b2_coef(    1) = 1.d0
+    List_all_comb_b2_expo(    1) = 0.d0
+    List_all_comb_b2_cent(1:3,1) = 0.d0
+    do i = 1, nucl_num
+      List_all_comb_b2_coef(  i+1) = -1.d0
+      List_all_comb_b2_expo(  i+1) = j1b_pen(   i)
+      List_all_comb_b2_cent(1,i+1) = nucl_coord(i,1)
+      List_all_comb_b2_cent(2,i+1) = nucl_coord(i,2)
+      List_all_comb_b2_cent(3,i+1) = nucl_coord(i,3)
+    enddo
+
+  else
+
+    print *, 'j1b_type = ', j1b_pen, 'is not implemented'
+    stop
+
+  endif
+
   !print *, ' coeff, expo & cent of list b2'
   !do i = 1, List_all_comb_b2_size
   !  print*, i, List_all_comb_b2_coef(i), List_all_comb_b2_expo(i)
@@ -115,9 +154,26 @@ END_PROVIDER
 BEGIN_PROVIDER [ integer, List_all_comb_b3_size]
 
   implicit none
+  double precision :: tmp
+
+  if(j1b_type .eq. 3) then
 
     List_all_comb_b3_size = 3**nucl_num
 
+  elseif(j1b_type .eq. 4) then
+
+    tmp                   = 0.5d0 * dble(nucl_num) * (dble(nucl_num) + 3.d0)
+    List_all_comb_b3_size = int(tmp) + 1
+
+  else
+
+    print *, 'j1b_type = ', j1b_pen, 'is not implemented'
+    stop
+
+  endif
+
+  print *, ' nb of linear terms in the square of the envelope is ', List_all_comb_b3_size
+
 END_PROVIDER
 
 ! ---
@@ -162,7 +218,11 @@ END_PROVIDER
 
   implicit none
   integer          :: i, j, k, phase
+  integer          :: ii
   double precision :: tmp_alphaj, tmp_alphak, facto
+  double precision :: tmp1, tmp2, tmp3, tmp4
+  double precision :: xi, yi, zi, xj, yj, zj
+  double precision :: dx, dy, dz, r2
 
   provide j1b_pen
 
@@ -170,6 +230,8 @@ END_PROVIDER
   List_all_comb_b3_expo = 0.d0
   List_all_comb_b3_cent = 0.d0
 
+  if(j1b_type .eq. 3) then
+
     do i = 1, List_all_comb_b3_size
 
       do j = 1, nucl_num
@@ -225,6 +287,70 @@ END_PROVIDER
       List_all_comb_b3_coef(i) = (-1.d0)**dble(phase) * facto * dexp(-List_all_comb_b3_coef(i))
     enddo
 
+  elseif(j1b_type .eq. 4) then
+
+    ii = 1
+    List_all_comb_b3_coef(    ii) = 1.d0
+    List_all_comb_b3_expo(    ii) = 0.d0
+    List_all_comb_b3_cent(1:3,ii) = 0.d0
+
+    do i = 1, nucl_num
+      ii = ii + 1
+      List_all_comb_b3_coef(  ii) = -2.d0
+      List_all_comb_b3_expo(  ii) = j1b_pen(   i)
+      List_all_comb_b3_cent(1,ii) = nucl_coord(i,1)
+      List_all_comb_b3_cent(2,ii) = nucl_coord(i,2)
+      List_all_comb_b3_cent(3,ii) = nucl_coord(i,3)
+    enddo
+
+    do i = 1, nucl_num
+      ii = ii + 1
+      List_all_comb_b3_coef(  ii) = 1.d0
+      List_all_comb_b3_expo(  ii) = 2.d0 * j1b_pen(i)
+      List_all_comb_b3_cent(1,ii) = nucl_coord(i,1)
+      List_all_comb_b3_cent(2,ii) = nucl_coord(i,2)
+      List_all_comb_b3_cent(3,ii) = nucl_coord(i,3)
+    enddo
+
+    do i = 1, nucl_num-1
+
+      tmp1 = j1b_pen(i)
+
+      xi = nucl_coord(i,1)
+      yi = nucl_coord(i,2)
+      zi = nucl_coord(i,3)
+
+      do j = i+1, nucl_num
+
+        tmp2 = j1b_pen(j)
+        tmp3 = tmp1 + tmp2
+        tmp4 = 1.d0 / tmp3
+
+        xj = nucl_coord(j,1)
+        yj = nucl_coord(j,2)
+        zj = nucl_coord(j,3)
+
+        dx = xi - xj
+        dy = yi - yj
+        dz = zi - zj
+        r2 = dx*dx + dy*dy + dz*dz
+        
+        ii = ii + 1
+        List_all_comb_b3_coef(  ii) = dexp(-tmp1*tmp2*tmp4*r2)
+        List_all_comb_b3_expo(  ii) = tmp3
+        List_all_comb_b3_cent(1,ii) = tmp4 * (tmp1 * xi + tmp2 * xj)
+        List_all_comb_b3_cent(2,ii) = tmp4 * (tmp1 * yi + tmp2 * yj)
+        List_all_comb_b3_cent(3,ii) = tmp4 * (tmp1 * zi + tmp2 * zj)
+      enddo
+    enddo
+
+  else
+
+    print *, 'j1b_type = ', j1b_pen, 'is not implemented'
+    stop
+
+  endif
+
   !print *, ' coeff, expo & cent of list b3'
   !do i = 1, List_all_comb_b3_size
   !  print*, i, List_all_comb_b3_coef(i), List_all_comb_b3_expo(i)

src/non_h_ints_mu/j12_nucl_utils.irp.f

@@ -8,6 +8,10 @@ BEGIN_PROVIDER [ double precision, v_1b, (n_points_final_grid)]
   double precision :: x, y, z, dx, dy, dz
   double precision :: a, d, e, fact_r
 
+  if(j1b_type .eq. 3) then
+
+    ! v(r) = \Pi_{a} [1 - \exp(-\alpha_a (r - r_a)^2)]
+
     do ipoint = 1, n_points_final_grid
 
       x = final_grid_points(1,ipoint)
@@ -29,6 +33,37 @@ BEGIN_PROVIDER [ double precision, v_1b, (n_points_final_grid)]
       v_1b(ipoint) = fact_r
     enddo
 
+  elseif(j1b_type .eq. 4) then
+
+    ! v(r) = 1 - \sum_{a} \exp(-\alpha_a (r - r_a)^2)
+
+    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_r = 1.d0
+      do j = 1, nucl_num
+        a  = j1b_pen(j)
+        dx = x - nucl_coord(j,1)
+        dy = y - nucl_coord(j,2)
+        dz = z - nucl_coord(j,3)
+        d  = dx*dx + dy*dy + dz*dz
+
+        fact_r = fact_r - dexp(-a*d)
+      enddo
+
+      v_1b(ipoint) = fact_r
+    enddo
+
+  else
+
+    print*, 'j1b_type = ', j1b_pen, 'is not implemented'
+    stop
+
+  endif
+
 END_PROVIDER
 
 ! ---
@@ -37,11 +72,17 @@ BEGIN_PROVIDER [ double precision, v_1b_grad, (3, n_points_final_grid)]
 
   implicit none
   integer          :: ipoint, i, j, phase
-  double precision :: x, y, z, dx, dy, dz
+  double precision :: x, y, z, dx, dy, dz, r2
   double precision :: a, d, e
   double precision :: fact_x, fact_y, fact_z
   double precision :: ax_der, ay_der, az_der, a_expo
 
+  PROVIDE j1b_type
+
+  if(j1b_type .eq. 3) then
+
+    ! v(r) = \Pi_{a} [1 - \exp(-\alpha_a (r - r_a)^2)]
+
     do ipoint = 1, n_points_final_grid
 
       x = final_grid_points(1,ipoint)
@@ -82,6 +123,46 @@ BEGIN_PROVIDER [ double precision, v_1b_grad, (3, n_points_final_grid)]
       v_1b_grad(3,ipoint) = fact_z
     enddo
 
+  elseif(j1b_type .eq. 4) then
+
+    ! v(r) = 1 - \sum_{a} \exp(-\alpha_a (r - r_a)^2)
+
+    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)
+
+      ax_der = 0.d0 
+      ay_der = 0.d0 
+      az_der = 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)
+        r2 = dx*dx + dy*dy + dz*dz
+
+        a = j1b_pen(j)
+        e = a * dexp(-a * r2)
+
+        ax_der += e * dx
+        ay_der += e * dy
+        az_der += e * dz
+      enddo
+
+      v_1b_grad(1,ipoint) = 2.d0 * ax_der
+      v_1b_grad(2,ipoint) = 2.d0 * ay_der
+      v_1b_grad(3,ipoint) = 2.d0 * az_der
+    enddo
+
+  else
+
+    print*, 'j1b_type = ', j1b_pen, 'is not implemented'
+    stop
+
+  endif
+
 END_PROVIDER
 
 ! ---

src/non_h_ints_mu/tc_integ.irp.f

@@ -68,7 +68,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
       !$OMP END DO
       !$OMP END PARALLEL
 
-    elseif(j1b_type .eq. 3) then
+    elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) 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
 
@@ -219,7 +219,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p
     !$OMP END DO
     !$OMP END PARALLEL
 
-  elseif(j1b_type .eq. 3) then
+  elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) then
 
     PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12