10
0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 20:34:58 +01:00

failure attempts in mu(r)

This commit is contained in:
eginer 2024-09-05 22:01:22 +02:00
parent 65bef0b266
commit 1221d48395
17 changed files with 877 additions and 120 deletions

View File

@ -1,8 +1,8 @@
[mu_of_r_tc] [mu_of_r_tc]
type: logical type: character*(32)
doc: If |true|, take the second formula for mu(r) doc: type of the mu(r): [ Standard | Erfmu | Erfmugauss ]
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: False default: Standard
[mu_of_r_av] [mu_of_r_av]
type: logical type: logical
@ -12,7 +12,7 @@ default: False
[j2e_type] [j2e_type]
type: character*(32) type: character*(32)
doc: type of the 2e-Jastrow: [ None | Mu | Mu_Nu | Mur | Boys | Boys_Handy | Qmckl ] doc: type of the 2e-Jastrow: [ None | Mu | Mugauss | Mu_Nu | Mur | Murgauss | Bump | Boys | Boys_Handy | Qmckl ]
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: Mu default: Mu

View File

@ -0,0 +1,28 @@
program test_j_mu_of_r
implicit none
double precision :: x,mu_min,dmu,mu_max, mu, mu_p, mu_m
double precision :: j_simple,j_p, j_m,numeric_d_mu,d_dx_mu
double precision :: accu
integer :: npt,i
npt = 1000
mu_min = 0.3d0
mu_max = 10.d0
dmu = (mu_max - mu_min)/dble(npt)
x = 0.7d0
mu = mu_min
do i = 1, npt
call get_deriv_mu_j12(x,mu,d_dx_mu)
mu_p = mu + dmu
mu_m = mu - dmu
j_p = j_simple(x,mu_p)
j_m = j_simple(x,mu_m)
numeric_d_mu = 0.5d0 * (j_p - j_m)/dmu
print*,mu
print*,numeric_d_mu,d_dx_mu,dabs(d_dx_mu-numeric_d_mu)
accu += dabs(d_dx_mu-numeric_d_mu)
mu += dmu
enddo
accu *= dmu
print*,'accu = ',accu
end

View File

