mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-12-13 07:03:29 +01:00
non_h_ints compiles
This commit is contained in:
parent
17d8197a67
commit
4472a6d9be
@ -46,142 +46,327 @@ double precision function NAI_pol_mult_erf_ao(i_ao,j_ao,mu_in,C_center)
|
||||
|
||||
end
|
||||
|
||||
double precision function NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in, mu_in)
|
||||
|
||||
|
||||
double precision function NAI_pol_mult_erf(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in,mu_in)
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Computes the following integral :
|
||||
!
|
||||
! .. math::
|
||||
!
|
||||
!
|
||||
! \int dr (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 )
|
||||
! \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
|
||||
! \frac{\erf(\mu |r - R_C |)}{| r - R_C |}$.
|
||||
!
|
||||
END_DOC
|
||||
|
||||
include 'utils/constants.include.F'
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: n_pt_in
|
||||
integer, intent(in) :: power_A(3), power_B(3)
|
||||
double precision, intent(in) :: C_center(3), A_center(3), B_center(3), alpha, beta, mu_in
|
||||
|
||||
integer :: i, n_pt, n_pt_out
|
||||
double precision :: P_center(3)
|
||||
double precision :: d(0:n_pt_in), coeff, dist, const, factor
|
||||
double precision :: const_factor, dist_integral
|
||||
double precision :: accu, p_inv, p, rho, p_inv_2
|
||||
double precision :: p_new
|
||||
|
||||
double precision :: rint
|
||||
|
||||
p = alpha + beta
|
||||
p_inv = 1.d0 / p
|
||||
p_inv_2 = 0.5d0 * p_inv
|
||||
rho = alpha * beta * p_inv
|
||||
|
||||
dist = 0.d0
|
||||
dist_integral = 0.d0
|
||||
do i = 1, 3
|
||||
P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv
|
||||
dist += (A_center(i) - B_center(i)) * (A_center(i) - B_center(i))
|
||||
dist_integral += (P_center(i) - C_center(i)) * (P_center(i) - C_center(i))
|
||||
enddo
|
||||
const_factor = dist * rho
|
||||
if(const_factor > 80.d0) then
|
||||
NAI_pol_mult_erf = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
p_new = mu_in / dsqrt(p + mu_in * mu_in)
|
||||
factor = dexp(-const_factor)
|
||||
coeff = dtwo_pi * factor * p_inv * p_new
|
||||
|
||||
n_pt = 2 * ( (power_A(1) + power_B(1)) + (power_A(2) + power_B(2)) + (power_A(3) + power_B(3)) )
|
||||
const = p * dist_integral * p_new * p_new
|
||||
if(n_pt == 0) then
|
||||
NAI_pol_mult_erf = coeff * rint(0, const)
|
||||
return
|
||||
endif
|
||||
|
||||
do i = 0, n_pt_in
|
||||
d(i) = 0.d0
|
||||
enddo
|
||||
! call give_polynomial_mult_center_one_e_erf(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in)
|
||||
p_new = p_new * p_new
|
||||
call give_polynomial_mult_center_one_e_erf_opt(A_center, B_center, power_A, power_B, C_center, n_pt_in, d, n_pt_out, p_inv_2, p_new, P_center)
|
||||
|
||||
if(n_pt_out < 0) then
|
||||
NAI_pol_mult_erf = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i
|
||||
accu = 0.d0
|
||||
do i = 0, n_pt_out, 2
|
||||
accu += d(i) * rint(i/2, const)
|
||||
enddo
|
||||
NAI_pol_mult_erf = accu * coeff
|
||||
|
||||
end function NAI_pol_mult_erf
|
||||
|
||||
! ---
|
||||
|
||||
|
||||
double precision function NAI_pol_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center)
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Computes the following integral :
|
||||
! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$.
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: n_pt_in
|
||||
double precision,intent(in) :: C_center(3),A_center(3),B_center(3),alpha,beta,mu_in
|
||||
integer, intent(in) :: power_A(3),power_B(3)
|
||||
integer :: i,j,k,l,n_pt
|
||||
double precision :: P_center(3)
|
||||
integer, intent(in) :: i_ao, j_ao
|
||||
double precision, intent(in) :: beta, B_center(3)
|
||||
double precision, intent(in) :: mu_in, C_center(3)
|
||||
|
||||
integer :: i, j, power_A1(3), power_A2(3), n_pt_in
|
||||
double precision :: A1_center(3), A2_center(3), alpha1, alpha2, coef12, coef1, integral
|
||||
|
||||
double precision, external :: NAI_pol_mult_erf_with1s, NAI_pol_mult_erf_ao
|
||||
|
||||
ASSERT(beta .ge. 0.d0)
|
||||
if(beta .lt. 1d-10) then
|
||||
NAI_pol_mult_erf_ao_with1s = NAI_pol_mult_erf_ao(i_ao, j_ao, mu_in, C_center)
|
||||
return
|
||||
endif
|
||||
|
||||
power_A1(1:3) = ao_power(i_ao,1:3)
|
||||
power_A2(1:3) = ao_power(j_ao,1:3)
|
||||
|
||||
A1_center(1:3) = nucl_coord(ao_nucl(i_ao),1:3)
|
||||
A2_center(1:3) = nucl_coord(ao_nucl(j_ao),1:3)
|
||||
|
||||
n_pt_in = n_pt_max_integrals
|
||||
|
||||
NAI_pol_mult_erf_ao_with1s = 0.d0
|
||||
do i = 1, ao_prim_num(i_ao)
|
||||
alpha1 = ao_expo_ordered_transp (i,i_ao)
|
||||
coef1 = ao_coef_normalized_ordered_transp(i,i_ao)
|
||||
|
||||
do j = 1, ao_prim_num(j_ao)
|
||||
alpha2 = ao_expo_ordered_transp(j,j_ao)
|
||||
coef12 = coef1 * ao_coef_normalized_ordered_transp(j,j_ao)
|
||||
if(dabs(coef12) .lt. 1d-14) cycle
|
||||
|
||||
integral = NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A1, power_A2, alpha1, alpha2 &
|
||||
, beta, B_center, C_center, n_pt_in, mu_in )
|
||||
|
||||
NAI_pol_mult_erf_ao_with1s += integral * coef12
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end function NAI_pol_mult_erf_ao_with1s
|
||||
|
||||
subroutine NAI_pol_mult_erf_with1s_v(A1_center, A2_center, power_A1, power_A2, alpha1, alpha2, beta, B_center, LD_B, C_center, LD_C, n_pt_in, mu_in, res_v, LD_resv, n_points)
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Computes the following integral :
|
||||
!
|
||||
! .. math ::
|
||||
!
|
||||
! \int dx (x - A1_x)^a_1 (x - B1_x)^a_2 \exp(-\alpha_1 (x - A1_x)^2 - \alpha_2 (x - A2_x)^2)
|
||||
! \int dy (y - A1_y)^b_1 (y - B1_y)^b_2 \exp(-\alpha_1 (y - A1_y)^2 - \alpha_2 (y - A2_y)^2)
|
||||
! \int dz (x - A1_z)^c_1 (z - B1_z)^c_2 \exp(-\alpha_1 (z - A1_z)^2 - \alpha_2 (z - A2_z)^2)
|
||||
! \exp(-\beta (r - B)^2)
|
||||
! \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$.
|
||||
!
|
||||
END_DOC
|
||||
|
||||
double precision :: d(0:n_pt_in),pouet,coeff,dist,const,pouet_2,factor
|
||||
double precision :: I_n_special_exact,integrate_bourrin,I_n_bibi
|
||||
double precision :: V_e_n,const_factor,dist_integral,tmp
|
||||
double precision :: accu,rint,p_inv,p,rho,p_inv_2
|
||||
integer :: n_pt_out,lmax
|
||||
include 'utils/constants.include.F'
|
||||
p = alpha + beta
|
||||
p_inv = 1.d0/p
|
||||
p_inv_2 = 0.5d0 * p_inv
|
||||
rho = alpha * beta * p_inv
|
||||
|
||||
dist = 0.d0
|
||||
dist_integral = 0.d0
|
||||
do i = 1, 3
|
||||
P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv
|
||||
dist += (A_center(i) - B_center(i))*(A_center(i) - B_center(i))
|
||||
dist_integral += (P_center(i) - C_center(i))*(P_center(i) - C_center(i))
|
||||
enddo
|
||||
const_factor = dist*rho
|
||||
if(const_factor > 80.d0)then
|
||||
NAI_pol_mult_erf = 0.d0
|
||||
return
|
||||
endif
|
||||
double precision :: p_new
|
||||
p_new = mu_in/dsqrt(p+ mu_in * mu_in)
|
||||
factor = dexp(-const_factor)
|
||||
coeff = dtwo_pi * factor * p_inv * p_new
|
||||
lmax = 20
|
||||
implicit none
|
||||
integer, intent(in) :: n_pt_in, LD_B, LD_C, LD_resv, n_points
|
||||
integer, intent(in) :: power_A1(3), power_A2(3)
|
||||
double precision, intent(in) :: A1_center(3), A2_center(3)
|
||||
double precision, intent(in) :: C_center(LD_C,3), B_center(LD_B,3)
|
||||
double precision, intent(in) :: alpha1, alpha2, beta, mu_in
|
||||
double precision, intent(out) :: res_v(LD_resv)
|
||||
|
||||
! print*, "b"
|
||||
do i = 0, n_pt_in
|
||||
d(i) = 0.d0
|
||||
enddo
|
||||
n_pt = 2 * ( (power_A(1) + power_B(1)) +(power_A(2) + power_B(2)) +(power_A(3) + power_B(3)) )
|
||||
const = p * dist_integral * p_new * p_new
|
||||
if (n_pt == 0) then
|
||||
pouet = rint(0,const)
|
||||
NAI_pol_mult_erf = coeff * pouet
|
||||
integer :: i, n_pt, n_pt_out, ipoint
|
||||
double precision :: alpha12, alpha12_inv, alpha12_inv_2, rho12, A12_center(3), dist12, const_factor12
|
||||
double precision :: p, p_inv, p_inv_2, rho, P_center(3), dist, const_factor
|
||||
double precision :: dist_integral
|
||||
double precision :: d(0:n_pt_in), coeff, const, factor
|
||||
double precision :: accu
|
||||
double precision :: p_new, p_new2, coef_tmp, cons_tmp
|
||||
|
||||
double precision :: rint
|
||||
|
||||
|
||||
res_V(1:LD_resv) = 0.d0
|
||||
|
||||
! e^{-alpha1 (r - A1)^2} e^{-alpha2 (r - A2)^2} = e^{-K12} e^{-alpha12 (r - A12)^2}
|
||||
alpha12 = alpha1 + alpha2
|
||||
alpha12_inv = 1.d0 / alpha12
|
||||
alpha12_inv_2 = 0.5d0 * alpha12_inv
|
||||
rho12 = alpha1 * alpha2 * alpha12_inv
|
||||
A12_center(1) = (alpha1 * A1_center(1) + alpha2 * A2_center(1)) * alpha12_inv
|
||||
A12_center(2) = (alpha1 * A1_center(2) + alpha2 * A2_center(2)) * alpha12_inv
|
||||
A12_center(3) = (alpha1 * A1_center(3) + alpha2 * A2_center(3)) * alpha12_inv
|
||||
dist12 = (A1_center(1) - A2_center(1)) * (A1_center(1) - A2_center(1))&
|
||||
+ (A1_center(2) - A2_center(2)) * (A1_center(2) - A2_center(2))&
|
||||
+ (A1_center(3) - A2_center(3)) * (A1_center(3) - A2_center(3))
|
||||
|
||||
const_factor12 = dist12 * rho12
|
||||
if(const_factor12 > 80.d0) then
|
||||
return
|
||||
endif
|
||||
|
||||
! call give_polynomial_mult_center_one_e_erf(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in)
|
||||
p_new = p_new * p_new
|
||||
call give_polynomial_mult_center_one_e_erf_opt(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in,p,p_inv,p_inv_2,p_new,P_center)
|
||||
! e^{-K12} e^{-alpha12 (r - A12)^2} e^{-beta (r - B)^2} = e^{-K} e^{-p (r - P)^2}
|
||||
p = alpha12 + beta
|
||||
p_inv = 1.d0 / p
|
||||
p_inv_2 = 0.5d0 * p_inv
|
||||
rho = alpha12 * beta * p_inv
|
||||
p_new = mu_in / dsqrt(p + mu_in * mu_in)
|
||||
p_new2 = p_new * p_new
|
||||
coef_tmp = dtwo_pi * p_inv * p_new
|
||||
cons_tmp = p * p_new2
|
||||
n_pt = 2 * (power_A1(1) + power_A2(1) + power_A1(2) + power_A2(2) + power_A1(3) + power_A2(3) )
|
||||
|
||||
if(n_pt == 0) then
|
||||
|
||||
do ipoint = 1, n_points
|
||||
|
||||
dist = (A12_center(1) - B_center(ipoint,1)) * (A12_center(1) - B_center(ipoint,1))&
|
||||
+ (A12_center(2) - B_center(ipoint,2)) * (A12_center(2) - B_center(ipoint,2))&
|
||||
+ (A12_center(3) - B_center(ipoint,3)) * (A12_center(3) - B_center(ipoint,3))
|
||||
const_factor = const_factor12 + dist * rho
|
||||
if(const_factor > 80.d0) cycle
|
||||
coeff = coef_tmp * dexp(-const_factor)
|
||||
|
||||
P_center(1) = (alpha12 * A12_center(1) + beta * B_center(ipoint,1)) * p_inv
|
||||
P_center(2) = (alpha12 * A12_center(2) + beta * B_center(ipoint,2)) * p_inv
|
||||
P_center(3) = (alpha12 * A12_center(3) + beta * B_center(ipoint,3)) * p_inv
|
||||
dist_integral = (P_center(1) - C_center(ipoint,1)) * (P_center(1) - C_center(ipoint,1))&
|
||||
+ (P_center(2) - C_center(ipoint,2)) * (P_center(2) - C_center(ipoint,2))&
|
||||
+ (P_center(3) - C_center(ipoint,3)) * (P_center(3) - C_center(ipoint,3))
|
||||
const = cons_tmp * dist_integral
|
||||
|
||||
res_v(ipoint) = coeff * rint(0, const)
|
||||
enddo
|
||||
|
||||
else
|
||||
|
||||
do ipoint = 1, n_points
|
||||
|
||||
dist = (A12_center(1) - B_center(ipoint,1)) * (A12_center(1) - B_center(ipoint,1))&
|
||||
+ (A12_center(2) - B_center(ipoint,2)) * (A12_center(2) - B_center(ipoint,2))&
|
||||
+ (A12_center(3) - B_center(ipoint,3)) * (A12_center(3) - B_center(ipoint,3))
|
||||
const_factor = const_factor12 + dist * rho
|
||||
if(const_factor > 80.d0) cycle
|
||||
coeff = coef_tmp * dexp(-const_factor)
|
||||
|
||||
P_center(1) = (alpha12 * A12_center(1) + beta * B_center(ipoint,1)) * p_inv
|
||||
P_center(2) = (alpha12 * A12_center(2) + beta * B_center(ipoint,2)) * p_inv
|
||||
P_center(3) = (alpha12 * A12_center(3) + beta * B_center(ipoint,3)) * p_inv
|
||||
dist_integral = (P_center(1) - C_center(ipoint,1)) * (P_center(1) - C_center(ipoint,1))&
|
||||
+ (P_center(2) - C_center(ipoint,2)) * (P_center(2) - C_center(ipoint,2))&
|
||||
+ (P_center(3) - C_center(ipoint,3)) * (P_center(3) - C_center(ipoint,3))
|
||||
const = cons_tmp * dist_integral
|
||||
|
||||
do i = 0, n_pt_in
|
||||
d(i) = 0.d0
|
||||
enddo
|
||||
!TODO: VECTORIZE HERE
|
||||
call give_polynomial_mult_center_one_e_erf_opt(A1_center, A2_center, power_A1, power_A2, C_center(ipoint,1:3), n_pt_in, d, n_pt_out, p_inv_2, p_new2, P_center)
|
||||
|
||||
if(n_pt_out < 0) then
|
||||
cycle
|
||||
endif
|
||||
|
||||
! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i
|
||||
accu = 0.d0
|
||||
do i = 0, n_pt_out, 2
|
||||
accu += d(i) * rint(i/2, const)
|
||||
enddo
|
||||
|
||||
res_v(ipoint) = accu * coeff
|
||||
enddo
|
||||
|
||||
if(n_pt_out<0)then
|
||||
NAI_pol_mult_erf = 0.d0
|
||||
return
|
||||
endif
|
||||
accu = 0.d0
|
||||
|
||||
! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i
|
||||
do i =0 ,n_pt_out,2
|
||||
accu += d(i) * rint(i/2,const)
|
||||
enddo
|
||||
NAI_pol_mult_erf = accu * coeff
|
||||
end subroutine NAI_pol_mult_erf_with1s_v
|
||||
|
||||
end
|
||||
! ---
|
||||
|
||||
subroutine give_polynomial_mult_center_one_e_erf_opt(A_center, B_center, power_A, power_B, C_center, n_pt_in, d, n_pt_out, p_inv_2, p_new, P_center)
|
||||
|
||||
subroutine give_polynomial_mult_center_one_e_erf_opt(A_center,B_center,alpha,beta,&
|
||||
power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in,p,p_inv,p_inv_2,p_new,P_center)
|
||||
BEGIN_DOC
|
||||
! Returns the explicit polynomial in terms of the $t$ variable of the
|
||||
! following polynomial:
|
||||
!
|
||||
! $I_{x1}(a_x, d_x,p,q) \times I_{x1}(a_y, d_y,p,q) \times I_{x1}(a_z, d_z,p,q)$.
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: n_pt_in
|
||||
integer,intent(out) :: n_pt_out
|
||||
double precision, intent(in) :: A_center(3), B_center(3),C_center(3),p,p_inv,p_inv_2,p_new,P_center(3)
|
||||
double precision, intent(in) :: alpha,beta,mu_in
|
||||
integer, intent(in) :: power_A(3), power_B(3)
|
||||
integer :: a_x,b_x,a_y,b_y,a_z,b_z
|
||||
double precision :: d(0:n_pt_in)
|
||||
double precision :: d1(0:n_pt_in)
|
||||
double precision :: d2(0:n_pt_in)
|
||||
double precision :: d3(0:n_pt_in)
|
||||
double precision :: accu
|
||||
integer, intent(in) :: n_pt_in
|
||||
integer, intent(in) :: power_A(3), power_B(3)
|
||||
double precision, intent(in) :: A_center(3), B_center(3), C_center(3), p_inv_2, p_new, P_center(3)
|
||||
integer, intent(out) :: n_pt_out
|
||||
double precision, intent(out) :: d(0:n_pt_in)
|
||||
|
||||
integer :: a_x, b_x, a_y, b_y, a_z, b_z
|
||||
integer :: n_pt1, n_pt2, n_pt3, dim, i
|
||||
integer :: n_pt_tmp
|
||||
double precision :: d1(0:n_pt_in)
|
||||
double precision :: d2(0:n_pt_in)
|
||||
double precision :: d3(0:n_pt_in)
|
||||
double precision :: accu
|
||||
double precision :: R1x(0:2), B01(0:2), R1xp(0:2), R2x(0:2)
|
||||
|
||||
accu = 0.d0
|
||||
ASSERT (n_pt_in > 1)
|
||||
|
||||
double precision :: R1x(0:2), B01(0:2), R1xp(0:2),R2x(0:2)
|
||||
R1x(0) = (P_center(1) - A_center(1))
|
||||
R1x(1) = 0.d0
|
||||
R1x(2) = -(P_center(1) - C_center(1))* p_new
|
||||
R1x(0) = (P_center(1) - A_center(1))
|
||||
R1x(1) = 0.d0
|
||||
R1x(2) = -(P_center(1) - C_center(1))* p_new
|
||||
! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
|
||||
R1xp(0) = (P_center(1) - B_center(1))
|
||||
R1xp(1) = 0.d0
|
||||
R1xp(2) =-(P_center(1) - C_center(1))* p_new
|
||||
R1xp(0) = (P_center(1) - B_center(1))
|
||||
R1xp(1) = 0.d0
|
||||
R1xp(2) =-(P_center(1) - C_center(1))* p_new
|
||||
!R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
|
||||
R2x(0) = p_inv_2
|
||||
R2x(1) = 0.d0
|
||||
R2x(2) = -p_inv_2* p_new
|
||||
R2x(0) = p_inv_2
|
||||
R2x(1) = 0.d0
|
||||
R2x(2) = -p_inv_2 * p_new
|
||||
!R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2
|
||||
do i = 0,n_pt_in
|
||||
d(i) = 0.d0
|
||||
enddo
|
||||
do i = 0,n_pt_in
|
||||
|
||||
do i = 0, n_pt_in
|
||||
d (i) = 0.d0
|
||||
d1(i) = 0.d0
|
||||
enddo
|
||||
do i = 0,n_pt_in
|
||||
d2(i) = 0.d0
|
||||
enddo
|
||||
do i = 0,n_pt_in
|
||||
d3(i) = 0.d0
|
||||
enddo
|
||||
integer :: n_pt1,n_pt2,n_pt3,dim,i
|
||||
|
||||
n_pt1 = n_pt_in
|
||||
n_pt2 = n_pt_in
|
||||
n_pt3 = n_pt_in
|
||||
a_x = power_A(1)
|
||||
b_x = power_B(1)
|
||||
call I_x1_pol_mult_one_e(a_x,b_x,R1x,R1xp,R2x,d1,n_pt1,n_pt_in)
|
||||
call I_x1_pol_mult_one_e(a_x, b_x, R1x, R1xp, R2x, d1, n_pt1, n_pt_in)
|
||||
if(n_pt1<0)then
|
||||
n_pt_out = -1
|
||||
do i = 0,n_pt_in
|
||||
@ -190,17 +375,17 @@ subroutine give_polynomial_mult_center_one_e_erf_opt(A_center,B_center,alpha,bet
|
||||
return
|
||||
endif
|
||||
|
||||
R1x(0) = (P_center(2) - A_center(2))
|
||||
R1x(1) = 0.d0
|
||||
R1x(2) = -(P_center(2) - C_center(2))* p_new
|
||||
R1x(0) = (P_center(2) - A_center(2))
|
||||
R1x(1) = 0.d0
|
||||
R1x(2) = -(P_center(2) - C_center(2))* p_new
|
||||
! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
|
||||
R1xp(0) = (P_center(2) - B_center(2))
|
||||
R1xp(1) = 0.d0
|
||||
R1xp(2) =-(P_center(2) - C_center(2))* p_new
|
||||
R1xp(0) = (P_center(2) - B_center(2))
|
||||
R1xp(1) = 0.d0
|
||||
R1xp(2) =-(P_center(2) - C_center(2))* p_new
|
||||
!R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
|
||||
a_y = power_A(2)
|
||||
b_y = power_B(2)
|
||||
call I_x1_pol_mult_one_e(a_y,b_y,R1x,R1xp,R2x,d2,n_pt2,n_pt_in)
|
||||
call I_x1_pol_mult_one_e(a_y, b_y, R1x, R1xp, R2x, d2, n_pt2, n_pt_in)
|
||||
if(n_pt2<0)then
|
||||
n_pt_out = -1
|
||||
do i = 0,n_pt_in
|
||||
@ -209,51 +394,151 @@ subroutine give_polynomial_mult_center_one_e_erf_opt(A_center,B_center,alpha,bet
|
||||
return
|
||||
endif
|
||||
|
||||
|
||||
R1x(0) = (P_center(3) - A_center(3))
|
||||
R1x(1) = 0.d0
|
||||
R1x(2) = -(P_center(3) - C_center(3))* p_new
|
||||
R1x(0) = (P_center(3) - A_center(3))
|
||||
R1x(1) = 0.d0
|
||||
R1x(2) = -(P_center(3) - C_center(3)) * p_new
|
||||
! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
|
||||
R1xp(0) = (P_center(3) - B_center(3))
|
||||
R1xp(1) = 0.d0
|
||||
R1xp(2) =-(P_center(3) - C_center(3))* p_new
|
||||
R1xp(0) = (P_center(3) - B_center(3))
|
||||
R1xp(1) = 0.d0
|
||||
R1xp(2) =-(P_center(3) - C_center(3)) * p_new
|
||||
!R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2
|
||||
a_z = power_A(3)
|
||||
b_z = power_B(3)
|
||||
|
||||
call I_x1_pol_mult_one_e(a_z,b_z,R1x,R1xp,R2x,d3,n_pt3,n_pt_in)
|
||||
if(n_pt3<0)then
|
||||
call I_x1_pol_mult_one_e(a_z, b_z, R1x, R1xp, R2x, d3, n_pt3, n_pt_in)
|
||||
if(n_pt3 < 0) then
|
||||
n_pt_out = -1
|
||||
do i = 0,n_pt_in
|
||||
d(i) = 0.d0
|
||||
enddo
|
||||
return
|
||||
endif
|
||||
integer :: n_pt_tmp
|
||||
|
||||
n_pt_tmp = 0
|
||||
call multiply_poly(d1,n_pt1,d2,n_pt2,d,n_pt_tmp)
|
||||
do i = 0,n_pt_tmp
|
||||
call multiply_poly(d1, n_pt1, d2, n_pt2, d, n_pt_tmp)
|
||||
do i = 0, n_pt_tmp
|
||||
d1(i) = 0.d0
|
||||
enddo
|
||||
n_pt_out = 0
|
||||
call multiply_poly(d ,n_pt_tmp ,d3,n_pt3,d1,n_pt_out)
|
||||
call multiply_poly(d, n_pt_tmp, d3, n_pt3, d1, n_pt_out)
|
||||
do i = 0, n_pt_out
|
||||
d(i) = d1(i)
|
||||
enddo
|
||||
|
||||
end
|
||||
end subroutine give_polynomial_mult_center_one_e_erf_opt
|
||||
|
||||
! ---
|
||||
subroutine NAI_pol_mult_erf_v(A_center, B_center, power_A, power_B, alpha, beta, C_center, LD_C, n_pt_in, mu_in, res_v, LD_resv, n_points)
|
||||
|
||||
|
||||
|
||||
subroutine give_polynomial_mult_center_one_e_erf(A_center,B_center,alpha,beta,&
|
||||
power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in)
|
||||
BEGIN_DOC
|
||||
! Returns the explicit polynomial in terms of the $t$ variable of the
|
||||
!
|
||||
! Computes the following integral :
|
||||
!
|
||||
! .. math::
|
||||
!
|
||||
! \int dr (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 )
|
||||
! \frac{\erf(\mu |r - R_C |)}{| r - R_C |}$.
|
||||
!
|
||||
END_DOC
|
||||
|
||||
include 'utils/constants.include.F'
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: n_pt_in, n_points, LD_C, LD_resv
|
||||
integer, intent(in) :: power_A(3), power_B(3)
|
||||
double precision, intent(in) :: A_center(3), B_center(3), alpha, beta, mu_in
|
||||
double precision, intent(in) :: C_center(LD_C,3)
|
||||
double precision, intent(out) :: res_v(LD_resv)
|
||||
|
||||
integer :: i, n_pt, n_pt_out, ipoint
|
||||
double precision :: P_center(3)
|
||||
double precision :: d(0:n_pt_in), coeff, dist, const, factor
|
||||
double precision :: const_factor, dist_integral
|
||||
double precision :: accu, p_inv, p, rho, p_inv_2
|
||||
double precision :: p_new, p_new2, coef_tmp
|
||||
|
||||
double precision :: rint
|
||||
|
||||
res_V(1:LD_resv) = 0.d0
|
||||
|
||||
p = alpha + beta
|
||||
p_inv = 1.d0 / p
|
||||
p_inv_2 = 0.5d0 * p_inv
|
||||
rho = alpha * beta * p_inv
|
||||
p_new = mu_in / dsqrt(p + mu_in * mu_in)
|
||||
p_new2 = p_new * p_new
|
||||
coef_tmp = p * p_new2
|
||||
|
||||
dist = 0.d0
|
||||
do i = 1, 3
|
||||
P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv
|
||||
dist += (A_center(i) - B_center(i)) * (A_center(i) - B_center(i))
|
||||
enddo
|
||||
|
||||
const_factor = dist * rho
|
||||
if(const_factor > 80.d0) then
|
||||
return
|
||||
endif
|
||||
factor = dexp(-const_factor)
|
||||
coeff = dtwo_pi * factor * p_inv * p_new
|
||||
|
||||
n_pt = 2 * ( power_A(1) + power_B(1) + power_A(2) + power_B(2) + power_A(3) + power_B(3) )
|
||||
|
||||
if(n_pt == 0) then
|
||||
|
||||
do ipoint = 1, n_points
|
||||
dist_integral = 0.d0
|
||||
do i = 1, 3
|
||||
dist_integral += (P_center(i) - C_center(ipoint,i)) * (P_center(i) - C_center(ipoint,i))
|
||||
enddo
|
||||
const = coef_tmp * dist_integral
|
||||
|
||||
res_v(ipoint) = coeff * rint(0, const)
|
||||
enddo
|
||||
|
||||
else
|
||||
|
||||
do ipoint = 1, n_points
|
||||
dist_integral = 0.d0
|
||||
do i = 1, 3
|
||||
dist_integral += (P_center(i) - C_center(ipoint,i)) * (P_center(i) - C_center(ipoint,i))
|
||||
enddo
|
||||
const = coef_tmp * dist_integral
|
||||
|
||||
do i = 0, n_pt_in
|
||||
d(i) = 0.d0
|
||||
enddo
|
||||
call give_polynomial_mult_center_one_e_erf_opt(A_center, B_center, power_A, power_B, C_center(ipoint,1:3), n_pt_in, d, n_pt_out, p_inv_2, p_new2, P_center)
|
||||
|
||||
if(n_pt_out < 0) then
|
||||
res_v(ipoint) = 0.d0
|
||||
cycle
|
||||
endif
|
||||
|
||||
! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i
|
||||
accu = 0.d0
|
||||
do i = 0, n_pt_out, 2
|
||||
accu += d(i) * rint(i/2, const)
|
||||
enddo
|
||||
|
||||
res_v(ipoint) = accu * coeff
|
||||
enddo
|
||||
|
||||
endif
|
||||
|
||||
end subroutine NAI_pol_mult_erf_v
|
||||
|
||||
|
||||
subroutine give_polynomial_mult_center_one_e_erf(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in)
|
||||
|
||||
BEGIN_DOC
|
||||
! Returns the explicit polynomial in terms of the $t$ variable of the
|
||||
! following polynomial:
|
||||
!
|
||||
! $I_{x1}(a_x, d_x,p,q) \times I_{x1}(a_y, d_y,p,q) \times I_{x1}(a_z, d_z,p,q)$.
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: n_pt_in
|
||||
integer,intent(out) :: n_pt_out
|
||||
@ -374,3 +659,113 @@ subroutine give_polynomial_mult_center_one_e_erf(A_center,B_center,alpha,beta,&
|
||||
|
||||
end
|
||||
|
||||
double precision function NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A1, power_A2, alpha1, alpha2 &
|
||||
, beta, B_center, C_center, n_pt_in, mu_in )
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Computes the following integral :
|
||||
!
|
||||
! .. math::
|
||||
!
|
||||
! \int dx (x - A1_x)^a_1 (x - B1_x)^a_2 \exp(-\alpha_1 (x - A1_x)^2 - \alpha_2 (x - A2_x)^2)
|
||||
! \int dy (y - A1_y)^b_1 (y - B1_y)^b_2 \exp(-\alpha_1 (y - A1_y)^2 - \alpha_2 (y - A2_y)^2)
|
||||
! \int dz (x - A1_z)^c_1 (z - B1_z)^c_2 \exp(-\alpha_1 (z - A1_z)^2 - \alpha_2 (z - A2_z)^2)
|
||||
! \exp(-\beta (r - B)^2)
|
||||
! \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$.
|
||||
!
|
||||
END_DOC
|
||||
|
||||
include 'utils/constants.include.F'
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: n_pt_in
|
||||
integer, intent(in) :: power_A1(3), power_A2(3)
|
||||
double precision, intent(in) :: C_center(3), A1_center(3), A2_center(3), B_center(3)
|
||||
double precision, intent(in) :: alpha1, alpha2, beta, mu_in
|
||||
|
||||
integer :: i, n_pt, n_pt_out
|
||||
double precision :: alpha12, alpha12_inv, alpha12_inv_2, rho12, A12_center(3), dist12, const_factor12
|
||||
double precision :: p, p_inv, p_inv_2, rho, P_center(3), dist, const_factor
|
||||
double precision :: dist_integral
|
||||
double precision :: d(0:n_pt_in), coeff, const, factor
|
||||
double precision :: accu
|
||||
double precision :: p_new
|
||||
|
||||
double precision :: rint
|
||||
|
||||
|
||||
! e^{-alpha1 (r - A1)^2} e^{-alpha2 (r - A2)^2} = e^{-K12} e^{-alpha12 (r - A12)^2}
|
||||
alpha12 = alpha1 + alpha2
|
||||
alpha12_inv = 1.d0 / alpha12
|
||||
alpha12_inv_2 = 0.5d0 * alpha12_inv
|
||||
rho12 = alpha1 * alpha2 * alpha12_inv
|
||||
A12_center(1) = (alpha1 * A1_center(1) + alpha2 * A2_center(1)) * alpha12_inv
|
||||
A12_center(2) = (alpha1 * A1_center(2) + alpha2 * A2_center(2)) * alpha12_inv
|
||||
A12_center(3) = (alpha1 * A1_center(3) + alpha2 * A2_center(3)) * alpha12_inv
|
||||
dist12 = (A1_center(1) - A2_center(1)) * (A1_center(1) - A2_center(1)) &
|
||||
+ (A1_center(2) - A2_center(2)) * (A1_center(2) - A2_center(2)) &
|
||||
+ (A1_center(3) - A2_center(3)) * (A1_center(3) - A2_center(3))
|
||||
|
||||
const_factor12 = dist12 * rho12
|
||||
if(const_factor12 > 80.d0) then
|
||||
NAI_pol_mult_erf_with1s = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
! ---
|
||||
|
||||
! e^{-K12} e^{-alpha12 (r - A12)^2} e^{-beta (r - B)^2} = e^{-K} e^{-p (r - P)^2}
|
||||
p = alpha12 + beta
|
||||
p_inv = 1.d0 / p
|
||||
p_inv_2 = 0.5d0 * p_inv
|
||||
rho = alpha12 * beta * p_inv
|
||||
P_center(1) = (alpha12 * A12_center(1) + beta * B_center(1)) * p_inv
|
||||
P_center(2) = (alpha12 * A12_center(2) + beta * B_center(2)) * p_inv
|
||||
P_center(3) = (alpha12 * A12_center(3) + beta * B_center(3)) * p_inv
|
||||
dist = (A12_center(1) - B_center(1)) * (A12_center(1) - B_center(1)) &
|
||||
+ (A12_center(2) - B_center(2)) * (A12_center(2) - B_center(2)) &
|
||||
+ (A12_center(3) - B_center(3)) * (A12_center(3) - B_center(3))
|
||||
|
||||
const_factor = const_factor12 + dist * rho
|
||||
if(const_factor > 80.d0) then
|
||||
NAI_pol_mult_erf_with1s = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
dist_integral = (P_center(1) - C_center(1)) * (P_center(1) - C_center(1)) &
|
||||
+ (P_center(2) - C_center(2)) * (P_center(2) - C_center(2)) &
|
||||
+ (P_center(3) - C_center(3)) * (P_center(3) - C_center(3))
|
||||
|
||||
! ---
|
||||
|
||||
p_new = mu_in / dsqrt(p + mu_in * mu_in)
|
||||
factor = dexp(-const_factor)
|
||||
coeff = dtwo_pi * factor * p_inv * p_new
|
||||
|
||||
n_pt = 2 * ( (power_A1(1) + power_A2(1)) + (power_A1(2) + power_A2(2)) + (power_A1(3) + power_A2(3)) )
|
||||
const = p * dist_integral * p_new * p_new
|
||||
if(n_pt == 0) then
|
||||
NAI_pol_mult_erf_with1s = coeff * rint(0, const)
|
||||
return
|
||||
endif
|
||||
|
||||
do i = 0, n_pt_in
|
||||
d(i) = 0.d0
|
||||
enddo
|
||||
p_new = p_new * p_new
|
||||
call give_polynomial_mult_center_one_e_erf_opt( A1_center, A2_center, power_A1, power_A2, C_center, n_pt_in, d, n_pt_out, p_inv_2, p_new, P_center)
|
||||
|
||||
if(n_pt_out < 0) then
|
||||
NAI_pol_mult_erf_with1s = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i
|
||||
accu = 0.d0
|
||||
do i = 0, n_pt_out, 2
|
||||
accu += d(i) * rint(i/2, const)
|
||||
enddo
|
||||
NAI_pol_mult_erf_with1s = accu * coeff
|
||||
|
||||
end function NAI_pol_mult_erf_with1s
|
||||
|
5
src/ao_tc_eff_map/NEED
Normal file
5
src/ao_tc_eff_map/NEED
Normal file
@ -0,0 +1,5 @@
|
||||
ao_two_e_erf_ints
|
||||
mo_one_e_ints
|
||||
ao_many_one_e_ints
|
||||
dft_utils_in_r
|
||||
tc_keywords
|
12
src/ao_tc_eff_map/README.rst
Normal file
12
src/ao_tc_eff_map/README.rst
Normal file
@ -0,0 +1,12 @@
|
||||
ao_tc_eff_map
|
||||
=============
|
||||
|
||||
This is a module to obtain the integrals on the AO basis of the SCALAR HERMITIAN
|
||||
effective potential defined in Eq. 32 of JCP 154, 084119 (2021)
|
||||
It also contains the modification by a one-body Jastrow factor.
|
||||
|
||||
The main routine/providers are
|
||||
|
||||
+) ao_tc_sym_two_e_pot_map : map of the SCALAR PART of total effective two-electron on the AO basis in PHYSICIST notations. It might contain the two-electron term coming from the one-e correlation factor.
|
||||
+) get_ao_tc_sym_two_e_pot(i,j,k,l,ao_tc_sym_two_e_pot_map) : routine to get the integrals from ao_tc_sym_two_e_pot_map.
|
||||
+) ao_tc_sym_two_e_pot(i,j,k,l) : FUNCTION that returns the scalar part of TC-potential EXCLUDING the erf(mu r12)/r12. See two_e_ints_gauss.irp.f for more details.
|
76
src/ao_tc_eff_map/compute_ints_eff_pot.irp.f
Normal file
76
src/ao_tc_eff_map/compute_ints_eff_pot.irp.f
Normal file
@ -0,0 +1,76 @@
|
||||
|
||||
|
||||
subroutine compute_ao_tc_sym_two_e_pot_jl(j, l, n_integrals, buffer_i, buffer_value)
|
||||
|
||||
use map_module
|
||||
|
||||
BEGIN_DOC
|
||||
! Parallel client for AO integrals
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: j, l
|
||||
integer,intent(out) :: n_integrals
|
||||
integer(key_kind),intent(out) :: buffer_i(ao_num*ao_num)
|
||||
real(integral_kind),intent(out) :: buffer_value(ao_num*ao_num)
|
||||
|
||||
integer :: i, k
|
||||
integer :: kk, m, j1, i1
|
||||
double precision :: cpu_1, cpu_2, wall_1, wall_2
|
||||
double precision :: integral, wall_0, integral_pot, integral_erf
|
||||
double precision :: thr
|
||||
|
||||
logical, external :: ao_two_e_integral_zero
|
||||
double precision :: ao_tc_sym_two_e_pot, ao_two_e_integral_erf
|
||||
double precision :: j1b_gauss_2e_j1, j1b_gauss_2e_j2
|
||||
|
||||
|
||||
PROVIDE j1b_type
|
||||
|
||||
thr = ao_integrals_threshold
|
||||
|
||||
n_integrals = 0
|
||||
|
||||
j1 = j+ishft(l*l-l,-1)
|
||||
do k = 1, ao_num ! r1
|
||||
i1 = ishft(k*k-k,-1)
|
||||
if (i1 > j1) then
|
||||
exit
|
||||
endif
|
||||
do i = 1, k
|
||||
i1 += 1
|
||||
if (i1 > j1) then
|
||||
exit
|
||||
endif
|
||||
|
||||
if (ao_two_e_integral_erf_schwartz(i,k)*ao_two_e_integral_erf_schwartz(j,l) < thr) then
|
||||
cycle
|
||||
endif
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
integral_pot = ao_tc_sym_two_e_pot (i, k, j, l) ! i,k : r1 j,l : r2
|
||||
integral_erf = ao_two_e_integral_erf(i, k, j, l)
|
||||
integral = integral_erf + integral_pot
|
||||
|
||||
if( j1b_type .eq. 1 ) then
|
||||
!print *, ' j1b type 1 is added'
|
||||
integral = integral + j1b_gauss_2e_j1(i, k, j, l)
|
||||
elseif( j1b_type .eq. 2 ) then
|
||||
!print *, ' j1b type 2 is added'
|
||||
integral = integral + j1b_gauss_2e_j2(i, k, j, l)
|
||||
endif
|
||||
|
||||
if(abs(integral) < thr) then
|
||||
cycle
|
||||
endif
|
||||
|
||||
n_integrals += 1
|
||||
!DIR$ FORCEINLINE
|
||||
call two_e_integrals_index(i, j, k, l, buffer_i(n_integrals))
|
||||
buffer_value(n_integrals) = integral
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end subroutine compute_ao_tc_sym_two_e_pot_jl
|
||||
|
510
src/ao_tc_eff_map/fit_j.irp.f
Normal file
510
src/ao_tc_eff_map/fit_j.irp.f
Normal file
@ -0,0 +1,510 @@
|
||||
BEGIN_PROVIDER [ double precision, expo_j_xmu_1gauss ]
|
||||
&BEGIN_PROVIDER [ double precision, coef_j_xmu_1gauss ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Upper bound long range fit of F(x) = x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2)
|
||||
!
|
||||
! with a single gaussian.
|
||||
!
|
||||
! Such a function can be used to screen integrals with F(x).
|
||||
END_DOC
|
||||
expo_j_xmu_1gauss = 0.5d0
|
||||
coef_j_xmu_1gauss = 1.d0
|
||||
END_PROVIDER
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, expo_erfc_gauss ]
|
||||
implicit none
|
||||
expo_erfc_gauss = 1.41211d0
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, expo_erfc_mu_gauss ]
|
||||
implicit none
|
||||
expo_erfc_mu_gauss = expo_erfc_gauss * mu_erf * mu_erf
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, expo_good_j_mu_1gauss ]
|
||||
&BEGIN_PROVIDER [ double precision, coef_good_j_mu_1gauss ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! exponent of Gaussian in order to obtain an upper bound of J(r12,mu)
|
||||
!
|
||||
! Can be used to scree integrals with J(r12,mu)
|
||||
END_DOC
|
||||
expo_good_j_mu_1gauss = 2.D0 * mu_erf * expo_j_xmu_1gauss
|
||||
coef_good_j_mu_1gauss = 0.5d0/mu_erf * coef_j_xmu_1gauss
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, expo_j_xmu, (n_fit_1_erf_x) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! F(x) = x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2) is fitted with a gaussian and a Slater
|
||||
!
|
||||
! \approx - 1/sqrt(pi) * exp(-alpha * x ) exp(-beta * x**2)
|
||||
!
|
||||
! where alpha = expo_j_xmu(1) and beta = expo_j_xmu(2)
|
||||
END_DOC
|
||||
expo_j_xmu(1) = 1.7477d0
|
||||
expo_j_xmu(2) = 0.668662d0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [double precision, expo_gauss_j_mu_x, (ng_fit_jast)]
|
||||
&BEGIN_PROVIDER [double precision, coef_gauss_j_mu_x, (ng_fit_jast)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! J(mu,r12) = 1/2 r12 * (1 - erf(mu*r12)) - 1/(2 sqrt(pi)*mu) exp(-(mu*r12)^2) is expressed as
|
||||
!
|
||||
! J(mu,r12) = 0.5/mu * F(r12*mu) where F(x) = x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2)
|
||||
!
|
||||
! F(x) is fitted by - 1/sqrt(pi) * exp(-alpha * x) exp(-beta * x^2) (see expo_j_xmu)
|
||||
!
|
||||
! The slater function exp(-alpha * x) is fitted with n_max_fit_slat gaussians
|
||||
!
|
||||
! See Appendix 2 of JCP 154, 084119 (2021)
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i
|
||||
double precision :: tmp
|
||||
double precision :: expos(ng_fit_jast), alpha, beta
|
||||
|
||||
if(ng_fit_jast .eq. 1) then
|
||||
|
||||
coef_gauss_j_mu_x = (/ -0.47947881d0 /)
|
||||
expo_gauss_j_mu_x = (/ 3.4987848d0 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_x(i) = tmp * expo_gauss_j_mu_x(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 2) then
|
||||
|
||||
coef_gauss_j_mu_x = (/ -0.18390742d0, -0.35512656d0 /)
|
||||
expo_gauss_j_mu_x = (/ 31.9279947d0 , 2.11428789d0 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_x(i) = tmp * expo_gauss_j_mu_x(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 3) then
|
||||
|
||||
coef_gauss_j_mu_x = (/ -0.07501725d0, -0.28499012d0, -0.1953932d0 /)
|
||||
expo_gauss_j_mu_x = (/ 206.74058566d0, 1.72974157d0, 11.18735164d0 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_x(i) = tmp * expo_gauss_j_mu_x(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 5) then
|
||||
|
||||
coef_gauss_j_mu_x = (/ -0.01832955d0 , -0.10188952d0 , -0.20710858d0 , -0.18975032d0 , -0.04641657d0 /)
|
||||
expo_gauss_j_mu_x = (/ 4.33116687d+03, 2.61292842d+01, 1.43447161d+00, 4.92767426d+00, 2.10654699d+02 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_x(i) = tmp * expo_gauss_j_mu_x(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 6) then
|
||||
|
||||
coef_gauss_j_mu_x = (/ -0.08783664d0 , -0.16088711d0 , -0.18464486d0 , -0.0368509d0 , -0.08130028d0 , -0.0126972d0 /)
|
||||
expo_gauss_j_mu_x = (/ 4.09729729d+01, 7.11620618d+00, 2.03692338d+00, 4.10831731d+02, 1.12480198d+00, 1.00000000d+04 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_x(i) = tmp * expo_gauss_j_mu_x(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 7) then
|
||||
|
||||
coef_gauss_j_mu_x = (/ -0.01756495d0 , -0.01023623d0 , -0.06548959d0 , -0.03539446d0 , -0.17150646d0 , -0.15071096d0 , -0.11326834d0 /)
|
||||
expo_gauss_j_mu_x = (/ 9.88572565d+02, 1.21363371d+04, 3.69794870d+01, 1.67364529d+02, 3.03962934d+00, 1.27854005d+00, 9.76383343d+00 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_x(i) = tmp * expo_gauss_j_mu_x(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 8) then
|
||||
|
||||
coef_gauss_j_mu_x = (/ -0.11489205d0 , -0.16008968d0 , -0.12892456d0 , -0.04250838d0 , -0.0718451d0 , -0.02394051d0 , -0.00913353d0 , -0.01285182d0 /)
|
||||
expo_gauss_j_mu_x = (/ 6.97632442d+00, 2.56010878d+00, 1.22760977d+00, 7.47697124d+01, 2.16104215d+01, 2.96549728d+02, 1.40773328d+04, 1.43335159d+03 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_x(i) = tmp * expo_gauss_j_mu_x(i)
|
||||
enddo
|
||||
|
||||
!elseif(ng_fit_jast .eq. 9) then
|
||||
|
||||
! coef_gauss_j_mu_x = (/ /)
|
||||
! expo_gauss_j_mu_x = (/ /)
|
||||
|
||||
! tmp = mu_erf * mu_erf
|
||||
! do i = 1, ng_fit_jast
|
||||
! expo_gauss_j_mu_x(i) = tmp * expo_gauss_j_mu_x(i)
|
||||
! enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 20) then
|
||||
|
||||
ASSERT(n_max_fit_slat == 20)
|
||||
|
||||
alpha = expo_j_xmu(1) * mu_erf
|
||||
call expo_fit_slater_gam(alpha, expos)
|
||||
beta = expo_j_xmu(2) * mu_erf * mu_erf
|
||||
|
||||
tmp = -1.0d0 / sqrt(dacos(-1.d0))
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_x(i) = expos(i) + beta
|
||||
coef_gauss_j_mu_x(i) = tmp * coef_fit_slat_gauss(i)
|
||||
enddo
|
||||
|
||||
else
|
||||
|
||||
print *, ' not implemented yet'
|
||||
stop
|
||||
|
||||
endif
|
||||
|
||||
tmp = 0.5d0 / mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
coef_gauss_j_mu_x(i) = tmp * coef_gauss_j_mu_x(i)
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [double precision, expo_gauss_j_mu_x_2, (ng_fit_jast)]
|
||||
&BEGIN_PROVIDER [double precision, coef_gauss_j_mu_x_2, (ng_fit_jast)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! J(mu,r12)^2 = 0.25/mu^2 F(r12*mu)^2
|
||||
!
|
||||
! F(x)^2 = 1/pi * exp(-2 * alpha * x) exp(-2 * beta * x^2)
|
||||
!
|
||||
! The slater function exp(-2 * alpha * x) is fitted with n_max_fit_slat gaussians
|
||||
!
|
||||
! See Appendix 2 of JCP 154, 084119 (2021)
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i
|
||||
double precision :: tmp
|
||||
double precision :: expos(ng_fit_jast), alpha, beta
|
||||
double precision :: alpha_opt, beta_opt
|
||||
|
||||
if(ng_fit_jast .eq. 1) then
|
||||
|
||||
coef_gauss_j_mu_x_2 = (/ 0.26699573d0 /)
|
||||
expo_gauss_j_mu_x_2 = (/ 11.71029824d0 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_x_2(i) = tmp * expo_gauss_j_mu_x_2(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 2) then
|
||||
|
||||
coef_gauss_j_mu_x_2 = (/ 0.11627934d0 , 0.18708824d0 /)
|
||||
expo_gauss_j_mu_x_2 = (/ 102.41386863d0, 6.36239771d0 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_x_2(i) = tmp * expo_gauss_j_mu_x_2(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 3) then
|
||||
|
||||
coef_gauss_j_mu_x_2 = (/ 0.04947216d0 , 0.14116238d0, 0.12276501d0 /)
|
||||
expo_gauss_j_mu_x_2 = (/ 635.29701766d0, 4.87696954d0, 33.36745891d0 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_x_2(i) = tmp * expo_gauss_j_mu_x_2(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 5) then
|
||||
|
||||
coef_gauss_j_mu_x_2 = (/ 0.01461527d0 , 0.03257147d0 , 0.08831354d0 , 0.11411794d0 , 0.06858783d0 /)
|
||||
expo_gauss_j_mu_x_2 = (/ 8.76554470d+03, 4.90224577d+02, 3.68267125d+00, 1.29663940d+01, 6.58240931d+01 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_x_2(i) = tmp * expo_gauss_j_mu_x_2(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 6) then
|
||||
|
||||
coef_gauss_j_mu_x_2 = (/ 0.01347632d0 , 0.03929124d0 , 0.06289468d0 , 0.10702493d0 , 0.06999865d0 , 0.02558191d0 /)
|
||||
expo_gauss_j_mu_x_2 = (/ 1.00000000d+04, 1.20900717d+02, 3.20346191d+00, 8.92157196d+00, 3.28119120d+01, 6.49045808d+02 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_x_2(i) = tmp * expo_gauss_j_mu_x_2(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 7) then
|
||||
|
||||
coef_gauss_j_mu_x_2 = (/ 0.05202849d0 , 0.01031081d0 , 0.04699157d0 , 0.01451002d0 , 0.07442576d0 , 0.02692033d0 , 0.09311842d0 /)
|
||||
expo_gauss_j_mu_x_2 = (/ 3.04469415d+00, 1.40682034d+04, 7.45960945d+01, 1.43067466d+03, 2.16815661d+01, 2.95750306d+02, 7.23471236d+00 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_x_2(i) = tmp * expo_gauss_j_mu_x_2(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 8) then
|
||||
|
||||
coef_gauss_j_mu_x_2 = (/ 0.00942115d0 , 0.07332421d0 , 0.0508308d0 , 0.08204949d0 , 0.0404099d0 , 0.03201288d0 , 0.01911313d0 , 0.01114732d0 /)
|
||||
expo_gauss_j_mu_x_2 = (/ 1.56957321d+04, 1.52867810d+01, 4.36016903d+01, 5.96818956d+00, 2.85535269d+00, 1.36064008d+02, 4.71968910d+02, 1.92022350d+03 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_x_2(i) = tmp * expo_gauss_j_mu_x_2(i)
|
||||
enddo
|
||||
|
||||
!elseif(ng_fit_jast .eq. 9) then
|
||||
|
||||
! coef_gauss_j_mu_x_2 = (/ /)
|
||||
! expo_gauss_j_mu_x_2 = (/ /)
|
||||
!
|
||||
! tmp = mu_erf * mu_erf
|
||||
! do i = 1, ng_fit_jast
|
||||
! expo_gauss_j_mu_x_2(i) = tmp * expo_gauss_j_mu_x_2(i)
|
||||
! enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 20) then
|
||||
|
||||
ASSERT(n_max_fit_slat == 20)
|
||||
|
||||
!alpha_opt = 2.d0 * expo_j_xmu(1)
|
||||
!beta_opt = 2.d0 * expo_j_xmu(2)
|
||||
|
||||
! direct opt
|
||||
alpha_opt = 3.52751759d0
|
||||
beta_opt = 1.26214809d0
|
||||
|
||||
alpha = alpha_opt * mu_erf
|
||||
call expo_fit_slater_gam(alpha, expos)
|
||||
beta = beta_opt * mu_erf * mu_erf
|
||||
|
||||
tmp = 1.d0 / dacos(-1.d0)
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_x_2(i) = expos(i) + beta
|
||||
coef_gauss_j_mu_x_2(i) = tmp * coef_fit_slat_gauss(i)
|
||||
enddo
|
||||
|
||||
else
|
||||
|
||||
print *, ' not implemented yet'
|
||||
stop
|
||||
|
||||
endif
|
||||
|
||||
tmp = 0.25d0 / (mu_erf * mu_erf)
|
||||
do i = 1, ng_fit_jast
|
||||
coef_gauss_j_mu_x_2(i) = tmp * coef_gauss_j_mu_x_2(i)
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [double precision, expo_gauss_j_mu_1_erf, (ng_fit_jast)]
|
||||
&BEGIN_PROVIDER [double precision, coef_gauss_j_mu_1_erf, (ng_fit_jast)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! J(mu,r12) x \frac{1 - erf(mu * r12)}{2} =
|
||||
!
|
||||
! - \frac{1}{4 \sqrt{\pi} \mu} \exp(-(alpha1 + alpha2) * mu * r12 - (beta1 + beta2) * mu^2 * r12^2)
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i
|
||||
double precision :: tmp
|
||||
double precision :: expos(ng_fit_jast), alpha, beta
|
||||
double precision :: alpha_opt, beta_opt
|
||||
|
||||
if(ng_fit_jast .eq. 1) then
|
||||
|
||||
coef_gauss_j_mu_1_erf = (/ -0.47742461d0 /)
|
||||
expo_gauss_j_mu_1_erf = (/ 8.72255696d0 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_1_erf(i) = tmp * expo_gauss_j_mu_1_erf(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 2) then
|
||||
|
||||
coef_gauss_j_mu_1_erf = (/ -0.19342649d0, -0.34563835d0 /)
|
||||
expo_gauss_j_mu_1_erf = (/ 78.66099999d0, 5.04324363d0 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_1_erf(i) = tmp * expo_gauss_j_mu_1_erf(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 3) then
|
||||
|
||||
coef_gauss_j_mu_1_erf = (/ -0.0802541d0 , -0.27019258d0, -0.20546681d0 /)
|
||||
expo_gauss_j_mu_1_erf = (/ 504.53350764d0, 4.01408169d0, 26.5758329d0 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_1_erf(i) = tmp * expo_gauss_j_mu_1_erf(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 5) then
|
||||
|
||||
coef_gauss_j_mu_1_erf = (/ -0.02330531d0 , -0.11888176d0 , -0.16476192d0 , -0.19874713d0 , -0.05889174d0 /)
|
||||
expo_gauss_j_mu_1_erf = (/ 1.00000000d+04, 4.66067922d+01, 3.04359857d+00, 9.54726649d+00, 3.59796835d+02 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_1_erf(i) = tmp * expo_gauss_j_mu_1_erf(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 6) then
|
||||
|
||||
coef_gauss_j_mu_1_erf = (/ -0.01865654d0 , -0.18319251d0 , -0.06543196d0 , -0.11522778d0 , -0.14825793d0 , -0.03327101d0 /)
|
||||
expo_gauss_j_mu_1_erf = (/ 1.00000000d+04, 8.05593848d+00, 1.27986190d+02, 2.92674319d+01, 2.93583623d+00, 7.65609148d+02 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_1_erf(i) = tmp * expo_gauss_j_mu_1_erf(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 7) then
|
||||
|
||||
coef_gauss_j_mu_1_erf = (/ -0.11853067d0 , -0.01522824d0 , -0.07419098d0 , -0.022202d0 , -0.12242283d0 , -0.04177571d0 , -0.16983107d0 /)
|
||||
expo_gauss_j_mu_1_erf = (/ 2.74057056d+00, 1.37626591d+04, 6.65578663d+01, 1.34693031d+03, 1.90547699d+01, 2.69445390d+02, 6.31845879d+00/)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_j_mu_1_erf(i) = tmp * expo_gauss_j_mu_1_erf(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 8) |