mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-12-22 12:23:43 +01:00
new attempt of mu(r)
This commit is contained in:
eginer
parent
cb8bef2ecd
commit
65bef0b266
@ -1,3 +1,14 @@
|
||||
[mu_of_r_tc]
|
||||
type: logical
|
||||
doc: If |true|, take the second formula for mu(r)
|
||||
interface: ezfio,provider,ocaml
|
||||
default: False
|
||||
|
||||
[mu_of_r_av]
|
||||
type: logical
|
||||
doc: If |true|, take the second formula for mu(r)
|
||||
interface: ezfio,provider,ocaml
|
||||
default: False
|
||||
|
||||
[j2e_type]
|
||||
type: character*(32)
|
||||
|
||||
@ -4,3 +4,4 @@ jastrow
|
||||
ao_tc_eff_map
|
||||
bi_ortho_mos
|
||||
trexio
|
||||
mu_of_r
|
||||
|
||||
BIN
plugins/local/non_h_ints_mu/deb_j_mu_of_r
Executable file
BIN
plugins/local/non_h_ints_mu/deb_j_mu_of_r
Executable file
Binary file not shown.
97
plugins/local/non_h_ints_mu/deb_j_mu_of_r.irp.f
Normal file
97
plugins/local/non_h_ints_mu/deb_j_mu_of_r.irp.f
Normal file
@ -0,0 +1,97 @@
|
||||
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)
|
||||
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)
|
||||
double precision :: norm
|
||||
norm = 0.D0
|
||||
do k = 1, 3
|
||||
r1bis= r1
|
||||
r1bis(k) += dr
|
||||
call get_j_sum_mu_of_r(r1bis,r2,jast_p)
|
||||
|
||||
r1bis= r1
|
||||
r1bis(k) -= dr
|
||||
call get_j_sum_mu_of_r(r1bis,r2,jast_m)
|
||||
|
||||
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
|
||||
norm = dsqrt(norm)
|
||||
if(norm.gt.1.d-10)then
|
||||
print*,'/////'
|
||||
print*,grad_jast_mu_r1
|
||||
print*,num_grad_jast_mu_r1
|
||||
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 :: r2(3), weight, dr, r1(3), r1bis(3)
|
||||
double precision :: accu_grad(3)
|
||||
double precision :: jast_mu_r1,grad_jast_mu_r1(3)
|
||||
double precision :: jast_mu_r1_p,jast_mu_r1_m,num_grad_jast_mu_r1(3)
|
||||
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
|
||||
r2 = 0.d0
|
||||
r2(1) = 2.5d0
|
||||
r2(2) = 2.25d0
|
||||
r2(3) = -2.5d0
|
||||
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_mu_of_r_1(r1,r2,jast_mu_r1,grad_jast_mu_r1, dm_r1, grad_dm_r1)
|
||||
do k = 1, 3
|
||||
r1bis= r1
|
||||
r1bis(k) += dr
|
||||
call grad_mu_of_r_mean_field(r1bis,mu_val_p, dm_r1, mu_der, grad_dm_r1)
|
||||
jast_mu_r1_p = j12_mu_input(r1bis, r2, mu_val_p)
|
||||
|
||||
r1bis= r1
|
||||
r1bis(k) -= dr
|
||||
call grad_mu_of_r_mean_field(r1bis,mu_val_m, dm_r1, mu_der, grad_dm_r1)
|
||||
jast_mu_r1_m = j12_mu_input(r1bis, r2, mu_val_m)
|
||||
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
|
||||
enddo
|
||||
print*,'/////'
|
||||
print*,grad_jast_mu_r1
|
||||
print*,num_grad_jast_mu_r1
|
||||
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
|
||||
164
plugins/local/non_h_ints_mu/jast_deriv_mu_of_r.irp.f
Normal file
164
plugins/local/non_h_ints_mu/jast_deriv_mu_of_r.irp.f
Normal file
@ -0,0 +1,164 @@
|
||||
subroutine get_j_sum_mu_of_r(r1,r2,jast)
|
||||
implicit none
|
||||
double precision, intent(in) :: r1(3),r2(3)
|
||||
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_r2, dm_r2, grad_mu_r2(3), grad_dm_r2(3), j_mu_r2
|
||||
double precision :: j12_mu_input
|
||||
jast = 0.d0
|
||||
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)
|
||||
j_mu_r1 = j12_mu_input(r1, r2, mu_r1)
|
||||
j_mu_r2 = j12_mu_input(r1, r2, mu_r2)
|
||||
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)
|
||||
|
||||
end
|
||||
|
||||
subroutine grad_j_sum_mu_of_r(r1,r2,jast,grad_jast)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
END_DOC
|
||||
double precision, intent(in) :: r1(3),r2(3)
|
||||
double precision, intent(out):: jast, grad_jast(3)
|
||||
jast = 0.d0
|
||||
grad_jast = 0.d0
|
||||
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 :: 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 :: j12_mu_input
|
||||
|
||||
if(murho_type==0)then
|
||||
! J(r1,r2) = [rho(r1) * j(mu(r1),r12) + rho(r2) * j(mu(r2),r12)] / [rho(r1) + rho(r2)]
|
||||
!
|
||||
! = num(r1,r2) / denom(r1,r2)
|
||||
!
|
||||
! d/dx1 J(r1,r2) = [denom(r1,r2) X d/dx1 num(r1,r2) - num(r1,r2) X d/dx1 denom(r1,r2) ] / denom(r1,r2)^2
|
||||
!
|
||||
! d/dx1 num(r1,r2) = j(mu(r1),r12)*d/dx1 rho(r1) + rho(r1) * d/dx1 j(mu(r1),r12)
|
||||
! + rho(r2) d/dx1 j(mu(r2),r12)
|
||||
! d/dx1 denom(r1,r2) = d/dx1 rho(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)
|
||||
j_r2 = j12_mu_input(r1, r2, mu_r2) ! j(mu(r2),r1,r2)
|
||||
num = dm_r1 * j_r1 + dm_r2 * j_r2
|
||||
denom = dm_r1 + dm_r2
|
||||
if(denom.lt.1.d-7)return
|
||||
jast = num / denom
|
||||
|
||||
grad_denom = grad_dm_r1
|
||||
|
||||
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_jast = (grad_num * denom - num * grad_denom)/(denom*denom)
|
||||
else if(murho_type==-1)then
|
||||
! J(r1,r2) = 0.5 * [j(mu(r1),r12) + j(mu(r2),r12)]
|
||||
!
|
||||
! d/dx1 J(r1,r2) = 0.5 * (d/dx1 j(mu(r1),r12) + d/dx1 j(mu(r2),r12))
|
||||
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)
|
||||
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)
|
||||
jast = 0.5d0 * (j_r1 + j_r2)
|
||||
grad_jast = 0.5d0 * (grad_j_r1 + grad_jmu_r2)
|
||||
|
||||
endif
|
||||
|
||||
end
|
||||
|
||||
subroutine grad_j_mu_of_r_1(r1,r2,jast, grad_jast, dm_r1, grad_dm_r1)
|
||||
implicit none
|
||||
include 'constants.include.F'
|
||||
BEGIN_DOC
|
||||
! 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)
|
||||
!
|
||||
! with
|
||||
!
|
||||
! j(mu,r12) = 1/2 r12 (1 - erf(mu r12)) - 1/2 (sqrt(pi) * mu) e^{-(mu*r12)^2}
|
||||
!
|
||||
! and d/dx1 j(mu,r12) = 0.5 * (1 - erf(mu *r12))/r12 * (x1 - x2)
|
||||
!
|
||||
! here mu(r1) is obtained by MU MEAN FIELD
|
||||
END_DOC
|
||||
double precision, intent(in) :: r1(3),r2(3)
|
||||
double precision, intent(out):: jast, grad_jast(3),dm_r1, grad_dm_r1(3)
|
||||
double precision :: dx, dy, dz, r12, mu_der(3)
|
||||
double precision :: mu_tmp, tmp, grad(3), mu_val
|
||||
jast = 0.d0
|
||||
grad = 0.d0
|
||||
|
||||
dx = r1(1) - r2(1)
|
||||
dy = r1(2) - r2(2)
|
||||
dz = r1(3) - r2(3)
|
||||
r12 = dsqrt(dx * dx + dy * dy + dz * dz)
|
||||
! get mu(r1) == mu_val and its gradient d/dx1 mu(r1) == mu_der
|
||||
call grad_mu_of_r_mean_field(r1,mu_val, dm_r1, mu_der, grad_dm_r1)
|
||||
mu_tmp = mu_val * r12
|
||||
! evalulation of the jastrow j(mu(r1),r12)
|
||||
jast = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu_val
|
||||
|
||||
! tmp = exp(-(mu(r1)*r12)**2) /(2 *sqrt(pi) * mu(r1)**2 )
|
||||
tmp = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val)
|
||||
! grad =
|
||||
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 ! d/dx1 j(mu(r1),r12)
|
||||
grad(1) = grad(1) + tmp * dx
|
||||
grad(2) = grad(2) + tmp * dy
|
||||
grad(3) = grad(3) + tmp * dz
|
||||
|
||||
grad_jast = grad
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
double precision function j12_mu_input(r1, r2, mu)
|
||||
|
||||
BEGIN_DOC
|
||||
! j(mu,r12) = 1/2 r12 (1 - erf(mu r12)) - 1/2 (sqrt(pi) * mu) e^{-(mu*r12)^2}
|
||||
END_DOC
|
||||
include 'constants.include.F'
|
||||
|
||||
implicit none
|
||||
double precision, intent(in) :: r1(3), r2(3), mu
|
||||
double precision :: mu_tmp, r12
|
||||
|
||||
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)) )
|
||||
mu_tmp = mu * r12
|
||||
|
||||
j12_mu_input = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu
|
||||
|
||||
end
|
||||
|
||||
subroutine grad_j12_mu_input(r1, r2, mu, grad_jmu)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! grad_jmu =
|
||||
END_DOC
|
||||
double precision, intent(in) :: r1(3), r2(3), mu
|
||||
double precision, intent(out):: grad_jmu(3)
|
||||
double precision :: mu_tmp, r12, dx, dy, dz, grad(3), tmp
|
||||
grad_jmu = 0.d0
|
||||
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
|
||||
mu_tmp = mu * r12
|
||||
tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12 ! d/dx1 j(mu(r1),r12)
|
||||
grad(1) = tmp * dx
|
||||
grad(2) = tmp * dy
|
||||
grad(3) = tmp * dz
|
||||
|
||||
grad_jmu = grad
|
||||
end
|
||||
@ -3,6 +3,9 @@
|
||||
|
||||
double precision function j12_mu(r1, r2)
|
||||
|
||||
BEGIN_DOC
|
||||
! j(mu,r12) = 1/2 r12 (1 - erf(mu r12)) - 1/2 (sqrt(pi) * mu) e^{-(mu*r12)^2}
|
||||
END_DOC
|
||||
include 'constants.include.F'
|
||||
|
||||
implicit none
|
||||
@ -73,24 +76,41 @@ subroutine grad1_j12_mu(r1, r2, grad)
|
||||
grad(3) = tmp * dz
|
||||
|
||||
elseif(j2e_type .eq. "Mur") then
|
||||
double precision :: jast
|
||||
call grad_j_sum_mu_of_r(r1,r2,jast,grad)
|
||||
|
||||
dx = r1(1) - r2(1)
|
||||
dy = r1(2) - r2(2)
|
||||
dz = r1(3) - r2(3)
|
||||
r12 = dsqrt(dx * dx + dy * dy + dz * dz)
|
||||
! 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
|
||||
!
|
||||
! tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12
|
||||
!
|
||||
! grad(1) = tmp * dx
|
||||
! grad(2) = tmp * dy
|
||||
! grad(3) = tmp * 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
|
||||
! 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
|
||||
|
||||
@ -369,7 +389,18 @@ end
|
||||
! ---
|
||||
|
||||
subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der)
|
||||
|
||||
BEGIN_DOC
|
||||
! various flavours of mu(r1,r2)
|
||||
! depends on essentially the density and other related quantities
|
||||
!
|
||||
! change the variable "murho_type" to change type
|
||||
!
|
||||
! murho_type == -1 :: mu(r1,r2) = (rho(r1) mu_mf(r1) + rho(r2) mu_mf(r2))/[rho(r1)+rho(r2)]
|
||||
!
|
||||
! == 0 :: mu(r1,r2) = (sqrt(rho(r1)) mu_mf(r1) + sqrt(rho(r2)) mu_mf(r2))/[sqrt(rho(r1))+sqrt(rho(r2))]
|
||||
!
|
||||
! == -2 :: mu(r1,r2) = 0.5(mu_mf(r1) + mu_mf(r2))
|
||||
END_DOC
|
||||
implicit none
|
||||
double precision, intent(in) :: r1(3), r2(3)
|
||||
double precision, intent(out) :: mu_val, mu_der(3)
|
||||
@ -379,11 +410,50 @@ subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der)
|
||||
double precision :: rho1, grad_rho1(3),rho2,rho_tot,inv_rho_tot
|
||||
double precision :: f_rho1, f_rho2, d_drho_f_rho1
|
||||
double precision :: d_dx1_f_rho1(3),d_dx_rho_f_rho(3),nume
|
||||
double precision :: mu_mf_r1, dm_r1, grad_mu_mf_r1(3), grad_dm_r1(3)
|
||||
double precision :: mu_mf_r2, dm_r2, grad_mu_mf_r2(3), grad_dm_r2(3)
|
||||
|
||||
double precision :: num, denom, grad_denom(3), grad_num(3)
|
||||
double precision :: dsqrt_dm_r1
|
||||
|
||||
PROVIDE murho_type
|
||||
PROVIDE mu_r_ct mu_erf
|
||||
|
||||
if(murho_type .eq. 1) then
|
||||
if(murho_type .eq. 0) then
|
||||
call grad_mu_of_r_mean_field(r1,mu_mf_r1, dm_r1, grad_mu_mf_r1, grad_dm_r1)
|
||||
call grad_mu_of_r_mean_field(r2,mu_mf_r2, dm_r2, grad_mu_mf_r2, grad_dm_r2)
|
||||
dsqrt_dm_r1 = dsqrt(dm_r1)
|
||||
denom = (dsqrt_dm_r1 + dsqrt(dm_r2) )
|
||||
if(denom.lt.1.d-7)then
|
||||
mu_val = 1.d+10
|
||||
mu_der = 0.d0
|
||||
return
|
||||
endif
|
||||
num = (dsqrt(dm_r1) * mu_mf_r1 + dsqrt(dm_r2) * mu_mf_r2)
|
||||
mu_val = num / denom
|
||||
grad_denom = grad_dm_r1/dsqrt_dm_r1
|
||||
grad_num = dsqrt(dm_r1) * grad_mu_mf_r1 + mu_mf_r1 * grad_dm_r1
|
||||
mu_der = (grad_num * denom - num * grad_denom)/(denom*denom)
|
||||
else if(murho_type .eq. -1) then
|
||||
call grad_mu_of_r_mean_field(r1,mu_mf_r1, dm_r1, grad_mu_mf_r1, grad_dm_r1)
|
||||
call grad_mu_of_r_mean_field(r2,mu_mf_r2, dm_r2, grad_mu_mf_r2, grad_dm_r2)
|
||||
denom = (dm_r1 + dm_r2 )
|
||||
if(denom.lt.1.d-7)then
|
||||
mu_val = 1.d+10
|
||||
mu_der = 0.d0
|
||||
return
|
||||
endif
|
||||
num = (dm_r1 * mu_mf_r1 + dm_r2 * mu_mf_r2)
|
||||
mu_val = num / denom
|
||||
grad_denom = grad_dm_r1
|
||||
grad_num = dm_r1 * grad_mu_mf_r1 + mu_mf_r1 * grad_dm_r1
|
||||
mu_der = (grad_num * denom - num * grad_denom)/(denom*denom)
|
||||
else if(murho_type .eq. -2) then
|
||||
call grad_mu_of_r_mean_field(r1,mu_mf_r1, dm_r1, grad_mu_mf_r1, grad_dm_r1)
|
||||
call grad_mu_of_r_mean_field(r2,mu_mf_r2, dm_r2, grad_mu_mf_r2, grad_dm_r2)
|
||||
mu_val = 0.5d0 * (mu_mf_r1 + mu_mf_r2)
|
||||
mu_der = 0.5d0 * grad_mu_mf_r1
|
||||
else if(murho_type .eq. 1) then
|
||||
|
||||
!
|
||||
! r = 0.5 (r1 + r2)
|
||||
|
||||
@ -216,31 +216,35 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz)
|
||||
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)
|
||||
|
||||
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)
|
||||
gradx(jpoint) = tmp * mu_der(1)
|
||||
grady(jpoint) = tmp * mu_der(2)
|
||||
gradz(jpoint) = tmp * mu_der(3)
|
||||
|
||||
if(r12 .lt. 1d-10) then
|
||||
gradx(jpoint) = 0.d0
|
||||
grady(jpoint) = 0.d0
|
||||
gradz(jpoint) = 0.d0
|
||||
cycle
|
||||
endif
|
||||
|
||||
tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12
|
||||
|
||||
gradx(jpoint) = gradx(jpoint) + tmp * dx
|
||||
grady(jpoint) = grady(jpoint) + tmp * dy
|
||||
gradz(jpoint) = gradz(jpoint) + tmp * dz
|
||||
double precision :: jast, grad_jast(3)
|
||||
call grad_j_sum_mu_of_r(r1,r2,jast,grad_jast)
|
||||
gradx(jpoint) = grad_jast(1)
|
||||
grady(jpoint) = grad_jast(2)
|
||||
gradz(jpoint) = grad_jast(3)
|
||||
! 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)
|
||||
! gradx(jpoint) = tmp * mu_der(1)
|
||||
! grady(jpoint) = tmp * mu_der(2)
|
||||
! gradz(jpoint) = tmp * mu_der(3)
|
||||
!
|
||||
! if(r12 .lt. 1d-10) then
|
||||
! gradx(jpoint) = 0.d0
|
||||
! grady(jpoint) = 0.d0
|
||||
! gradz(jpoint) = 0.d0
|
||||
! cycle
|
||||
! endif
|
||||
!
|
||||
! tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12
|
||||
!
|
||||
! gradx(jpoint) = gradx(jpoint) + tmp * dx
|
||||
! grady(jpoint) = grady(jpoint) + tmp * dy
|
||||
! gradz(jpoint) = gradz(jpoint) + tmp * dz
|
||||
enddo
|
||||
|
||||
elseif(j2e_type .eq. "Boys") then
|
||||
|
||||
@ -127,8 +127,8 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
|
||||
! TODO combine 1shot & int2_grad1_u12_ao_num
|
||||
PROVIDE int2_grad1_u12_ao_num
|
||||
int2_grad1_u12_ao = int2_grad1_u12_ao_num
|
||||
!PROVIDE int2_grad1_u12_ao_num_1shot
|
||||
!int2_grad1_u12_ao = int2_grad1_u12_ao_num_1shot
|
||||
! PROVIDE int2_grad1_u12_ao_num_1shot
|
||||
! int2_grad1_u12_ao = int2_grad1_u12_ao_num_1shot
|
||||
endif
|
||||
|
||||
elseif(tc_integ_type .eq. "semi-analytic") then
|
||||
|
||||
@ -230,7 +230,7 @@ default: 70
|
||||
type: character*(32)
|
||||
doc: approach used to evaluate TC integrals [ analytic | numeric | semi-analytic ]
|
||||
interface: ezfio,ocaml,provider
|
||||
default: semi-analytic
|
||||
default: numeric
|
||||
|
||||
[minimize_lr_angles]
|
||||
type: logical
|
||||
|
||||
@ -146,26 +146,18 @@ subroutine mu_of_r_mean_field(r,mu_mf, dm)
|
||||
endif
|
||||
end
|
||||
|
||||
subroutine grad_mu_of_r_mean_field(r,mu_mf, dm, grad_mu_mf, grad_dm)
|
||||
subroutine mu_of_r_mean_field_tc(r,mu_mf, dm)
|
||||
implicit none
|
||||
include 'constants.include.F'
|
||||
double precision, intent(in) :: r(3)
|
||||
double precision, intent(out):: grad_mu_mf(3), grad_dm(3)
|
||||
double precision, intent(out):: mu_mf, dm
|
||||
double precision :: grad_f_mf_ab(3), grad_two_bod_dens(3),grad_dm_a(3), grad_dm_b(3)
|
||||
double precision :: f_mf_ab,two_bod_dens, dm_a, dm_b
|
||||
call 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)
|
||||
|
||||
call get_f_mf_ab(r,f_mf_ab,two_bod_dens, dm_a, dm_b)
|
||||
dm = dm_a + dm_b
|
||||
grad_dm(1:3) = grad_dm_a(1:3) + grad_dm_b(1:3)
|
||||
|
||||
if(dabs(two_bod_dens).lt.1.d-10)then
|
||||
mu_mf = 1.d+10
|
||||
grad_mu_mf = 0.d0
|
||||
else
|
||||
mu_mf = 0.5d0 * sqpi * f_mf_ab/two_bod_dens
|
||||
grad_mu_mf(1:3) = 0.5d0 * sqpi * (grad_f_mf_ab(1:3) * two_bod_dens - f_mf_ab * grad_two_bod_dens(1:3))&
|
||||
/(two_bod_dens*two_bod_dens)
|
||||
mu_mf = 0.3333333333d0 * sqpi * (f_mf_ab/two_bod_dens + 0.25d0)
|
||||
endif
|
||||
|
||||
end
|
||||
|
||||
|
||||
@ -20,7 +20,8 @@ program projected_operators
|
||||
! call test
|
||||
! call test_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
|
||||
end
|
||||
|
||||
|
||||
@ -174,3 +175,33 @@ subroutine test_grad_mu_mf
|
||||
print*, accu_grad_mu_mf
|
||||
|
||||
end
|
||||
|
||||
subroutine plot_mu_of_r_mf
|
||||
implicit none
|
||||
include 'constants.include.F'
|
||||
integer :: ipoint,npoint
|
||||
double precision :: dx,r(3),xmax,xmin
|
||||
double precision :: accu_mu,accu_nelec,mu_mf, dm,mu_mf_tc
|
||||
character*(128) :: output
|
||||
integer :: i_unit_output,getUnitAndOpen
|
||||
output=trim(ezfio_filename)//'.mu_mf'
|
||||
i_unit_output = getUnitAndOpen(output,'w')
|
||||
xmax = 5.D0
|
||||
xmin = 0.d0
|
||||
npoint = 10000
|
||||
dx = (xmax - xmin)/dble(npoint)
|
||||
r = 0.d0
|
||||
r(1) = xmin
|
||||
accu_mu = 0.d0
|
||||
accu_nelec = 0.d0
|
||||
do ipoint = 1, npoint
|
||||
call mu_of_r_mean_field(r,mu_mf, dm)
|
||||
call mu_of_r_mean_field_tc(r,mu_mf_tc, dm)
|
||||
write(i_unit_output,'(100(F16.10,X))')r(1),mu_mf,mu_mf_tc,dm
|
||||
accu_mu += mu_mf * dm * r(1)**2*dx*4.D0*pi
|
||||
accu_nelec += dm * r(1)**2*dx*4.D0*pi
|
||||
r(1) += dx
|
||||
enddo
|
||||
print*,'nelec = ',accu_nelec
|
||||
print*,'mu average = ',accu_mu/accu_nelec
|
||||
end
|
||||
|
||||
Loading…
Reference in New Issue
Block a user