@ -0,0 +1,98 @@
program test_j_mu_of_r
implicit none
! call routine_test_mu_of_r
call routine_test_mu_of_r_tot
end
subroutine routine_test_mu_of_r_tot
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_mu_r1(3),j_bump
double precision :: jast_p,jast_m,num_grad_jast_mu_r1(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 grad_j_sum_mu_of_r(r1,r2,jast,grad_jast_mu_r1)
call get_grad_j_bump_mu_of_r(r1,r2,grad_jast_mu_r1)
double precision :: norm,error
norm = 0.D0
do k = 1, 3
r1bis= r1
r1bis(k) += dr
jast_p = j_bump(r1bis,r2,a_boys)
r1bis= r1
r1bis(k) -= dr
jast_m = j_bump(r1bis,r2,a_boys)
num_grad_jast_mu_r1(k) = (jast_p - jast_m)/(2.d0* dr)
norm += num_grad_jast_mu_r1(k)*num_grad_jast_mu_r1(k)
enddo
error = 0.d0
do k = 1, 3
error += dabs(grad_jast_mu_r1(k) - num_grad_jast_mu_r1(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_mu_r1
print*,num_grad_jast_mu_r1
endif
endif
do k = 1,3
accu_grad(k) += weight * dabs(grad_jast_mu_r1(k) - num_grad_jast_mu_r1(k))
enddo
enddo
print*,'accu_grad = '
print*, accu_grad
end
subroutine routine_test_mu_of_r
implicit none
integer :: ipoint,k
double precision :: weight, dr, r1(3), r1bis(3),accu_grad(3),num_grad_mu_r1(3)
double precision :: mu_r1,dm_r1, mu_der_r1(3), grad_dm_r1(3)
double precision :: mu_der_rp(3), grad_dm_rp(3),mu_rp
double precision :: mu_der_rm(3), grad_dm_rm(3),mu_rm
dr = 0.0001d0
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 grad_mu_of_r_mean_field(r1,mu_r1,dm_r1, mu_der_r1, grad_dm_r1)
do k = 1, 3
r1bis= r1
r1bis(k) += dr
call grad_mu_of_r_mean_field(r1bis,mu_rp, dm_r1, mu_der_rp, grad_dm_r1)
r1bis= r1
r1bis(k) -= dr
call grad_mu_of_r_mean_field(r1bis,mu_rm, dm_r1, mu_der_rm, grad_dm_r1)
num_grad_mu_r1(k) = (mu_rp - mu_rm)/(2.d0* dr)
! print*,jast_mu_r1_p,jast_mu_r1_m
enddo
print*,'/////'
print*,mu_der_r1
print*,num_grad_mu_r1
do k = 1,3
accu_grad(k) += weight * dabs(mu_der_r1(k) - num_grad_mu_r1(k))
enddo
enddo
print*,'accu_grad = '
print*, accu_grad
end

View File

@ -0,0 +1,62 @@
program test_j_mu_of_r
implicit none
! call routine_test_mu_of_r
call routine_test_mu_of_r_tot
end
subroutine routine_test_mu_of_r_tot
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,j12_mu
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 grad1_j12_mu(r1, r2, grad_jast)
grad_jast = - grad_jast
double precision :: norm,error
norm = 0.D0
do k = 1, 3
r1bis= r1
r1bis(k) += dr
jast_p = j12_mu(r1bis, r2)
r1bis= r1
r1bis(k) -= dr
jast_m = j12_mu(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

View File

@ -12,7 +12,7 @@ subroutine routine_test_mu_of_r_tot
double precision :: jast,grad_jast_mu_r1(3) double precision :: jast,grad_jast_mu_r1(3)
double precision :: jast_p,jast_m,num_grad_jast_mu_r1(3) double precision :: jast_p,jast_m,num_grad_jast_mu_r1(3)
dr = 0.00001d0 dr = 0.000001d0
r2 = 0.d0 r2 = 0.d0
r2(1) = 0.5d0 r2(1) = 0.5d0
r2(2) = -0.1d0 r2(2) = -0.1d0
@ -22,7 +22,7 @@ subroutine routine_test_mu_of_r_tot
r1(1:3) = final_grid_points(1:3,ipoint) r1(1:3) = final_grid_points(1:3,ipoint)
weight = final_weight_at_r_vector(ipoint) weight = final_weight_at_r_vector(ipoint)
call grad_j_sum_mu_of_r(r1,r2,jast,grad_jast_mu_r1) call grad_j_sum_mu_of_r(r1,r2,jast,grad_jast_mu_r1)
double precision :: norm double precision :: norm,error
norm = 0.D0 norm = 0.D0
do k = 1, 3 do k = 1, 3
r1bis= r1 r1bis= r1
@ -36,12 +36,20 @@ subroutine routine_test_mu_of_r_tot
num_grad_jast_mu_r1(k) = (jast_p - jast_m)/(2.d0* dr) num_grad_jast_mu_r1(k) = (jast_p - jast_m)/(2.d0* dr)
norm += num_grad_jast_mu_r1(k)*num_grad_jast_mu_r1(k) norm += num_grad_jast_mu_r1(k)*num_grad_jast_mu_r1(k)
enddo enddo
error = 0.d0
do k = 1, 3
error += dabs(grad_jast_mu_r1(k) - num_grad_jast_mu_r1(k))
enddo
error *= 0.33333333d0
norm = dsqrt(norm) norm = dsqrt(norm)
if(norm.gt.1.d-10)then if(norm.gt.1.d-05)then
if(dabs(error/norm).gt.10.d0*dr)then
print*,'/////' print*,'/////'
print*,error,norm,dabs(error/norm)
print*,grad_jast_mu_r1 print*,grad_jast_mu_r1
print*,num_grad_jast_mu_r1 print*,num_grad_jast_mu_r1
endif endif
endif
do k = 1,3 do k = 1,3
accu_grad(k) += weight * dabs(grad_jast_mu_r1(k) - num_grad_jast_mu_r1(k)) accu_grad(k) += weight * dabs(grad_jast_mu_r1(k) - num_grad_jast_mu_r1(k))
enddo enddo
@ -54,41 +62,33 @@ end
subroutine routine_test_mu_of_r subroutine routine_test_mu_of_r
implicit none implicit none
integer :: ipoint,k integer :: ipoint,k
double precision :: r2(3), weight, dr, r1(3), r1bis(3) double precision :: weight, dr, r1(3), r1bis(3),accu_grad(3),num_grad_mu_r1(3)
double precision :: accu_grad(3) double precision :: mu_r1,dm_r1, mu_der_r1(3), grad_dm_r1(3)
double precision :: jast_mu_r1,grad_jast_mu_r1(3) double precision :: mu_der_rp(3), grad_dm_rp(3),mu_rp
double precision :: jast_mu_r1_p,jast_mu_r1_m,num_grad_jast_mu_r1(3) double precision :: mu_der_rm(3), grad_dm_rm(3),mu_rm
double precision :: j12_mu_input
double precision :: mu_val_p, mu_val_m, dm_r1, mu_der(3), grad_dm_r1(3)
dr = 0.0001d0 dr = 0.0001d0
r2 = 0.d0
r2(1) = 2.5d0
r2(2) = 2.25d0
r2(3) = -2.5d0
accu_grad = 0.d0 accu_grad = 0.d0
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
r1(1:3) = final_grid_points(1:3,ipoint) r1(1:3) = final_grid_points(1:3,ipoint)
weight = final_weight_at_r_vector(ipoint) weight = final_weight_at_r_vector(ipoint)
call grad_j_mu_of_r_1(r1,r2,jast_mu_r1,grad_jast_mu_r1, dm_r1, grad_dm_r1) call grad_mu_of_r_mean_field(r1,mu_r1,dm_r1, mu_der_r1, grad_dm_r1)
do k = 1, 3 do k = 1, 3
r1bis= r1 r1bis= r1
r1bis(k) += dr r1bis(k) += dr
call grad_mu_of_r_mean_field(r1bis,mu_val_p, dm_r1, mu_der, grad_dm_r1) call grad_mu_of_r_mean_field(r1bis,mu_rp, dm_r1, mu_der_rp, grad_dm_r1)
jast_mu_r1_p = j12_mu_input(r1bis, r2, mu_val_p)
r1bis= r1 r1bis= r1
r1bis(k) -= dr r1bis(k) -= dr
call grad_mu_of_r_mean_field(r1bis,mu_val_m, dm_r1, mu_der, grad_dm_r1) call grad_mu_of_r_mean_field(r1bis,mu_rm, dm_r1, mu_der_rm, grad_dm_r1)
jast_mu_r1_m = j12_mu_input(r1bis, r2, mu_val_m) num_grad_mu_r1(k) = (mu_rp - mu_rm)/(2.d0* dr)
num_grad_jast_mu_r1(k) = -(jast_mu_r1_p - jast_mu_r1_m)/(2.d0* dr)
! print*,jast_mu_r1_p,jast_mu_r1_m ! print*,jast_mu_r1_p,jast_mu_r1_m
enddo enddo
print*,'/////' print*,'/////'
print*,grad_jast_mu_r1 print*,mu_der_r1
print*,num_grad_jast_mu_r1 print*,num_grad_mu_r1
do k = 1,3 do k = 1,3
accu_grad(k) += weight * dabs(grad_jast_mu_r1(k) - num_grad_jast_mu_r1(k)) accu_grad(k) += weight * dabs(mu_der_r1(k) - num_grad_mu_r1(k))
enddo enddo
enddo enddo
print*,'accu_grad = ' print*,'accu_grad = '

View File

@ -0,0 +1,90 @@
double precision function wigner_radius(rho)
implicit none
include 'constants.include.F'
double precision, intent(in) :: rho
wigner_radius = 4.d0 * pi * rho * 0.333333333333d0
wigner_radius = wigner_radius**(-0.3333333d0)
end
double precision function j_bump(r1,r2,a)
implicit none
include 'constants.include.F'
double precision, intent(in) :: r1(3),r2(3),a
double precision :: inv_a,factor,x_scaled,scalar
double precision :: r12
r12 = (r1(1) - r2(1))*(r1(1) - r2(1))
r12 += (r1(2) - r2(2))*(r1(2) - r2(2))
r12 += (r1(3) - r2(3))*(r1(3) - r2(3))
r12 = dsqrt(r12)
inv_a = 1.d0/a
x_scaled = r12*inv_a*inv_sq_pi
x_scaled*= x_scaled
j_bump = 0.5d0 * (r12-a) * dexp(-x_scaled)
end
subroutine get_grad_j_bump(x,a,grad)
implicit none
BEGIN_DOC
! gradient of the Jastrow with a bump
!
! j(x,a) = 1/2 * (x-a)* exp[-(x/(a*sqrt(pi)))^2]
!
! d/dx j(x,a) = 1/(2 pi a^2) * exp[-(x/(a*sqrt(pi)))^2] * (pi a^2 + 2 a x - 2x^2)
END_DOC
include 'constants.include.F'
double precision, intent(in) :: x,a
double precision, intent(out) :: grad
double precision :: inv_a,factor,x_scaled,scalar
inv_a = 1.d0/a
factor = 0.5d0*inv_pi*inv_a*inv_a
x_scaled = x*inv_a*inv_sq_pi
x_scaled*= x_scaled
grad = factor * dexp(-x_scaled) * (pi*a*a + 2.d0 * a*x - 2.d0*x*x)
end
subroutine get_d_da_j_bump(x,a,d_da)
implicit none
BEGIN_DOC
! Derivative with respect by to the parameter "a" of the Jastrow with a bump
!
! j(x,a) = 1/2 * (x-a)* exp[-(x/(a*sqrt(pi)))^2]
!
! d/da j(x,a) = - 1/(pi*a^3) * exp[-(x/(a*sqrt(pi)))^2] * (-2 x^3 + 2 a x^2 + pi a^x3)
END_DOC
include 'constants.include.F'
double precision, intent(in) :: x,a
double precision, intent(out) :: d_da
double precision :: factor, inv_a,x_scaled,scalar
inv_a = 1.d0/a
factor = inv_a*inv_a*inv_a*inv_pi
x_scaled = x*inv_a*inv_sq_pi
x_scaled*= x_scaled
d_da = factor * dexp(-x_scaled) * (-2.d0 * x*x*x + 2.d0*x*x*a+pi*a*a*a)
end
subroutine get_grad_j_bump_mu_of_r(r1,r2,grad_j_bump)
implicit none
BEGIN_DOC
! d/dx1 j(x,a(r1,r2)) where j(x,a) is the Jastrow with a bump
!
! j(x,a) = 1/2 * (x-a)* exp[-(x/(a*sqrt(pi)))^2]
!
! a(r1,r2) = [rho(r1) a(r1) + rho(r2) a(r2)]/[rho(r1) + rho(r2)]
!
! d/dx1 j(x,a) = d/dx1 j(x,a(r1,r2))
END_DOC
double precision, intent(in) :: r1(3),r2(3)
double precision, intent(out):: grad_j_bump(3)
double precision :: r12,r12_vec(3),grad_scal,inv_r12
r12_vec = r1 - r2
r12 = (r1(1) - r2(1))*(r1(1) - r2(1))
r12 += (r1(2) - r2(2))*(r1(2) - r2(2))
r12 += (r1(3) - r2(3))*(r1(3) - r2(3))
r12 = dsqrt(r12)
call get_grad_j_bump(r12,a_boys,grad_scal)
if(r12.lt.1.d-10)then
grad_j_bump = 0.d0
else
grad_j_bump = grad_scal * r12_vec/r12
endif
end

View File

@ -31,7 +31,7 @@
grad1_u12_squared_num = 0.d0 grad1_u12_squared_num = 0.d0
if( ((j2e_type .eq. "Mu") .and. (env_type .eq. "None")) .or. & if( ((j2e_type .eq. "Mu") .and. (env_type .eq. "None")) .or. &
(j2e_type .eq. "Mur") ) then (j2e_type .eq. "Mur").or.(j2e_type .eq. "Mugauss") .or. (j2e_type .eq. "Murgauss")) then
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &

View File

@ -4,19 +4,48 @@ subroutine get_j_sum_mu_of_r(r1,r2,jast)
double precision, intent(out):: jast double precision, intent(out):: jast
double precision :: mu_r1, dm_r1, grad_mu_r1(3), grad_dm_r1(3), j_mu_r1 double precision :: mu_r1, dm_r1, grad_mu_r1(3), grad_dm_r1(3), j_mu_r1
double precision :: mu_r2, dm_r2, grad_mu_r2(3), grad_dm_r2(3), j_mu_r2 double precision :: mu_r2, dm_r2, grad_mu_r2(3), grad_dm_r2(3), j_mu_r2
double precision :: j12_mu_input double precision :: j12_mu_input,mu_tot,r12,j_simple
jast = 0.d0 jast = 0.d0
if(murho_type==0)then
! J(r1,r2) = [rho(r1) * j(mu(r1),r12) + rho(r2) * j(mu(r2),r12)] / [rho(r1) + rho(r2)]
call grad_mu_of_r_mean_field(r1,mu_r1, dm_r1, grad_mu_r1, grad_dm_r1) call grad_mu_of_r_mean_field(r1,mu_r1, dm_r1, grad_mu_r1, grad_dm_r1)
call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2) call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2)
j_mu_r1 = j12_mu_input(r1, r2, mu_r1) j_mu_r1 = j12_mu_input(r1, r2, mu_r1)
j_mu_r2 = j12_mu_input(r1, r2, mu_r2) j_mu_r2 = j12_mu_input(r1, r2, mu_r2)
if(dm_r1 + dm_r2.lt.1.d-7)return if(dm_r1 + dm_r2.lt.1.d-7)return
jast = (dm_r1 * j_mu_r1 + dm_r2 * j_mu_r2) / (dm_r1 + dm_r2) jast = (dm_r1 * j_mu_r1 + dm_r2 * j_mu_r2) / (dm_r1 + dm_r2)
else if(murho_type==1)then
! J(r1,r2) = j(0.5 * (mu(r1)+mu(r2)),r12), MU(r1,r2) = 0.5 *(mu(r1)+mu(r2))
call grad_mu_of_r_mean_field(r1,mu_r1, dm_r1, grad_mu_r1, grad_dm_r1)
call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2)
mu_tot = 0.5d0 * (mu_r1 + mu_r2)
jast = j12_mu_input(r1, r2, mu_tot)
else if(murho_type==2)then
! MU(r1,r2) = (rho(1) * mu(r1)+ rho(2) * mu(r2))/(rho(1)+rho(2))
! J(r1,r2) = j(MU(r1,r2),r12)
call grad_mu_of_r_mean_field(r1,mu_r1, dm_r1, grad_mu_r1, grad_dm_r1)
call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2)
double precision :: mu_tmp, dm_tot, dm_tot_inv
dm_tot = dm_r1**a_boys + dm_r2**a_boys ! rho(1)**alpha+rho(2)**alpha
if(dm_tot.lt.1.d-12)then
dm_tot_inv = 1.d+12
else
dm_tot_inv = 1.d0/dm_tot
endif
mu_tmp = dm_r1**a_boys * mu_r1 + dm_r2**a_boys * mu_r2 !rho(1)**alpha * mu(r1)+ rho(2)**alpha * mu(r2)
mu_tot = nu_erf * mu_tmp*dm_tot_inv !
r12 = (r1(1) - r2(1)) * (r1(1) - r2(1))
r12 += (r1(2) - r2(2)) * (r1(2) - r2(2))
r12 += (r1(3) - r2(3)) * (r1(3) - r2(3))
r12 = dsqrt(r12)
jast = j_simple(r12,mu_tot)
endif
end end
subroutine grad_j_sum_mu_of_r(r1,r2,jast,grad_jast) subroutine grad_j_sum_mu_of_r(r1,r2,jast,grad_jast)
implicit none implicit none
include 'constants.include.F'
BEGIN_DOC BEGIN_DOC
END_DOC END_DOC
double precision, intent(in) :: r1(3),r2(3) double precision, intent(in) :: r1(3),r2(3)
@ -26,8 +55,25 @@ subroutine grad_j_sum_mu_of_r(r1,r2,jast,grad_jast)
double precision :: num, denom, grad_num(3), grad_denom(3) double precision :: num, denom, grad_num(3), grad_denom(3)
double precision :: j_r1, grad_j_r1(3),j_r2, grad_j_r2(3) double precision :: j_r1, grad_j_r1(3),j_r2, grad_j_r2(3)
double precision :: dm_r1, grad_dm_r1(3), grad_jmu_r2(3) double precision :: dm_r1, grad_dm_r1(3), grad_jmu_r2(3)
double precision :: dm_r2, grad_dm_r2(3),mu_r2, grad_mu_r2(3) double precision :: dm_r2, grad_dm_r2(3),mu_r2, grad_mu_r2(3),mu_r1
double precision :: j12_mu_input double precision :: j12_mu_input,r12,grad_mu_r1(3),grad_jmu_r1(3)
double precision :: mu_tot,dx,dy,dz,r12_vec(3),d_dmu_j,d_dr12_j
dx = r1(1) - r2(1)
dy = r1(2) - r2(2)
dz = r1(3) - r2(3)
r12 = dsqrt(dx * dx + dy * dy + dz * dz)
if(r12.gt.1.d-10)then
r12_vec(1) = dx
r12_vec(2) = dy
r12_vec(3) = dz
r12_vec *= 1.d0/r12
! r12_vec = grad_r1 (r12)
else
r12 = 1.d-10
r12_vec = 0.d0
endif
if(murho_type==0)then if(murho_type==0)then
! J(r1,r2) = [rho(r1) * j(mu(r1),r12) + rho(r2) * j(mu(r2),r12)] / [rho(r1) + rho(r2)] ! J(r1,r2) = [rho(r1) * j(mu(r1),r12) + rho(r2) * j(mu(r2),r12)] / [rho(r1) + rho(r2)]
@ -48,11 +94,65 @@ subroutine grad_j_sum_mu_of_r(r1,r2,jast,grad_jast)
jast = num / denom jast = num / denom
grad_denom = grad_dm_r1 grad_denom = grad_dm_r1
call grad_j12_mu_input(r1, r2, mu_r2, grad_jmu_r2,r12)
call grad_j12_mu_input(r1, r2, mu_r2, grad_jmu_r2)
grad_num = j_r1 * grad_dm_r1 + dm_r1 * grad_j_r1 + dm_r2 * grad_jmu_r2 grad_num = j_r1 * grad_dm_r1 + dm_r1 * grad_j_r1 + dm_r2 * grad_jmu_r2
grad_jast = (grad_num * denom - num * grad_denom)/(denom*denom) grad_jast = (grad_num * denom - num * grad_denom)/(denom*denom)
else if(murho_type==1)then
! J(r1,r2) = j(0.5 * (mu(r1)+mu(r2)),r12), MU(r1,r2) = 0.5 *(mu(r1)+mu(r2))
!
! d/dx1 J(r1,r2) = d/dx1 j(MU(r1,r2),r12)|MU=cst
! + d/dMU [j(MU,r12)]
! x d/d(mu(r1)) MU(r1,r2)
! x d/dx1 mu(r1)
! = 0.5 * (1 - erf(MU(r1,r2) *r12))/r12 * (x1 - x2) == grad_jmu_r1
! + e^{-(r12*MU(r1,r2))^2}/(2 sqrt(pi) * MU(r1,r2)^2)
! x 0.5
! x d/dx1 mu(r1)
call grad_mu_of_r_mean_field(r1,mu_r1, dm_r1, grad_mu_r1, grad_dm_r1)
call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2)
mu_tot = 0.5d0 * (mu_r1 + mu_r2)
call grad_j12_mu_input(r1, r2, mu_tot, grad_jmu_r1,r12)
grad_jast = grad_jmu_r1
grad_jast+= dexp(-r12*mu_tot*r12*mu_tot) * inv_sq_pi_2 /(mu_tot* mu_tot) * 0.5d0 * grad_mu_r1
else if(murho_type==2)then
! MU(r1,r2) = beta * (rho(1)**alpha * mu(r1)+ rho(2)**alpha * mu(r2))/(rho(1)**alpha+rho(2)**alpha)
! J(r1,r2) = j(MU(r1,r2),r12)
!
! d/dx1 J(r1,r2) = d/dx1 j(MU(r1,r2),r12)|MU=cst
! + d/dMU [j(MU,r12)]
! x d/d(mu(r1)) MU(r1,r2)
! x d/dx1 mu(r1)
! = 0.5 * (1 - erf(MU(r1,r2) *r12))/r12 * (x1 - x2) == grad_jmu_r1
! + 0.5 e^{-(r12*MU(r1,r2))^2}/(2 sqrt(pi) * MU(r1,r2)^2)
! x d/dx1 MU(r1,r2)
! with d/dx1 MU(r1,r2) = beta * {[mu(1) d/dx1 [rho(1)**alpha] + rho(1)**alpha * d/dx1 mu(1)](rho(1)**alpha+rho(2)**alpha)
! - MU(1,2) d/dx1 [rho(1)]**alpha}/(rho(1)**alpha+rho(2)**alpha)^2
! d/dx1 [rho(1)]**alpha = alpha [rho(1)]**(alpha-1) d/dx1 rho(1)
!
call grad_mu_of_r_mean_field(r1,mu_r1, dm_r1, grad_mu_r1, grad_dm_r1)
call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2)
double precision :: dm_tot,dm_tot_inv,grad_mu_tot(3),mu_tmp,grad_dm_r1_alpha(3),d_dx_j
dm_tot = dm_r1**a_boys + dm_r2**a_boys ! rho(1)**alpha+rho(2)**alpha
grad_dm_r1_alpha = a_boys * dm_r1**(a_boys-1) * grad_dm_r1
if(dm_tot.lt.1.d-12)then
dm_tot_inv = 1.d+12
else
dm_tot_inv = 1.d0/dm_tot
endif
mu_tmp = dm_r1**a_boys * mu_r1 + dm_r2**a_boys * mu_r2 !rho(1)**alpha * mu(r1)+ rho(2)**alpha * mu(r2)
mu_tot = nu_erf * mu_tmp*dm_tot_inv !
grad_mu_tot = ( mu_r1 * grad_dm_r1_alpha + dm_r1**a_boys * grad_mu_r1 ) * dm_tot &
- mu_tmp * grad_dm_r1_alpha
grad_mu_tot *= dm_tot_inv * dm_tot_inv * nu_erf
call get_deriv_r12_j12(r12,mu_tot,d_dr12_j) ! d/dr12 j(MU(r1,r2,r12)
! d/dx1 j(MU(r1,r2),r12) | MU(r1,r2) = cst
! d/dr12 j(MU(r1,r2,r12) x d/dx1 r12
grad_jmu_r1 = d_dr12_j * r12_vec
! call grad_j12_mu_input(r1, r2, mu_tot, grad_jmu_r1,r12)
grad_jast = grad_jmu_r1
! d/dMU j(MU(r1,r2),r12)
call get_deriv_mu_j12(r12,mu_tot,d_dmu_j)
grad_jast+= d_dmu_j * grad_mu_tot
else if(murho_type==-1)then else if(murho_type==-1)then
! J(r1,r2) = 0.5 * [j(mu(r1),r12) + j(mu(r2),r12)] ! J(r1,r2) = 0.5 * [j(mu(r1),r12) + j(mu(r2),r12)]
! !
@ -60,7 +160,7 @@ subroutine grad_j_sum_mu_of_r(r1,r2,jast,grad_jast)
call grad_j_mu_of_r_1(r1,r2,j_r1, grad_j_r1,dm_r1, grad_dm_r1) call grad_j_mu_of_r_1(r1,r2,j_r1, grad_j_r1,dm_r1, grad_dm_r1)
call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2) call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2)
j_r2 = j12_mu_input(r1, r2, mu_r2) ! j(mu(r2),r1,r2) j_r2 = j12_mu_input(r1, r2, mu_r2) ! j(mu(r2),r1,r2)
call grad_j12_mu_input(r1, r2, mu_r2, grad_jmu_r2) call grad_j12_mu_input(r1, r2, mu_r2, grad_jmu_r2,r12)
jast = 0.5d0 * (j_r1 + j_r2) jast = 0.5d0 * (j_r1 + j_r2)
grad_jast = 0.5d0 * (grad_j_r1 + grad_jmu_r2) grad_jast = 0.5d0 * (grad_j_r1 + grad_jmu_r2)
@ -73,6 +173,7 @@ subroutine grad_j_mu_of_r_1(r1,r2,jast, grad_jast, dm_r1, grad_dm_r1)
include 'constants.include.F' include 'constants.include.F'
BEGIN_DOC BEGIN_DOC
! grad_r1 of j(mu(r1),r12) ! grad_r1 of j(mu(r1),r12)
!
! !
! d/dx1 j(mu(r1),r12) = exp(-(mu(r1)*r12)**2) /(2 *sqrt(pi) * mu(r1)**2 ) d/dx1 mu(r1) ! d/dx1 j(mu(r1),r12) = exp(-(mu(r1)*r12)**2) /(2 *sqrt(pi) * mu(r1)**2 ) d/dx1 mu(r1)
! + d/dx1 j(mu(r1),r12) ! + d/dx1 j(mu(r1),r12)
@ -83,6 +184,8 @@ subroutine grad_j_mu_of_r_1(r1,r2,jast, grad_jast, dm_r1, grad_dm_r1)
! !
! and d/dx1 j(mu,r12) = 0.5 * (1 - erf(mu *r12))/r12 * (x1 - x2) ! and d/dx1 j(mu,r12) = 0.5 * (1 - erf(mu *r12))/r12 * (x1 - x2)
! !
! d/d mu j(mu,r12) = e^{-(r12*mu)^2}/(2 sqrt(pi) * mu^2)
!
! here mu(r1) is obtained by MU MEAN FIELD ! here mu(r1) is obtained by MU MEAN FIELD
END_DOC END_DOC
double precision, intent(in) :: r1(3),r2(3) double precision, intent(in) :: r1(3),r2(3)
@ -140,14 +243,16 @@ double precision function j12_mu_input(r1, r2, mu)
end end
subroutine grad_j12_mu_input(r1, r2, mu, grad_jmu) subroutine grad_j12_mu_input(r1, r2, mu, grad_jmu,r12)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! grad_jmu = ! grad_jmu = d/dx1 j(mu,r12) assuming mu=cst(r1)
!
! = 0.5/r_12 * (x_1 - x_2) * [1 - erf(mu*r12)]
END_DOC END_DOC
double precision, intent(in) :: r1(3), r2(3), mu double precision, intent(in) :: r1(3), r2(3), mu
double precision, intent(out):: grad_jmu(3) double precision, intent(out):: grad_jmu(3),r12
double precision :: mu_tmp, r12, dx, dy, dz, grad(3), tmp double precision :: mu_tmp, dx, dy, dz, grad(3), tmp
grad_jmu = 0.d0 grad_jmu = 0.d0
dx = r1(1) - r2(1) dx = r1(1) - r2(1)
dy = r1(2) - r2(2) dy = r1(2) - r2(2)
@ -162,3 +267,40 @@ subroutine grad_j12_mu_input(r1, r2, mu, grad_jmu)
grad_jmu = grad grad_jmu = grad
end end
subroutine j12_and_grad_j12_mu_input(r1, r2, mu, jmu, grad_jmu)
implicit none
include 'constants.include.F'
BEGIN_DOC
! jmu = j(mu,r12)
! grad_jmu = d/dx1 j(mu,r12) assuming mu=cst(r1)
!
! = 0.5/r_12 * (x_1 - x_2) * [1 - erf(mu*r12)]
END_DOC
double precision, intent(in) :: r1(3), r2(3), mu
double precision, intent(out):: grad_jmu(3), jmu
double precision :: mu_tmp, r12, dx, dy, dz, grad(3), tmp
double precision :: erfc_mur12,inv_mu
inv_mu = 1.d0/mu
grad_jmu = 0.d0 ! initialization when r12 --> 0
jmu = - inv_sq_pi_2 * inv_mu ! initialization when r12 --> 0
dx = r1(1) - r2(1)
dy = r1(2) - r2(2)
dz = r1(3) - r2(3)
r12 = dsqrt(dx * dx + dy * dy + dz * dz)
if(r12 .lt. 1d-10) return
erfc_mur12 = (1.d0 - derf(mu_tmp))
mu_tmp = mu * r12
tmp = 0.5d0 * erfc_mur12 / r12 ! d/dx1 j(mu(r1),r12)
grad(1) = tmp * dx
grad(2) = tmp * dy
grad(3) = tmp * dz
grad_jmu = grad
jmu= 0.5d0 * r12 * erfc_mur12 - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) * inv_mu
end

