added Boys & Handy's Jastrow

plugins/local/jastrow/EZFIO.cfg

@@ -1,13 +1,13 @@
 
 [j2e_type]
 type: character*(32)
-doc: type of the 2e-Jastrow: [ None | Mu | Mur | Qmckl ]
+doc: type of the 2e-Jastrow: [ None | Mu | Mu_Nu | Mur | Boys | Boys_Handy | Qmckl ]
 interface: ezfio,provider,ocaml
 default: Mu
 
 [j1e_type]
 type: character*(32)
-doc: type of the 1e-Jastrow: [ None | Gauss | Charge_Harmonizer ]
+doc: type of the 1e-Jastrow: [ None | Gauss | Charge_Harmonizer | Charge_Harmonizer_AO ]
 interface: ezfio,provider,ocaml
 default: None
 
@@ -151,3 +151,45 @@ interface: ezfio,provider,ocaml
 default: 1.0
 ezfio_name: nu_erf
 
+[jBH_size]
+type: integer
+doc: number of terms per atom in Boys-Handy-Jastrow
+interface: ezfio,provider,ocaml
+default: 1
+
+[jBH_c]
+type: double precision
+doc: coefficients of terms in Boys-Handy-Jastrow
+interface: ezfio
+size: (jastrow.jBH_size,nuclei.nucl_num)
+
+[jBH_m]
+type: integer
+doc: powers of terms in Boys-Handy-Jastrow
+interface: ezfio
+size: (jastrow.jBH_size,nuclei.nucl_num)
+
+[jBH_n]
+type: integer
+doc: powers of terms in Boys-Handy-Jastrow
+interface: ezfio
+size: (jastrow.jBH_size,nuclei.nucl_num)
+
+[jBH_o]
+type: integer
+doc: powers of terms in Boys-Handy-Jastrow
+interface: ezfio
+size: (jastrow.jBH_size,nuclei.nucl_num)
+
+[jBH_ee]
+type: double precision
+doc: parameters of e-e terms in Boys-Handy-Jastrow
+interface: ezfio
+size: (nuclei.nucl_num)
+
+[jBH_en]
+type: double precision
+doc: parameters of e-n terms in Boys-Handy-Jastrow
+interface: ezfio
+size: (nuclei.nucl_num)
+

plugins/local/jastrow/bh_param.irp.f

