starting jpsiking

plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f

@@ -155,6 +155,7 @@ subroutine run_stochastic_cipsi
   call pt2_alloc(pt2_data_err, N_states)
   call ZMQ_pt2(E_tc, pt2_data, pt2_data_err, relative_error,0) ! Stochastic PT2 and selection
   call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm)
+  call print_summary_tc(psi_energy_with_nucl_rep, pt2_data, pt2_data_err, N_det, N_configuration, N_states, psi_s2)
   call pt2_dealloc(pt2_data)
   call pt2_dealloc(pt2_data_err)
 

plugins/local/jastrow/EZFIO.cfg

@@ -1,3 +1,10 @@
+[log_jpsi]
+type: logical
+doc: If |true|, the Jpsi is taken as log(1+psi_cor)
+interface: ezfio,provider,ocaml
+default: False
+
+
 [mu_of_r_tc]
 type: character*(32)
 doc: type of the mu(r): [ Standard | Erfmu | Erfmugauss ]

plugins/local/non_h_ints_mu/deb_j_psi.irp.f

@@ -0,0 +1,131 @@
+program test_j_mu_of_r
+ implicit none
+! call routine_deb_j_psi
+ call routine_deb_denom
+end
+
+subroutine routine_deb_j_psi
+ implicit none
+ integer :: ipoint,k 
+ double precision :: r2(3), weight, dr, r1(3), r1bis(3)
+ double precision :: accu_grad(3)
+ double precision :: jast,grad_jast(3),j_bump,jastrow_psi
+ double precision :: jast_p,jast_m,num_grad_jast(3)
+
+ dr = 0.00001d0
+ r2 = 0.d0
+ r2(1) =  0.5d0
+ r2(2) = -0.1d0
+ r2(3) =  1.0d0
+ accu_grad = 0.d0
+  do ipoint = 1, n_points_final_grid
+   r1(1:3) = final_grid_points(1:3,ipoint)
+   weight = final_weight_at_r_vector(ipoint)
+   call get_grad_r1_jastrow_psi(r1,r2,grad_jast,jast)
+!   grad_jast = - grad_jast
+   double precision :: norm,error
+   norm = 0.D0
+   do k = 1, 3
+    r1bis= r1
+    r1bis(k) += dr
+    jast_p = jastrow_psi(r1bis, r2)
+
+    r1bis= r1
+    r1bis(k) -= dr
+    jast_m = jastrow_psi(r1bis, r2)
+
+    num_grad_jast(k) = (jast_p - jast_m)/(2.d0* dr)
+    norm += num_grad_jast(k)*num_grad_jast(k)
+   enddo
+   error = 0.d0
+   do k = 1, 3
+    error += dabs(grad_jast(k) - num_grad_jast(k))
+   enddo
+   error *= 0.33333333d0
+   norm = dsqrt(norm)
+   if(norm.gt.1.d-05)then
+    if(dabs(error/norm).gt.dr)then
+     print*,'/////'
+     print*,error,norm
+     print*,grad_jast
+     print*,num_grad_jast
+    endif
+   endif
+   do k = 1,3
+    accu_grad(k) += weight * dabs(grad_jast(k) - num_grad_jast(k))
+   enddo
+  enddo
+ print*,'accu_grad = '
+ print*, accu_grad
+
+end
+
+
+subroutine routine_deb_denom
+ implicit none
+ integer :: ipoint,k,i,j
+ double precision :: r2(3), weight, dr, r1(3), r1bis(3)
+ double precision :: accu_grad(3)
+ double precision :: jast,grad_jast(3),j_bump,jastrow_psi,grad_jast_bis(3)
+ double precision :: jast_p,jast_m,num_grad_jast(3)
+
+ dr = 0.00001d0
+ r2 = 0.d0
+ r2(1) =  0.5d0
+ r2(2) = -0.1d0
+ r2(3) =  1.0d0
+ double precision, allocatable :: mos_array_r1(:), mos_array_r2(:)
+ double precision, allocatable :: mos_grad_array_r1(:,:),mos_grad_array_r2(:,:)
+ allocate(mos_array_r1(mo_num), mos_array_r2(mo_num))
+ allocate(mos_grad_array_r1(3,mo_num), mos_grad_array_r2(3,mo_num))
+ do i = 1, 1
+  do j = 1, 1
+   accu_grad = 0.d0
+   call give_all_mos_and_grad_at_r(r2,mos_array_r2,mos_grad_array_r2)
+    do ipoint = 1, n_points_final_grid
+     r1(1:3) = final_grid_points(1:3,ipoint)
+     weight = final_weight_at_r_vector(ipoint)
+     call give_all_mos_and_grad_at_r(r1,mos_array_r1,mos_grad_array_r1)
+     call denom_jpsi(i,j,mos_array_r1,mos_grad_array_r1,mos_array_r2,jast, grad_jast)
+     double precision :: norm,error
+     norm = 0.D0
+     do k = 1, 3
+      r1bis= r1
+      r1bis(k) += dr
+      call give_all_mos_and_grad_at_r(r1bis,mos_array_r1,mos_grad_array_r1)
+      call denom_jpsi(i,j,mos_array_r1,mos_grad_array_r1,mos_array_r2,jast_p, grad_jast_bis)
+  
+      r1bis= r1
+      r1bis(k) -= dr
+      call give_all_mos_and_grad_at_r(r1bis,mos_array_r1,mos_grad_array_r1)
+      call denom_jpsi(i,j,mos_array_r1,mos_grad_array_r1,mos_array_r2,jast_m, grad_jast_bis)
+  
+      num_grad_jast(k) = (jast_p - jast_m)/(2.d0* dr)
+      norm += num_grad_jast(k)*num_grad_jast(k)
+     enddo
+     error = 0.d0
+     do k = 1, 3
+      error += dabs(grad_jast(k) - num_grad_jast(k))
+     enddo
+     error *= 0.33333333d0
+     norm = dsqrt(norm)
+     if(norm.gt.1.d-05)then
+      if(dabs(error/norm).gt.dr)then
+       print*,'/////'
+       print*,error,norm
+       print*,grad_jast
+       print*,num_grad_jast
+      endif
+     endif
+     do k = 1,3
+      accu_grad(k) += weight * dabs(grad_jast(k) - num_grad_jast(k))
+     enddo
+    enddo
+   print*,'i,j = ',i,j
+   print*,'accu_grad = '
+   print*, accu_grad
+  enddo
+ enddo
+
+end
+

plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f

@@ -35,6 +35,7 @@ subroutine get_grad1_u12_withsq_r1_seq(ipoint, n_grid2, resx, resy, resz, res)
 
   if( (j2e_type .eq. "Mu")        .or. &
       (j2e_type .eq. "Mur")       .or. &
+      (j2e_type .eq. "Jpsi")   .or. &
       (j2e_type .eq. "Mugauss")   .or. &
       (j2e_type .eq. "Murgauss")  .or. &
       (j2e_type .eq. "Bump")      .or. &
@@ -417,6 +418,17 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz)
       enddo ! i_nucl
     enddo ! jpoint
 
+  elseif(j2e_type .eq. "Jpsi") then
+   double precision :: grad_j_psi_r1(3),jast_psi
+   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 get_grad_r1_jastrow_psi(r1,r2,grad_j_psi_r1,jast_psi)
+     gradx(jpoint) = grad_j_psi_r1(1)
+     grady(jpoint) = grad_j_psi_r1(2)
+     gradz(jpoint) = grad_j_psi_r1(3)
+   enddo
   else
 
     print *, ' Error in grad1_j12_r1_seq: Unknown j2e_type = ', j2e_type
@@ -723,6 +735,7 @@ subroutine get_grad1_u12_2e_r1_seq(ipoint, n_grid2, resx, resy, resz)
   if( (j2e_type .eq. "Mu")        .or. &
       (j2e_type .eq. "Mugauss")   .or. &
       (j2e_type .eq. "Mur")       .or. &
+      (j2e_type .eq. "Jpsi")       .or. &
       (j2e_type .eq. "Murgauss")  .or. &
       (j2e_type .eq. "Bump")      .or. &
       (j2e_type .eq. "Boys") ) then