View File

@ -1,3 +1,65 @@
subroutine get_deriv_r12_j12(x,mu,d_dx_j)
implicit none
include 'constants.include.F'
BEGIN_DOC
! d/dr12 j(mu,r12)
END_DOC
double precision, intent(in) :: x,mu
double precision, intent(out) :: d_dx_j
d_dx_j = 0.d0
if(x .lt. 1d-10) return
if(j2e_type .eq. "Mu" .or. j2e_type .eq. "Mur") then
d_dx_j = 0.5d0 * (1.d0 - derf(mu * x))
else if(j2e_type .eq. "Mugauss" .or. j2e_type .eq. "Murgauss" ) then
double precision :: x_tmp
x_tmp = mu * x
! gradient of j(mu,x)
d_dx_j = 0.5d0 * (1.d0 - derf(x_tmp))
! gradient of gaussian additional term
x_tmp *= alpha_mu_gauss
x_tmp *= x_tmp
d_dx_j += -0.5d0 * mu * c_mu_gauss * x * dexp(-x_tmp)
else
print *, ' Error in get_deriv_r12_j12: Unknown j2e_type = ', j2e_type
stop
endif
end
subroutine get_deriv_mu_j12(x,mu,d_d_mu)
implicit none
BEGIN_DOC
! d/dmu j(mu,r12)
END_DOC
include 'constants.include.F'
double precision, intent(in) :: x,mu
double precision, intent(out) :: d_d_mu
double precision :: x_tmp,inv_mu_2,inv_alpha_2
d_d_mu = 0.d0
if(x .lt. 1d-10) return
x_tmp = x*mu
if(mu.lt.1.d-10) return
inv_mu_2 = mu*mu
inv_mu_2 = 1.d0/inv_mu_2
if(j2e_type .eq. "Mu" .or. j2e_type .eq. "Mur") then
! e^{-(r12*mu)^2}/(2 sqrt(pi) * mu^2)
d_d_mu = dexp(-x_tmp*x_tmp) * inv_sq_pi_2 * inv_mu_2
else if(j2e_type .eq. "Mugauss" .or. j2e_type .eq. "Murgauss" ) then
d_d_mu = dexp(-x_tmp*x_tmp) * inv_sq_pi_2 * inv_mu_2
inv_alpha_2 = 1.d0/alpha_mu_gauss
inv_alpha_2 *= inv_alpha_2
x_tmp *= alpha_mu_gauss
x_tmp *= x_tmp
d_d_mu += -0.25d0 * c_mu_gauss*inv_alpha_2*dexp(-x_tmp) * (1.d0 + 2.d0 * x_tmp) * inv_mu_2
else
print *, ' Error in get_deriv_r12_j12: Unknown j2e_type = ', j2e_type
stop
endif
end
! --- ! ---
@ -21,6 +83,18 @@ double precision function j12_mu(r1, r2)
j12_mu = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu_erf j12_mu = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu_erf
else if(j2e_type .eq. "Mugauss") then
r12 = 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)) )
double precision :: r12_tmp
r12_tmp = mu_erf * r12
j12_mu = 0.5d0 * r12 * (1.d0 - derf(r12_tmp)) - inv_sq_pi_2 * dexp(-r12_tmp*r12_tmp) / mu_erf
r12_tmp *= alpha_mu_gauss
j12_mu += 0.25d0 * c_mu_gauss / (alpha_mu_gauss*alpha_mu_gauss*mu_erf) * dexp(-r12_tmp*r12_tmp)
else else
print *, ' Error in j12_mu: Unknown j2e_type = ', j2e_type print *, ' Error in j12_mu: Unknown j2e_type = ', j2e_type
@ -60,7 +134,7 @@ subroutine grad1_j12_mu(r1, r2, grad)
grad = 0.d0 grad = 0.d0
if(j2e_type .eq. "Mu") then if(j2e_type .eq. "Mu".or.j2e_type .eq. "Mugauss") then
dx = r1(1) - r2(1) dx = r1(1) - r2(1)
dy = r1(2) - r2(2) dy = r1(2) - r2(2)
@ -69,48 +143,42 @@ subroutine grad1_j12_mu(r1, r2, grad)
r12 = dsqrt(dx * dx + dy * dy + dz * dz) r12 = dsqrt(dx * dx + dy * dy + dz * dz)
if(r12 .lt. 1d-10) return if(r12 .lt. 1d-10) return
tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12 call get_deriv_r12_j12(r12,mu_erf,tmp)
! tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12
grad(1) = tmp * dx grad(1) = tmp * dx
grad(2) = tmp * dy grad(2) = tmp * dy
grad(3) = tmp * dz grad(3) = tmp * dz
grad *= 1.d0/r12
elseif(j2e_type .eq. "Mur") then elseif(j2e_type .eq. "Mur" .or. j2e_type .eq. "Murgauss") then
double precision :: jast double precision :: jast
call grad_j_sum_mu_of_r(r1,r2,jast,grad) call grad_j_sum_mu_of_r(r1,r2,jast,grad)
! dx = r1(1) - r2(1) elseif(j2e_type .eq. "Bump") then
! dy = r1(2) - r2(2) double precision ::grad_jast(3)
! dz = r1(3) - r2(3) call get_grad_j_bump_mu_of_r(r1,r2,grad_jast)
! dx = r1(1) - r2(1)
! r12 = dsqrt(dx * dx + dy * dy + dz * dz) dy = r1(2) - r2(2)
! if(r12 .lt. 1d-10) return dz = r1(3) - r2(3)
!
! tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12 r12 = dsqrt(dx * dx + dy * dy + dz * dz)
! if(r12 .lt. 1d-10) then
! grad(1) = tmp * dx grad(1) = 0.d0
! grad(2) = tmp * dy grad(2) = 0.d0
! grad(3) = tmp * dz grad(3) = 0.d0
return
endif
tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12
grad(1) = 0.5d0 * tmp * dx
grad(2) = 0.5d0 * tmp * dy
grad(3) = 0.5d0 * tmp * dz
grad(1) += 0.5d0 * grad_jast(1)
grad(2) += 0.5d0 * grad_jast(2)
grad(3) += 0.5d0 * grad_jast(3)
! elseif(j2e_type .eq. "Mur") then
!
! dx = r1(1) - r2(1)
! dy = r1(2) - r2(2)
! dz = r1(3) - r2(3)
! r12 = dsqrt(dx * dx + dy * dy + dz * dz)
!
! call mu_r_val_and_grad(r1, r2, mu_val, mu_der)
! mu_tmp = mu_val * r12
! tmp = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val)
! grad(1) = tmp * mu_der(1)
! grad(2) = tmp * mu_der(2)
! grad(3) = tmp * mu_der(3)
!
! if(r12 .lt. 1d-10) return
! tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12
! grad(1) = grad(1) + tmp * dx
! grad(2) = grad(2) + tmp * dy
! grad(3) = grad(3) + tmp * dz
else else

