mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-12-23 20:13:42 +01:00
91 lines
2.6 KiB
Fortran
91 lines
2.6 KiB
Fortran
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
|