@@ -0,0 +1,252 @@
+
+ BEGIN_PROVIDER [double precision, jBH_ee,           (nucl_num)]
+&BEGIN_PROVIDER [double precision, jBH_en,           (nucl_num)]
+&BEGIN_PROVIDER [double precision, jBH_c , (jBH_size, nucl_num)]
+&BEGIN_PROVIDER [integer         , jBH_m , (jBH_size, nucl_num)]
+&BEGIN_PROVIDER [integer         , jBH_n , (jBH_size, nucl_num)]
+&BEGIN_PROVIDER [integer         , jBH_o , (jBH_size, nucl_num)]
+
+  BEGIN_DOC
+  !
+  ! parameters of Boys-Handy-Jastrow
+  !
+  END_DOC
+
+  implicit none
+  logical :: exists
+  integer :: i_nucl, p
+  integer :: ierr
+
+  PROVIDE ezfio_filename
+
+  ! ---
+
+  if(mpi_master) then
+    call ezfio_has_jastrow_jBH_ee(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(jBH_ee, (nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+    if(ierr /= MPI_SUCCESS) then
+      stop 'Unable to read Boys-Handy e-e param with MPI'
+    endif
+  IRP_ENDIF
+
+  if(exists) then
+    if(mpi_master) then
+      write(6,'(A)') '.. >>>>> [ IO READ: jBH_ee ] <<<<< ..'
+      call ezfio_get_jastrow_jBH_ee(jBH_ee)
+      IRP_IF MPI
+        call MPI_BCAST(jBH_ee, (nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+        if(ierr /= MPI_SUCCESS) then
+          stop 'Unable to read jBH_ee with MPI'
+        endif
+      IRP_ENDIF
+    endif
+  else
+
+    jBH_ee = 1.d0
+    call ezfio_set_jastrow_jBH_ee(jBH_ee)
+  endif
+
+  ! ---
+
+  if(mpi_master) then
+    call ezfio_has_jastrow_jBH_en(exists)
+  endif
+
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+
+  IRP_IF MPI
+    call MPI_BCAST(jBH_en, (nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+    if(ierr /= MPI_SUCCESS) then
+      stop 'Unable to read Boys-Handy e-n param with MPI'
+    endif
+  IRP_ENDIF
+
+  if(exists) then
+    if(mpi_master) then
+      write(6,'(A)') '.. >>>>> [ IO READ: jBH_en ] <<<<< ..'
+      call ezfio_get_jastrow_jBH_en(jBH_en)
+      IRP_IF MPI
+        call MPI_BCAST(jBH_en, (nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+        if (ierr /= MPI_SUCCESS) then
+          stop 'Unable to read jBH_en with MPI'
+        endif
+      IRP_ENDIF
+    endif
+  else
+
+    jBH_en = 1.d0
+    call ezfio_set_jastrow_jBH_en(jBH_en)
+  endif
+
+  ! ---
+
+  if(mpi_master) then
+    call ezfio_has_jastrow_jBH_c(exists)
+  endif
+
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+
+  IRP_IF MPI
+    call MPI_BCAST(jBH_c, (jBH_size*nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+    if(ierr /= MPI_SUCCESS) then
+      stop 'Unable to read Boys-Handy coeff with MPI'
+    endif
+  IRP_ENDIF
+
+  if(exists) then
+    if(mpi_master) then
+      write(6,'(A)') '.. >>>>> [ IO READ: jBH_c ] <<<<< ..'
+      call ezfio_get_jastrow_jBH_c(jBH_c)
+      IRP_IF MPI
+        call MPI_BCAST(jBH_c, (jBH_size*nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+        if(ierr /= MPI_SUCCESS) then
+          stop 'Unable to read jBH_c with MPI'
+        endif
+      IRP_ENDIF
+    endif
+  else
+
+    jBH_c = 0.d0
+    call ezfio_set_jastrow_jBH_c(jBH_c)
+  endif
+
+  ! ---
+
+  if(mpi_master) then
+    call ezfio_has_jastrow_jBH_m(exists)
+  endif
+
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+
+  IRP_IF MPI
+    call MPI_BCAST(jBH_m, (jBH_size*nucl_num), MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+    if(ierr /= MPI_SUCCESS) then
+      stop 'Unable to read Boys-Handy m powers with MPI'
+    endif
+  IRP_ENDIF
+
+  if(exists) then
+    if(mpi_master) then
+      write(6,'(A)') '.. >>>>> [ IO READ: jBH_m ] <<<<< ..'
+      call ezfio_get_jastrow_jBH_m(jBH_m)
+      IRP_IF MPI
+        call MPI_BCAST(jBH_m, (jBH_size*nucl_num), MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+        if(ierr /= MPI_SUCCESS) then
+          stop 'Unable to read jBH_m with MPI'
+        endif
+      IRP_ENDIF
+    endif
+  else
+
+    jBH_m = 0
+    call ezfio_set_jastrow_jBH_m(jBH_m)
+  endif
+
+  ! ---
+
+  if(mpi_master) then
+    call ezfio_has_jastrow_jBH_n(exists)
+  endif
+
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+
+  IRP_IF MPI
+    call MPI_BCAST(jBH_n, (jBH_size*nucl_num), MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+    if(ierr /= MPI_SUCCESS) then
+      stop 'Unable to read Boys-Handy n powers with MPI'
+    endif
+  IRP_ENDIF
+
+  if(exists) then
+    if(mpi_master) then
+      write(6,'(A)') '.. >>>>> [ IO READ: jBH_n ] <<<<< ..'
+      call ezfio_get_jastrow_jBH_n(jBH_n)
+      IRP_IF MPI
+        call MPI_BCAST(jBH_n, (jBH_size*nucl_num), MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+        if(ierr /= MPI_SUCCESS) then
+          stop 'Unable to read jBH_n with MPI'
+        endif
+      IRP_ENDIF
+    endif
+  else
+
+    jBH_n = 0
+    call ezfio_set_jastrow_jBH_n(jBH_n)
+  endif
+
+  ! ---
+
+  if(mpi_master) then
+    call ezfio_has_jastrow_jBH_o(exists)
+  endif
+
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+
+  IRP_IF MPI
+    call MPI_BCAST(jBH_o, (jBH_size*nucl_num), MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+    if(ierr /= MPI_SUCCESS) then
+      stop 'Unable to read Boys-Handy o powers with MPI'
+    endif
+  IRP_ENDIF
+
+  if(exists) then
+    if(mpi_master) then
+      write(6,'(A)') '.. >>>>> [ IO READ: jBH_o ] <<<<< ..'
+      call ezfio_get_jastrow_jBH_o(jBH_o)
+      IRP_IF MPI
+        call MPI_BCAST(jBH_o, (jBH_size*nucl_num), MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+        if(ierr /= MPI_SUCCESS) then
+          stop 'Unable to read jBH_o with MPI'
+        endif
+      IRP_ENDIF
+    endif
+  else
+
+    jBH_o = 0
+    call ezfio_set_jastrow_jBH_o(jBH_o)
+  endif
+
+  ! ---
+
+  print *, ' parameters for Boys-Handy Jastrow'
+  print *, ' nb of terms per nucleus = ', jBH_size
+
+  do i_nucl = 1, nucl_num
+    print *, ' i_nucl  = ', i_nucl
+    print *, ' ee-term = ', jBH_ee(i_nucl)
+    print *, ' en-term = ', jBH_en(i_nucl)
+    print *, 'm     n     o     c'
+    do p = 1, jBH_size
+      write(*,'(3(I4,2x), E15.7)') jBH_m(p,i_nucl), jBH_n(p,i_nucl), jBH_o(p,i_nucl), jBH_c(p,i_nucl)
+    enddo
+  enddo
+
+
+END_PROVIDER
+
+! ---
+

plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f

@@ -109,6 +109,16 @@ subroutine get_grad1_u12_withsq_r1_seq(ipoint, n_grid2, resx, resy, resz, res)
 
     endif ! env_type
 
+  elseif(j2e_type .eq. "Boys_Handy") then
+
+    PROVIDE jBH_size jBH_en jBH_ee jBH_m jBH_n jBH_o jBH_c
+
+    if(env_type .ne. "None") then
+
+      call grad1_j12_r1_seq(r1, n_grid2, resx, resy, resz)
+
+    endif ! env_type
+
   else
 
     print *, ' Error in get_grad1_u12_withsq_r1_seq: Unknown Jastrow'
@@ -157,9 +167,13 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz)
   double precision, intent(out) :: gradz(n_grid2)
 
   integer                       :: jpoint
+  integer                       :: i_nucl, p, mpA, npA, opA
   double precision              :: r2(3)
   double precision              :: dx, dy, dz, r12, tmp
   double precision              :: mu_val, mu_tmp, mu_der(3)
+  double precision              :: rn(3), f1A, gard1_f1A(3), f2A, gard2_f2A(3), g12, gard1_g12(3)
+  double precision              :: tmp1, tmp2
+
 
   PROVIDE j2e_type
 
@@ -267,6 +281,57 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz)
       gradz(jpoint) = tmp * dz
     enddo
 
+  elseif(j2e_type .eq. "Boys_Handy") then
+
+    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)
+
+      gradx(jpoint) = 0.d0 
+      grady(jpoint) = 0.d0 
+      gradz(jpoint) = 0.d0 
+      do i_nucl = 1, nucl_num 
+
+        rn(1) = nucl_coord(i_nucl,1)
+        rn(2) = nucl_coord(i_nucl,2)
+        rn(3) = nucl_coord(i_nucl,3)
+
+        call jBH_elem_fct_grad(jBH_en(i_nucl), r1, rn, f1A, gard1_f1A)
+        call jBH_elem_fct_grad(jBH_en(i_nucl), r2, rn, f2A, gard2_f2A)
+        call jBH_elem_fct_grad(jBH_ee(i_nucl), r1, r2, g12, gard1_g12)
+
+        do p = 1, jBH_size
+          mpA = jBH_m(p,i_nucl)
+          npA = jBH_n(p,i_nucl)
+          opA = jBH_o(p,i_nucl)
+          tmp = jBH_c(p,i_nucl)
+          if(mpA .eq. npA) then
+            tmp = tmp * 0.5d0
+          endif
+
+          tmp1 = 0.d0
+          if(mpA .gt. 0) then
+            tmp1 = tmp1 + dble(mpA) * f1A**dble(mpA-1) * f2A**dble(npA)
+          endif
+          if(npA .gt. 0) then
+            tmp1 = tmp1 + dble(npA) * f1A**dble(npA-1) * f2A**dble(mpA)
+          endif
+          tmp1 = tmp1 * g12**dble(opA)
+
+          tmp2 = 0.d0
+          if(opA .gt. 0) then
+            tmp2 = tmp2 + dble(opA) * g12**dble(opA-1) * (f1A**dble(mpA) * f2A**dble(npA) + f1A**dble(npA) * f2A**dble(mpA))
+          endif
+
+          gradx(jpoint) = gradx(jpoint) + tmp * (tmp1 * gard1_f1A(1) + tmp2 * gard1_g12(1))
+          grady(jpoint) = grady(jpoint) + tmp * (tmp1 * gard1_f1A(2) + tmp2 * gard1_g12(2))
+          gradz(jpoint) = gradz(jpoint) + tmp * (tmp1 * gard1_f1A(3) + tmp2 * gard1_g12(3))
+        enddo ! p
+      enddo ! i_nucl
+    enddo ! jpoint
+
   else
 
     print *, ' Error in grad1_j12_r1_seq: Unknown j2e_type = ', j2e_type
@@ -757,3 +822,34 @@ end
 
 ! ---
 
+subroutine jBH_elem_fct_grad(alpha, r1, r2, fct, gard1_fct)
+
+  implicit none
+  double precision, intent(in)  :: alpha, r1(3), r2(3)
+  double precision, intent(out) :: fct, gard1_fct(3)
+  double precision              :: dist, tmp1, tmp2
+
+  dist = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) &
+              + (r1(2) - r2(2)) * (r1(2) - r2(2)) &
+              + (r1(3) - r2(3)) * (r1(3) - r2(3)) )
+
+  tmp1 = 1.d0 / (1.d0 + alpha * dist)
+
+  fct = alpha * dist * tmp1
+
+  if(dist .lt. 1d-10) then
+    gard1_fct(1) = 0.d0
+    gard1_fct(2) = 0.d0
+    gard1_fct(3) = 0.d0
+  else
+    tmp2 = alpha * tmp1 * tmp1 / dist
+    gard1_fct(1) = tmp2 * (r1(1) - r2(1))
+    gard1_fct(2) = tmp2 * (r1(2) - r2(2))
+    gard1_fct(3) = tmp2 * (r1(3) - r2(3))
+  endif
+
+  return
+end 
+
+! ---
+