j1b_type 4/104 modif

src/ao_many_one_e_ints/listj1b.irp.f

@@ -62,6 +62,7 @@ END_PROVIDER
   double precision :: tmp_cent_x, tmp_cent_y, tmp_cent_z
 
   provide j1b_pen
+  provide j1b_pen_coef
 
   List_all_comb_b2_coef = 0.d0
   List_all_comb_b2_expo = 0.d0
@@ -127,8 +128,8 @@ END_PROVIDER
     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_coef(  i+1) = -1.d0 * j1b_pen_coef(i)
-      List_all_comb_b2_expo(  i+1) = j1b_pen(   i)
+      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)
@@ -225,6 +226,7 @@ END_PROVIDER
   double precision :: dx, dy, dz, r2
 
   provide j1b_pen
+  provide j1b_pen_coef
 
   List_all_comb_b3_coef = 0.d0
   List_all_comb_b3_expo = 0.d0
@@ -296,8 +298,8 @@ END_PROVIDER
 
     do i = 1, nucl_num
       ii = ii + 1
-      List_all_comb_b3_coef(  ii) = -2.d0
+      List_all_comb_b3_coef(  ii) = -2.d0 * j1b_pen_coef(i)
-      List_all_comb_b3_expo(  ii) = j1b_pen(   i)
+      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)
@@ -305,7 +307,7 @@ END_PROVIDER
 
     do i = 1, nucl_num
       ii = ii + 1
-      List_all_comb_b3_coef(  ii) = 1.d0
+      List_all_comb_b3_coef(  ii) = 1.d0 * j1b_pen_coef(i) * j1b_pen_coef(i)
       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)
@@ -337,7 +339,7 @@ END_PROVIDER
         
         ii = ii + 1
         ! x 2 to avoid doing integrals twice
-        List_all_comb_b3_coef(  ii) = 2.d0 * dexp(-tmp1*tmp2*tmp4*r2)
+        List_all_comb_b3_coef(  ii) = 2.d0 * dexp(-tmp1*tmp2*tmp4*r2) * j1b_pen_coef(i) * j1b_pen_coef(j)
         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)

src/non_h_ints_mu/j12_nucl_utils.irp.f

@@ -35,7 +35,7 @@ BEGIN_PROVIDER [ double precision, v_1b, (n_points_final_grid)]
 
   elseif(j1b_type .eq. 4) then
 
-    ! v(r) = 1 - \sum_{a} \exp(-\alpha_a (r - r_a)^2)
+    ! v(r) = 1 - \sum_{a} \beta_a \exp(-\alpha_a (r - r_a)^2)
 
     do ipoint = 1, n_points_final_grid
 
@@ -51,7 +51,7 @@ BEGIN_PROVIDER [ double precision, v_1b, (n_points_final_grid)]
         dz = z - nucl_coord(j,3)
         d  = dx*dx + dy*dy + dz*dz
 
-        fact_r = fact_r - dexp(-a*d)
+        fact_r = fact_r - j1b_pen_coef(j) * dexp(-a*d)
       enddo
 
       v_1b(ipoint) = fact_r
@@ -125,7 +125,7 @@ BEGIN_PROVIDER [double precision, v_1b_grad, (3, n_points_final_grid)]
 
   elseif(j1b_type .eq. 4) then
 
-    ! v(r) = 1 - \sum_{a} \exp(-\alpha_a (r - r_a)^2)
+    ! v(r) = 1 - \sum_{a} \beta_a \exp(-\alpha_a (r - r_a)^2)
 
     do ipoint = 1, n_points_final_grid
 
@@ -144,7 +144,7 @@ BEGIN_PROVIDER [double precision, v_1b_grad, (3, n_points_final_grid)]
         r2 = dx*dx + dy*dy + dz*dz
 
         a = j1b_pen(j)
-        e = a * dexp(-a * r2)
+        e = a * j1b_pen_coef(j) * dexp(-a * r2)
 
         ax_der += e * dx
         ay_der += e * dy

src/non_h_ints_mu/jast_deriv.irp.f

@@ -296,7 +296,7 @@ double precision function j1b_nucl(r)
       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 = j1b_nucl - dexp(-a*d)
+      j1b_nucl = j1b_nucl - j1b_pen_coef(i) * dexp(-a*d)
     enddo
 
   elseif((j1b_type .eq. 5) .or. (j1b_type .eq. 105)) then
@@ -363,7 +363,7 @@ double precision function j1b_nucl_square(r)
       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)
+      j1b_nucl_square = j1b_nucl_square - j1b_pen_coef(i) * dexp(-a*d)
     enddo
     j1b_nucl_square = j1b_nucl_square * j1b_nucl_square
 