View File

@ -35,6 +35,9 @@ subroutine get_grad1_u12_withsq_r1_seq(ipoint, n_grid2, resx, resy, resz, res)
if( (j2e_type .eq. "Mu") .or. & if( (j2e_type .eq. "Mu") .or. &
(j2e_type .eq. "Mur") .or. & (j2e_type .eq. "Mur") .or. &
(j2e_type .eq. "Mugauss") .or. &
(j2e_type .eq. "Murgauss") .or. &
(j2e_type .eq. "Bump") .or. &
(j2e_type .eq. "Boys") ) then (j2e_type .eq. "Boys") ) then
if(env_type .eq. "None") then if(env_type .eq. "None") then
@ -206,7 +209,43 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz)
gradz(jpoint) = tmp * dz gradz(jpoint) = tmp * dz
enddo enddo
elseif(j2e_type .eq. "Mur") then else if(j2e_type .eq. "Mugauss") then
! d/dx1 j(mu,r12) = 0.5 * [(1 - erf(mu * r12)) / r12 - mu*c*r12*exp(-(mu*alpha*r12)^2] * (x1 - x2)
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)
dx = r1(1) - r2(1)
dy = r1(2) - r2(2)
dz = r1(3) - r2(3)
r12 = dsqrt(dx * dx + dy * dy + dz * dz)
if(r12 .lt. 1d-10) then
gradx(jpoint) = 0.d0
grady(jpoint) = 0.d0
gradz(jpoint) = 0.d0
cycle
endif
double precision :: r12_tmp
r12_tmp = mu_erf * r12
! gradient of j(mu,r12)
tmp = 0.5d0 * (1.d0 - derf(r12_tmp)) / r12
! gradient of gaussian additional term
r12_tmp *= alpha_mu_gauss
r12_tmp *= r12_tmp
tmp += -0.5d0 * mu_erf * c_mu_gauss * r12 * dexp(-r12_tmp)/r12
gradx(jpoint) = tmp * dx
grady(jpoint) = tmp * dy
gradz(jpoint) = tmp * dz
enddo
elseif(j2e_type .eq. "Mur".or.j2e_type .eq. "Murgauss") then
! d/dx1 j(mu(r1,r2),r12) = exp(-(mu(r1,r2)*r12)**2) /(2 *sqrt(pi) * mu(r1,r2)**2 ) d/dx1 mu(r1,r2) ! d/dx1 j(mu(r1,r2),r12) = exp(-(mu(r1,r2)*r12)**2) /(2 *sqrt(pi) * mu(r1,r2)**2 ) d/dx1 mu(r1,r2)
! + 0.5 * (1 - erf(mu(r1,r2) *r12))/r12 * (x1 - x2) ! + 0.5 * (1 - erf(mu(r1,r2) *r12))/r12 * (x1 - x2)
@ -221,30 +260,41 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz)
gradx(jpoint) = grad_jast(1) gradx(jpoint) = grad_jast(1)
grady(jpoint) = grad_jast(2) grady(jpoint) = grad_jast(2)
gradz(jpoint) = grad_jast(3) gradz(jpoint) = grad_jast(3)
! dx = r1(1) - r2(1) enddo
! dy = r1(2) - r2(2) elseif(j2e_type .eq. "Bump") then
! dz = r1(3) - r2(3)
! r12 = dsqrt(dx * dx + dy * dy + dz * dz) ! d/dx1 jbump(r1,r2)
!
! call mu_r_val_and_grad(r1, r2, mu_val, mu_der) do jpoint = 1, n_points_extra_final_grid ! r2
! mu_tmp = mu_val * r12
! tmp = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val) r2(1) = final_grid_points_extra(1,jpoint)
! gradx(jpoint) = tmp * mu_der(1) r2(2) = final_grid_points_extra(2,jpoint)
! grady(jpoint) = tmp * mu_der(2) r2(3) = final_grid_points_extra(3,jpoint)
! gradz(jpoint) = tmp * mu_der(3) call get_grad_j_bump_mu_of_r(r1,r2,grad_jast)
!
! if(r12 .lt. 1d-10) then dx = r1(1) - r2(1)
! gradx(jpoint) = 0.d0 dy = r1(2) - r2(2)
! grady(jpoint) = 0.d0 dz = r1(3) - r2(3)
! gradz(jpoint) = 0.d0
! cycle r12 = dsqrt(dx * dx + dy * dy + dz * dz)
! endif if(r12 .lt. 1d-10) then
! gradx(jpoint) = 0.d0
! tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12 grady(jpoint) = 0.d0
! gradz(jpoint) = 0.d0
! gradx(jpoint) = gradx(jpoint) + tmp * dx cycle
! grady(jpoint) = grady(jpoint) + tmp * dy endif
! gradz(jpoint) = gradz(jpoint) + tmp * dz
tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12
gradx(jpoint) = 0.5d0 * tmp * dx
grady(jpoint) = 0.5d0 * tmp * dy
gradz(jpoint) = 0.5d0 * tmp * dz
gradx(jpoint) += 0.5d0 * grad_jast(1)
grady(jpoint) += 0.5d0 * grad_jast(2)
gradz(jpoint) += 0.5d0 * grad_jast(3)
! gradx(jpoint) = grad_jast(1)
! grady(jpoint) = grad_jast(2)
! gradz(jpoint) = grad_jast(3)
enddo enddo
elseif(j2e_type .eq. "Boys") then elseif(j2e_type .eq. "Boys") then
@ -674,7 +724,10 @@ subroutine get_grad1_u12_2e_r1_seq(ipoint, n_grid2, resx, resy, resz)
r1(3) = final_grid_points(3,ipoint) r1(3) = final_grid_points(3,ipoint)
if( (j2e_type .eq. "Mu") .or. & if( (j2e_type .eq. "Mu") .or. &
(j2e_type .eq. "Mugauss") .or. &
(j2e_type .eq. "Mur") .or. & (j2e_type .eq. "Mur") .or. &
(j2e_type .eq. "Murgauss") .or. &
(j2e_type .eq. "Bump") .or. &
(j2e_type .eq. "Boys") ) then (j2e_type .eq. "Boys") ) then
if(env_type .eq. "None") then if(env_type .eq. "None") then
@ -793,6 +846,9 @@ subroutine get_u12_2e_r1_seq(ipoint, n_grid2, res)
if( (j2e_type .eq. "Mu") .or. & if( (j2e_type .eq. "Mu") .or. &
(j2e_type .eq. "Mur") .or. & (j2e_type .eq. "Mur") .or. &
(j2e_type .eq. "Mugauss") .or. &
(j2e_type .eq. "Murgauss") .or. &
(j2e_type .eq. "Mugauss") .or. &
(j2e_type .eq. "Boys") ) then (j2e_type .eq. "Boys") ) then
if(env_type .eq. "None") then if(env_type .eq. "None") then

