From 869efdd0772cb1bf7b94463023cd5e2974a04ed2 Mon Sep 17 00:00:00 2001 From: eginer Date: Sun, 8 Sep 2024 19:14:15 +0200 Subject: [PATCH] new jpsiking is working on He --- plugins/local/non_h_ints_mu/deb_j_psi.irp.f | 16 ++++---- plugins/local/non_h_ints_mu/jastrow_psi.irp.f | 40 ++++++++++--------- 2 files changed, 29 insertions(+), 27 deletions(-) diff --git a/plugins/local/non_h_ints_mu/deb_j_psi.irp.f b/plugins/local/non_h_ints_mu/deb_j_psi.irp.f index 248d2f31..1b034684 100644 --- a/plugins/local/non_h_ints_mu/deb_j_psi.irp.f +++ b/plugins/local/non_h_ints_mu/deb_j_psi.irp.f @@ -1,7 +1,7 @@ program test_j_mu_of_r implicit none -! call routine_deb_j_psi - call routine_deb_denom + call routine_deb_j_psi +! call routine_deb_denom end subroutine routine_deb_j_psi @@ -9,7 +9,7 @@ subroutine routine_deb_j_psi 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,grad_jast(3),j_bump,jastrow_psi,grad_jast_bis(3) double precision :: jast_p,jast_m,num_grad_jast(3) dr = 0.00001d0 @@ -28,11 +28,11 @@ subroutine routine_deb_j_psi do k = 1, 3 r1bis= r1 r1bis(k) += dr - jast_p = jastrow_psi(r1bis, r2) + call get_grad_r1_jastrow_psi(r1bis,r2,grad_jast_bis,jast_p) r1bis= r1 r1bis(k) -= dr - jast_m = jastrow_psi(r1bis, r2) + call get_grad_r1_jastrow_psi(r1bis,r2,grad_jast_bis,jast_m) num_grad_jast(k) = (jast_p - jast_m)/(2.d0* dr) norm += num_grad_jast(k)*num_grad_jast(k) @@ -86,19 +86,19 @@ subroutine routine_deb_denom 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) + call denom_jpsi(i,j,a_boys, 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) + call denom_jpsi(i,j,a_boys, 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) + call denom_jpsi(i,j,a_boys, 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) diff --git a/plugins/local/non_h_ints_mu/jastrow_psi.irp.f b/plugins/local/non_h_ints_mu/jastrow_psi.irp.f index a3f0b524..4c88e793 100644 --- a/plugins/local/non_h_ints_mu/jastrow_psi.irp.f +++ b/plugins/local/non_h_ints_mu/jastrow_psi.irp.f @@ -56,29 +56,25 @@ subroutine get_grad_r1_jastrow_psi(r1,r2,grad_j_psi_r1,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(:,:) + double precision :: num_j, denom_j, num_j_grad(3), denom_j_grad(3),delta,coef + double precision :: inv_denom_j 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)) + delta = a_boys 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 + call denom_jpsi(i,j,delta,mos_array_r1,mos_grad_array_r1,mos_array_r2,denom_j, denom_j_grad) + inv_denom_j = 1.d0/denom_j do a = elec_beta_num+1, mo_num ! r1 do b = elec_alpha_num+1, mo_num ! r2 + call numerator_psi(a,b,mos_array_r1,mos_grad_array_r1,mos_array_r2,num_j, num_j_grad) 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 + jast += coef * num_j * inv_denom_j + grad_j_psi_r1 += coef * (num_j_grad * denom_j - num_j * denom_j_grad) * inv_denom_j * inv_denom_j enddo enddo enddo @@ -100,14 +96,12 @@ subroutine get_grad_r1_jastrow_psi(r1,r2,grad_j_psi_r1,jast) end -subroutine denom_jpsi(i,j,mos_array_r1,mos_grad_array_r1,mos_array_r2,denom, grad_denom) +subroutine denom_jpsi(i,j,delta,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(in) :: mos_array_r1(mo_num),mos_grad_array_r1(3,mo_num),mos_array_r2(mo_num),delta 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 @@ -116,7 +110,15 @@ subroutine denom_jpsi(i,j,mos_array_r1,mos_grad_array_r1,mos_array_r2,denom, gra 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) - + denom = phi_i_phi_j + delta * inv_phi_i_phi_j + grad_denom(:) = (1.d0 - delta*inv_phi_i_phi_j_2) * mos_array_r2(j) * mos_grad_array_r1(:,i) +end + +subroutine numerator_psi(a,b,mos_array_r1,mos_grad_array_r1,mos_array_r2,num, grad_num) + implicit none + integer, intent(in) :: a,b + 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) :: num, grad_num(3) + num = mos_array_r1(a) * mos_array_r2(b) + grad_num(:) = mos_array_r2(b) * mos_grad_array_r1(:,a) end