@@ -475,7 +475,7 @@ subroutine grad1_j1b_nucl(r, grad)
       y = r(2) - nucl_coord(i,2)
       z = r(3) - nucl_coord(i,3)
       d = x*x + y*y + z*z
-      e = a * dexp(-a*d)
+      e = a * j1b_pen_coef(i) * dexp(-a*d)
 
       fact_x += e * x
       fact_y += e * y

src/tc_keywords/EZFIO.cfg

@@ -130,6 +130,12 @@ doc: exponents of the 1-body Jastrow
 interface: ezfio
 size: (nuclei.nucl_num)
 
+[j1b_pen_coef]
+type: double precision
+doc: coefficients of the 1-body Jastrow
+interface: ezfio
+size: (nuclei.nucl_num)
+
 [j1b_coeff]
 type: double precision
 doc: coeff of the 1-body Jastrow

src/tc_keywords/j1b_pen.irp.f

@@ -1,17 +1,22 @@
 
 ! ---
 
-BEGIN_PROVIDER [ double precision, j1b_pen, (nucl_num) ]
+ BEGIN_PROVIDER [ double precision, j1b_pen     , (nucl_num) ]
+&BEGIN_PROVIDER [ double precision, j1b_pen_coef, (nucl_num) ]
 
   BEGIN_DOC
-  ! exponents of the 1-body Jastrow
+  ! parameters of the 1-body Jastrow
   END_DOC
 
   implicit none
   logical :: exists
+  integer :: i
+  integer :: ierr
 
   PROVIDE ezfio_filename
 
+  ! ---
+
   if (mpi_master) then
     call ezfio_has_tc_keywords_j1b_pen(exists)
   endif
@@ -23,7 +28,6 @@ BEGIN_PROVIDER [ double precision, j1b_pen, (nucl_num) ]
 
   IRP_IF MPI
     include 'mpif.h'
-    integer :: ierr
     call MPI_BCAST(j1b_pen, (nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
     if (ierr /= MPI_SUCCESS) then
       stop 'Unable to read j1b_pen with MPI'
@@ -31,7 +35,6 @@ BEGIN_PROVIDER [ double precision, j1b_pen, (nucl_num) ]
   IRP_ENDIF
 
   if (exists) then
-
     if (mpi_master) then
       write(6,'(A)') '.. >>>>> [ IO READ: j1b_pen ] <<<<< ..'
       call ezfio_get_tc_keywords_j1b_pen(j1b_pen)
@@ -42,19 +45,55 @@ BEGIN_PROVIDER [ double precision, j1b_pen, (nucl_num) ]
         endif
       IRP_ENDIF
     endif
-
   else
- 
-    integer :: i
     do i = 1, nucl_num
       j1b_pen(i) = 1d5
     enddo
-
   endif
- print*,'parameters for nuclei jastrow'
+
- do i = 1, nucl_num
+  ! ---
-  print*,'i,Z,j1b_pen(i)',i,nucl_charge(i),j1b_pen(i)
+
- enddo
+  if (mpi_master) then
+    call ezfio_has_tc_keywords_j1b_pen_coef(exists)
+  endif
+
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+
+  IRP_IF MPI
+    include 'mpif.h'
+    call MPI_BCAST(j1b_pen_coef, (nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+    if (ierr /= MPI_SUCCESS) then
+      stop 'Unable to read j1b_pen_coef with MPI'
+    endif
+  IRP_ENDIF
+
+  if (exists) then
+    if (mpi_master) then
+      write(6,'(A)') '.. >>>>> [ IO READ: j1b_pen_coef ] <<<<< ..'
+      call ezfio_get_tc_keywords_j1b_pen_coef(j1b_pen_coef)
+      IRP_IF MPI
+        call MPI_BCAST(j1b_pen_coef, (nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+        if (ierr /= MPI_SUCCESS) then
+          stop 'Unable to read j1b_pen_coef with MPI'
+        endif
+      IRP_ENDIF
+    endif
+  else
+    do i = 1, nucl_num
+      j1b_pen_coef(i) = 1d0
+    enddo
+  endif
+
+  ! ---
+
+  print *, ' parameters for nuclei jastrow'
+  print *, ' i, Z, j1b_pen, j1b_pen_coef'
+  do i = 1, nucl_num
+    print *, i, nucl_charge(i), j1b_pen(i), j1b_pen_coef(i)
+  enddo
 
 END_PROVIDER
 
@@ -114,3 +153,4 @@ BEGIN_PROVIDER [ double precision, j1b_coeff, (nucl_num) ]
 END_PROVIDER
 
 ! ---
+