plugins/local/non_h_ints_mu/jastrow_psi.irp.f

@@ -0,0 +1,122 @@
+BEGIN_PROVIDER [ double precision, c_ij_ab_jastrow, (mo_num, mo_num, elec_alpha_num, elec_beta_num)]
+ implicit none
+ integer :: iunit, getUnitAndOpen
+ c_ij_ab_jastrow = 0.d0
+ iunit = getUnitAndOpen(trim(ezfio_work_dir)//'c_ij_ab', 'R')                                                     
+ read(iunit) c_ij_ab_jastrow
+ close(iunit)
+ print*,'c_ij_ab_jastrow = '
+ integer :: i,j,a,b
+ do i = 1, elec_beta_num ! r2
+  do j = 1, elec_alpha_num ! r1
+   do a = elec_beta_num+1, mo_num ! r2
+    do b = elec_alpha_num+1, mo_num ! r1 
+!     print*,b,a,j,i
+     print*,c_ij_ab_jastrow(b,a,j,i),b,a,j,i
+     if(dabs(c_ij_ab_jastrow(b,a,j,i)).lt.1.d-12)then
+      c_ij_ab_jastrow(b,a,j,i) = 0.d0
+     endif
+    enddo
+   enddo
+  enddo
+ enddo
+END_PROVIDER 
+
+double precision function jastrow_psi(r1,r2)
+ implicit none
+ double precision, intent(in) :: r1(3), r2(3)
+ integer :: i,j,a,b
+ double precision, allocatable :: mos_array_r1(:), mos_array_r2(:)
+ allocate(mos_array_r1(mo_num), mos_array_r2(mo_num))
+ call give_all_mos_at_r(r1,mos_array_r1)
+ call give_all_mos_at_r(r2,mos_array_r2)
+ double precision :: eps,coef, numerator,denominator
+ double precision :: phi_i_phi_j
+ eps = a_boys
+ jastrow_psi= 0.d0
+ do i = 1, elec_beta_num ! r1
+  do j = 1, elec_alpha_num ! r2
+   phi_i_phi_j = mos_array_r1(i) * mos_array_r2(j) + eps
+   denominator = 1.d0/phi_i_phi_j
+   do a = elec_beta_num+1, mo_num ! r1
+    do b = elec_alpha_num+1, mo_num ! r2 
+     coef = c_ij_ab_jastrow(b,a,j,i)
+     numerator = mos_array_r2(b) * mos_array_r1(a)
+     jastrow_psi += coef * numerator*denominator
+    enddo
+   enddo
+  enddo
+ enddo
+end
+
+subroutine get_grad_r1_jastrow_psi(r1,r2,grad_j_psi_r1,jast)
+ implicit none
+ double precision, intent(in) :: r1(3), r2(3)
+ double precision, intent(out):: grad_j_psi_r1(3),jast
+ integer :: i,j,a,b
+ double precision, allocatable :: mos_array_r1(:), mos_array_r2(:)
+ double precision, allocatable :: mos_grad_array_r1(:,:),mos_grad_array_r2(:,:)
+ allocate(mos_array_r1(mo_num), mos_array_r2(mo_num))
+ allocate(mos_grad_array_r1(3,mo_num), mos_grad_array_r2(3,mo_num))
+ call give_all_mos_and_grad_at_r(r1,mos_array_r1,mos_grad_array_r1)
+ call give_all_mos_and_grad_at_r(r2,mos_array_r2,mos_grad_array_r2)
+ double precision :: eps,coef, numerator(3),denominator
+ double precision :: phi_i_phi_j
+ eps = a_boys
+ grad_j_psi_r1 = 0.d0
+ jast = 0.d0
+ do i = 1, elec_beta_num ! r1
+  do j = 1, elec_alpha_num ! r2
+   phi_i_phi_j = mos_array_r1(i) * mos_array_r2(j) + eps
+   denominator = 1.d0/phi_i_phi_j
+   denominator *= denominator
+   do a = elec_beta_num+1, mo_num ! r1
+    do b = elec_alpha_num+1, mo_num ! r2 
+     coef = c_ij_ab_jastrow(b,a,j,i)
+     jast += phi_i_phi_j * mos_array_r2(b) * mos_array_r1(a) * coef
+!     print*,b,a,j,i,c_ij_ab_jastrow(b,a,j,i)
+!     print*,'jast = ',jast
+     numerator(:) = mos_array_r2(b) * mos_grad_array_r1(:,a) & 
+      * phi_i_phi_j - mos_array_r1(a) * mos_array_r2(b) * mos_array_r2(j) * mos_grad_array_r1(:,i)
+     grad_j_psi_r1 += coef * numerator*denominator
+    enddo
+   enddo
+  enddo
+ enddo
+ if(jast.lt.-1.d0.or.dabs(jast).gt.1.d0)then
+  print*,'pb ! '
+  print*,jast
+  print*,dsqrt(r1(1)**2+r1(2)**2+r1(3)**2),dsqrt(r2(1)**2+r2(2)**2+r2(3)**2)
+  print*,r1
+!  print*,mos_array_r1(1:2)
+  print*,r2
+!  print*,mos_array_r2(1:2)
+  stop
+ endif
+ if(log_jpsi)then
+  grad_j_psi_r1 = grad_j_psi_r1/(1.d0 + jast)
+ endif
+
+end
+
+
+subroutine denom_jpsi(i,j,mos_array_r1,mos_grad_array_r1,mos_array_r2,denom, grad_denom)
+ implicit none
+ integer, intent(in)           :: i,j
+ double precision, intent(in)  :: mos_array_r1(mo_num),mos_grad_array_r1(3,mo_num),mos_array_r2(mo_num)
+ double precision, intent(out) :: denom, grad_denom(3)
+ double precision :: coef,phi_i_phi_j,inv_phi_i_phi_j,inv_phi_i_phi_j_2
+ denom = 0.d0
+ grad_denom = 0.d0
+ phi_i_phi_j = mos_array_r1(i) * mos_array_r2(j)
+ if(phi_i_phi_j /= 0.d0)then
+  inv_phi_i_phi_j   = 1.d0/phi_i_phi_j
+  inv_phi_i_phi_j_2 = 1.d0/(phi_i_phi_j * phi_i_phi_j)
+ else
+  inv_phi_i_phi_j   = huge(1.0)
+  inv_phi_i_phi_j_2 = huge(1.d0)
+ endif
+ denom += phi_i_phi_j + a_boys * inv_phi_i_phi_j
+ grad_denom(:) += (1.d0 - a_boys*inv_phi_i_phi_j_2) * mos_array_r2(j) * mos_grad_array_r1(:,i)
+
+end

plugins/local/non_h_ints_mu/plot_mo.irp.f

@@ -0,0 +1,19 @@
+program plot_mo
+ implicit none
+ integer :: i,npt
+ double precision :: xmin,xmax,dx,r(3)
+ double precision,allocatable :: mos_array(:)
+ allocate(mos_array(mo_num))
+ npt = 10000
+ xmin =0.d0
+ xmax =10.d0
+ dx=(xmax-xmin)/dble(npt)
+ r=0.d0
+ r(1) = xmin
+ do i = 1, npt
+  call give_all_mos_at_r(r,mos_array)
+  write(33,'(100(F16.10,X))')r(1),mos_array(1),mos_array(2),mos_array(3)
+  r(1) += dx
+ enddo
+
+end

plugins/local/non_h_ints_mu/print_jastrow_psi.irp.f

@@ -0,0 +1,15 @@
+program print_j_psi
+ implicit none
+ integer :: i,j,a,b
+ do i = 1, elec_beta_num ! r2
+  do j = 1, elec_alpha_num ! r1
+   do a = elec_beta_num+1, mo_num ! r2
+    do b = elec_alpha_num+1, mo_num ! r1 
+     print*,b,a,j,i
+     print*,c_ij_ab_jastrow(b,a,j,i)
+    enddo
+   enddo
+  enddo
+ enddo
+
+end