View File

@ -0,0 +1,59 @@
program plot_j_gauss
implicit none
double precision :: xmin, xmax, x, dx
double precision :: mu_min, mu_max, mu, d_mu
double precision :: pot_j_gauss,j_mu_simple,j_gauss_simple,pot_j_mu
double precision, allocatable :: mu_tab(:),j_mu(:),j_mu_gauss(:)
double precision, allocatable :: w_mu(:), w_mu_gauss(:)
character*(128) :: output
integer :: getUnitAndOpen
integer :: i_unit_output_wee_gauss,i_unit_output_wee_mu
integer :: i_unit_output_j_gauss,i_unit_output_j_mu
output=trim(ezfio_filename)//'.w_ee_mu_gauss'
i_unit_output_wee_gauss = getUnitAndOpen(output,'w')
output=trim(ezfio_filename)//'.w_ee_mu'
i_unit_output_wee_mu = getUnitAndOpen(output,'w')
output=trim(ezfio_filename)//'.j_mu_gauss'
i_unit_output_j_gauss = getUnitAndOpen(output,'w')
output=trim(ezfio_filename)//'.j_mu'
i_unit_output_j_mu = getUnitAndOpen(output,'w')
integer :: npt, i, j, n_mu
n_mu = 3
allocate(mu_tab(n_mu),j_mu(n_mu),j_mu_gauss(n_mu),w_mu(n_mu), w_mu_gauss(n_mu))
mu_min = 0.5d0
mu_max = 2.d0
d_mu = (mu_max - mu_min)/dble(n_mu)
mu = mu_min
do i = 1, n_mu
mu_tab(i) = mu
print*,'mu = ',mu
mu += d_mu
enddo
mu_tab(1) = 0.9d0
mu_tab(2) = 0.95d0
mu_tab(3) = 1.d0
xmin = 0.01d0
xmax = 10.d0
npt = 1000
dx = (xmax - xmin)/dble(npt)
x = xmin
do i = 1, npt
do j = 1, n_mu
mu = mu_tab(j)
w_mu_gauss(j) = pot_j_gauss(x,mu)
w_mu(j) = pot_j_mu(x,mu)
j_mu(j) = j_mu_simple(x,mu)
j_mu_gauss(j) = j_gauss_simple(x,mu) + j_mu(j)
enddo
write(i_unit_output_wee_gauss,'(100(F16.10,X))')x,w_mu_gauss(:)
write(i_unit_output_wee_mu,'(100(F16.10,X))')x,w_mu(:)
write(i_unit_output_j_gauss,'(100(F16.10,X))')x,j_mu_gauss(:)
write(i_unit_output_j_mu,'(100(F16.10,X))')x,j_mu(:)
x += dx
enddo
end

View File

@ -16,15 +16,16 @@ subroutine routine_print
integer :: ipoint,nx,i integer :: ipoint,nx,i
double precision :: xmax,xmin,r(3),dx,sigma double precision :: xmax,xmin,r(3),dx,sigma
double precision :: mu_val, mu_der(3),dm_a,dm_b,grad,grad_dm_a(3), grad_dm_b(3) double precision :: mu_val, mu_der(3),dm_a,dm_b,grad,grad_dm_a(3), grad_dm_b(3)
xmax = 5.D0 xmax = 3.9D0
xmin = -5.D0 xmin = -3.9D0
nx = 10000 nx = 10000
dx = (xmax - xmin)/dble(nx) dx = (xmax - xmin)/dble(nx)
r = 0.d0 r = 0.d0
r(1) = xmin r(1) = xmin
do ipoint = 1, nx do ipoint = 1, nx
call mu_r_val_and_grad(r, r, mu_val, mu_der) ! call mu_r_val_and_grad(r, r, mu_val, mu_der)
call density_and_grad_alpha_beta(r,dm_a,dm_b, grad_dm_a, grad_dm_b) call grad_mu_of_r_mean_field(r,mu_val, dm_a, mu_der, grad_dm_a)
! call density_and_grad_alpha_beta(r,dm_a,dm_b, grad_dm_a, grad_dm_b)
sigma = 0.d0 sigma = 0.d0
do i = 1,3 do i = 1,3
sigma += grad_dm_a(i)**2 sigma += grad_dm_a(i)**2
@ -32,7 +33,8 @@ subroutine routine_print
sigma=dsqrt(sigma) sigma=dsqrt(sigma)
grad = mu_der(1)**2 + mu_der(2)**2 + mu_der(3)**2 grad = mu_der(1)**2 + mu_der(2)**2 + mu_der(3)**2
grad = dsqrt(grad) grad = dsqrt(grad)
write(i_unit_output,'(100(F16.7,X))')r(1),mu_val,dm_a+dm_b,grad,sigma/dm_a print*,r(1),mu_val
write(i_unit_output,'(100(F16.7,X))')r(1),mu_val,dm_a,grad,sigma/dm_a
r(1) += dx r(1) += dx
enddo enddo
end end

View File

@ -0,0 +1,146 @@
double precision function j_simple(x,mu)
implicit none
double precision, intent(in) :: x,mu
double precision :: j_mu_simple,j_gauss_simple
if(j2e_type .eq. "Mu".or.j2e_type .eq. "Mur") then
j_simple = j_mu_simple(x,mu)
else if(j2e_type .eq. "Mugauss".or.j2e_type .eq. "Murgauss") then
j_simple = j_gauss_simple(x,mu) + j_mu_simple(x,mu)
endif
end
double precision function j_mu_simple(x,mu)
implicit none
double precision, intent(in):: x,mu
include 'constants.include.F'
BEGIN_DOC
! j_mu(mu,x) = 0.5 x (1 - erf(mu x)) - 1/[2 sqrt(pi)mu] exp(-(x*mu)^2)
END_DOC
j_mu_simple = 0.5d0 * x * (1.D0 - derf(mu*x)) - 0.5d0 * inv_sq_pi/mu * dexp(-x*mu*x*mu)
end
double precision function j_gauss_simple(x,mu)
implicit none
double precision, intent(in):: x,mu
include 'constants.include.F'
BEGIN_DOC
! j_mu(mu,x) = c/[4 alpha^2 mu] exp(-(alpha * mu * x)^2)
! with c = 27/(8 sqrt(pi)), alpha=3/2
END_DOC
double precision :: x_tmp
x_tmp = alpha_mu_gauss * mu * x
j_gauss_simple = 0.25d0 * c_mu_gauss / (alpha_mu_gauss*alpha_mu_gauss*mu) * dexp(-x_tmp*x_tmp)
end
double precision function j_mu_deriv(x,mu)
implicit none
BEGIN_DOC
! d/dx j_mu(mu,x) = d/dx 0.5 x (1 - erf(mu x)) - 1/[2 sqrt(pi)mu] exp(-(x*mu)^2)
! = 0.5*(1 - erf(mu x))
END_DOC
include 'constants.include.F'
double precision, intent(in) :: x,mu
j_mu_deriv = 0.5d0 * (1.d0 - derf(mu*x))
end
double precision function j_mu_deriv_2(x,mu)
implicit none
BEGIN_DOC
! d^2/dx^2 j_mu(mu,x) = d^2/dx^2 0.5 x (1 - erf(mu x)) - 1/[2 sqrt(pi)mu] exp(-(x*mu)^2)
! = -mu/sqrt(pi) * exp(-(mu x)^2)
END_DOC
include 'constants.include.F'
double precision, intent(in) :: x,mu
j_mu_deriv_2 = - mu * inv_sq_pi * dexp(-x*mu*x*mu)
end
double precision function j_gauss_deriv(x,mu)
implicit none
include 'constants.include.F'
double precision, intent(in) :: x,mu
BEGIN_DOC
! d/dx j_gauss(mu,x) = d/dx c/[4 alpha^2 mu] exp(-(alpha * mu * x)^2)
! with c = 27/(8 sqrt(pi)), alpha=3/2
! = -0.5 * mu * c * x * exp(-(alpha * mu * x)^2)
END_DOC
double precision :: x_tmp
x_tmp = alpha_mu_gauss * mu * x
j_gauss_deriv = -0.5d0 * mu * c_mu_gauss * x * exp(-x_tmp*x_tmp)
end
double precision function j_gauss_deriv_2(x,mu)
implicit none
include 'constants.include.F'
double precision, intent(in) :: x,mu
BEGIN_DOC
! d/dx j_gauss(mu,x) = d/dx c/[4 alpha^2 mu] exp(-(alpha * mu * x)^2)
! with c = 27/(8 sqrt(pi)), alpha=3/2
! = 0.5 * mu * c * exp(-(alpha * mu * x)^2) * (2 (alpha*mu*x)^2 - 1)
END_DOC
double precision :: x_tmp
x_tmp = alpha_mu_gauss * mu * x
x_tmp = x_tmp * x_tmp
j_gauss_deriv_2 = 0.5d0 * mu * c_mu_gauss * exp(-x_tmp) * (2.d0*x_tmp - 1.d0)
end
double precision function j_erf_gauss_deriv(x,mu)
implicit none
double precision, intent(in) :: x,mu
BEGIN_DOC
! d/dx (j_gauss(mu,x)+j_mu(mu,x))
END_DOC
double precision :: j_gauss_deriv,j_mu_deriv
j_erf_gauss_deriv = j_gauss_deriv(x,mu)+j_mu_deriv(x,mu)
end
double precision function j_erf_gauss_deriv_2(x,mu)
implicit none
double precision, intent(in) :: x,mu
BEGIN_DOC
! d^2/dx^2 (j_gauss(mu,x)+j_mu(mu,x))
END_DOC
double precision :: j_gauss_deriv_2,j_mu_deriv_2
j_erf_gauss_deriv_2 = j_gauss_deriv_2(x,mu)+j_mu_deriv_2(x,mu)
end
double precision function pot_j_gauss(x,mu)
implicit none
double precision, intent(in) :: x,mu
BEGIN_DOC
! effective scalar potential associated with the erf_gauss correlation factor
!
! 1/x( 1 - 2 * d/dx j_erf_gauss(x,mu)) - d^2/dx^2 j_erf_gauss(x,mu)) - d/dx d/dx (j_erf_gauss(x,mu))^2
END_DOC
double precision :: j_erf_gauss_deriv_2,j_erf_gauss_deriv
double precision :: deriv_1, deriv_2
pot_j_gauss = 0.d0
if(x.ne.0.d0)then
deriv_1 = j_erf_gauss_deriv(x,mu)
deriv_2 = j_erf_gauss_deriv_2(x,mu)
pot_j_gauss = 1.d0/x * (1.d0 - 2.d0 * deriv_1) - deriv_1 * deriv_1 - deriv_2
endif
end
double precision function pot_j_mu(x,mu)
implicit none
double precision, intent(in) :: x,mu
BEGIN_DOC
! effective scalar potential associated with the correlation factor
!
! 1/x( 1 - 2 * d/dx j_erf(x,mu)) - d^2/dx^2 j_erf(x,mu)) - d/dx d/dx (j_erf(x,mu))^2
END_DOC
double precision :: j_mu_deriv_2,j_mu_deriv
double precision :: deriv_1, deriv_2
pot_j_mu = 0.d0
if(x.ne.0.d0)then
deriv_1 = j_mu_deriv(x,mu)
deriv_2 = j_mu_deriv_2(x,mu)
pot_j_mu= 1.d0/x * (1.d0 - 2.d0 * deriv_1) - deriv_1 * deriv_1 - deriv_2
endif
end

View File

@ -57,6 +57,9 @@ end
subroutine get_grad_f_mf_ab(r,grad_f_mf_ab, grad_two_bod_dens,f_mf_ab,two_bod_dens, dm_a, dm_b,grad_dm_a, grad_dm_b) subroutine get_grad_f_mf_ab(r,grad_f_mf_ab, grad_two_bod_dens,f_mf_ab,two_bod_dens, dm_a, dm_b,grad_dm_a, grad_dm_b)
implicit none implicit none
BEGIN_DOC
! gradient of mu(r) mean field, together with the gradient of the one- and two-body HF density.
END_DOC
double precision, intent(in) :: r(3) double precision, intent(in) :: r(3)
double precision, intent(out) :: f_mf_ab, two_bod_dens double precision, intent(out) :: f_mf_ab, two_bod_dens
double precision, intent(out) :: grad_two_bod_dens(3), grad_f_mf_ab(3) double precision, intent(out) :: grad_two_bod_dens(3), grad_f_mf_ab(3)

View File

@ -19,9 +19,9 @@ program projected_operators
! call test_f_ii_ia_aa_valence_ab ! call test_f_ii_ia_aa_valence_ab
! call test ! call test
! call test_f_mean_field ! call test_f_mean_field
! call test_grad_f_mean_field call test_grad_f_mean_field
! call test_grad_mu_mf call test_grad_mu_mf
call plot_mu_of_r_mf ! call plot_mu_of_r_mf
end end

View File

@ -9,6 +9,9 @@ double precision, parameter :: pi_5_2 = 34.9868366552d0
double precision, parameter :: dfour_pi = 4.d0*dacos(-1.d0) double precision, parameter :: dfour_pi = 4.d0*dacos(-1.d0)
double precision, parameter :: dtwo_pi = 2.d0*dacos(-1.d0) double precision, parameter :: dtwo_pi = 2.d0*dacos(-1.d0)
double precision, parameter :: inv_sq_pi = 1.d0/dsqrt(dacos(-1.d0)) double precision, parameter :: inv_sq_pi = 1.d0/dsqrt(dacos(-1.d0))
double precision, parameter :: c_mu_gauss = 27.d0/(8.d0*dsqrt(dacos(-1.d0)))
double precision, parameter :: c_mu_gauss_tot = 1.5d0*27.d0/(8.d0*dsqrt(dacos(-1.d0)))+3.d0/dsqrt(dacos(-1.d0))
double precision, parameter :: alpha_mu_gauss = 1.5d0
double precision, parameter :: inv_sq_pi_2 = 0.5d0/dsqrt(dacos(-1.d0)) double precision, parameter :: inv_sq_pi_2 = 0.5d0/dsqrt(dacos(-1.d0))
double precision, parameter :: thresh = 1.d-15 double precision, parameter :: thresh = 1.d-15
double precision, parameter :: cx_lda = -0.73855876638202234d0 double precision, parameter :: cx_lda = -0.73855876638202234d0