mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-08-24 21:21:47 +02:00
new keywords for Jastrow
This commit is contained in:
parent
ef60141fbf
commit
b4ba0eda6f
@ -1,94 +0,0 @@
|
||||
BEGIN_PROVIDER [integer, n_max_fit_slat]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! number of gaussian to fit exp(-x)
|
||||
!
|
||||
! I took 20 gaussians from the program bassto.f
|
||||
END_DOC
|
||||
n_max_fit_slat = 20
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, coef_fit_slat_gauss, (n_max_fit_slat)]
|
||||
&BEGIN_PROVIDER [double precision, expo_fit_slat_gauss, (n_max_fit_slat)]
|
||||
implicit none
|
||||
include 'constants.include.F'
|
||||
BEGIN_DOC
|
||||
! fit the exp(-x) as
|
||||
!
|
||||
! \sum_{i = 1, n_max_fit_slat} coef_fit_slat_gauss(i) * exp(-expo_fit_slat_gauss(i) * x**2)
|
||||
!
|
||||
! The coefficient are taken from the program bassto.f
|
||||
END_DOC
|
||||
|
||||
|
||||
expo_fit_slat_gauss(01)=30573.77073000000
|
||||
coef_fit_slat_gauss(01)=0.00338925525
|
||||
expo_fit_slat_gauss(02)=5608.45238100000
|
||||
coef_fit_slat_gauss(02)=0.00536433869
|
||||
expo_fit_slat_gauss(03)=1570.95673400000
|
||||
coef_fit_slat_gauss(03)=0.00818702846
|
||||
expo_fit_slat_gauss(04)=541.39785110000
|
||||
coef_fit_slat_gauss(04)=0.01202047655
|
||||
expo_fit_slat_gauss(05)=212.43469630000
|
||||
coef_fit_slat_gauss(05)=0.01711289568
|
||||
expo_fit_slat_gauss(06)=91.31444574000
|
||||
coef_fit_slat_gauss(06)=0.02376001022
|
||||
expo_fit_slat_gauss(07)=42.04087246000
|
||||
coef_fit_slat_gauss(07)=0.03229121736
|
||||
expo_fit_slat_gauss(08)=20.43200443000
|
||||
coef_fit_slat_gauss(08)=0.04303646818
|
||||
expo_fit_slat_gauss(09)=10.37775161000
|
||||
coef_fit_slat_gauss(09)=0.05624657578
|
||||
expo_fit_slat_gauss(10)=5.46880754500
|
||||
coef_fit_slat_gauss(10)=0.07192311571
|
||||
expo_fit_slat_gauss(11)=2.97373529200
|
||||
coef_fit_slat_gauss(11)=0.08949389001
|
||||
expo_fit_slat_gauss(12)=1.66144190200
|
||||
coef_fit_slat_gauss(12)=0.10727599240
|
||||
expo_fit_slat_gauss(13)=0.95052560820
|
||||
coef_fit_slat_gauss(13)=0.12178961750
|
||||
expo_fit_slat_gauss(14)=0.55528683970
|
||||
coef_fit_slat_gauss(14)=0.12740141870
|
||||
expo_fit_slat_gauss(15)=0.33043360020
|
||||
coef_fit_slat_gauss(15)=0.11759168160
|
||||
expo_fit_slat_gauss(16)=0.19982303230
|
||||
coef_fit_slat_gauss(16)=0.08953504394
|
||||
expo_fit_slat_gauss(17)=0.12246840760
|
||||
coef_fit_slat_gauss(17)=0.05066721317
|
||||
expo_fit_slat_gauss(18)=0.07575825322
|
||||
coef_fit_slat_gauss(18)=0.01806363869
|
||||
expo_fit_slat_gauss(19)=0.04690146243
|
||||
coef_fit_slat_gauss(19)=0.00305632563
|
||||
expo_fit_slat_gauss(20)=0.02834749861
|
||||
coef_fit_slat_gauss(20)=0.00013317513
|
||||
|
||||
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
double precision function slater_fit_gam(x,gam)
|
||||
implicit none
|
||||
double precision, intent(in) :: x,gam
|
||||
BEGIN_DOC
|
||||
! fit of the function exp(-gam * x) with gaussian functions
|
||||
END_DOC
|
||||
integer :: i
|
||||
slater_fit_gam = 0.d0
|
||||
do i = 1, n_max_fit_slat
|
||||
slater_fit_gam += coef_fit_slat_gauss(i) * dexp(-expo_fit_slat_gauss(i) * gam * gam * x * x)
|
||||
enddo
|
||||
end
|
||||
|
||||
subroutine expo_fit_slater_gam(gam,expos)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! returns the array of the exponents of the gaussians to fit exp(-gam*x)
|
||||
END_DOC
|
||||
double precision, intent(in) :: gam
|
||||
double precision, intent(out) :: expos(n_max_fit_slat)
|
||||
integer :: i
|
||||
do i = 1, n_max_fit_slat
|
||||
expos(i) = expo_fit_slat_gauss(i) * gam * gam
|
||||
enddo
|
||||
end
|
||||
|
@ -1,335 +0,0 @@
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [integer, n_gauss_eff_pot]
|
||||
|
||||
BEGIN_DOC
|
||||
! number of gaussians to represent the effective potential :
|
||||
!
|
||||
! V(mu,r12) = -0.25 * (1 - erf(mu*r12))^2 + 1/(\sqrt(pi)mu) * exp(-(mu*r12)^2)
|
||||
!
|
||||
! Here (1 - erf(mu*r12))^2 is expanded in Gaussians as Eqs A11-A20 in JCP 154, 084119 (2021)
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
n_gauss_eff_pot = ng_fit_jast + 1
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [integer, n_gauss_eff_pot_deriv]
|
||||
|
||||
BEGIN_DOC
|
||||
! V(r12) = -(1 - erf(mu*r12))^2 is expanded in Gaussians as Eqs A11-A20 in JCP 154, 084119 (2021)
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
n_gauss_eff_pot_deriv = ng_fit_jast
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [double precision, expo_gauss_eff_pot, (n_gauss_eff_pot)]
|
||||
&BEGIN_PROVIDER [double precision, coef_gauss_eff_pot, (n_gauss_eff_pot)]
|
||||
|
||||
BEGIN_DOC
|
||||
! Coefficients and exponents of the Fit on Gaussians of V(X) = -(1 - erf(mu*X))^2 + 1/(\sqrt(pi)mu) * exp(-(mu*X)^2)
|
||||
!
|
||||
! V(X) = \sum_{i=1,n_gauss_eff_pot} coef_gauss_eff_pot(i) * exp(-expo_gauss_eff_pot(i) * X^2)
|
||||
!
|
||||
! Relies on the fit proposed in Eqs A11-A20 in JCP 154, 084119 (2021)
|
||||
END_DOC
|
||||
|
||||
include 'constants.include.F'
|
||||
|
||||
implicit none
|
||||
integer :: i
|
||||
|
||||
! fit of the -0.25 * (1 - erf(mu*x))^2 with n_max_fit_slat gaussians
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_eff_pot(i) = expo_gauss_1_erf_x_2(i)
|
||||
coef_gauss_eff_pot(i) = -0.25d0 * coef_gauss_1_erf_x_2(i) ! -1/4 * (1 - erf(mu*x))^2
|
||||
enddo
|
||||
|
||||
! Analytical Gaussian part of the potential: + 1/(\sqrt(pi)mu) * exp(-(mu*x)^2)
|
||||
expo_gauss_eff_pot(ng_fit_jast+1) = mu_erf * mu_erf
|
||||
coef_gauss_eff_pot(ng_fit_jast+1) = 1.d0 * mu_erf * inv_sq_pi
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
double precision function eff_pot_gauss(x, mu)
|
||||
|
||||
BEGIN_DOC
|
||||
! V(mu,r12) = -0.25 * (1 - erf(mu*r12))^2 + 1/(\sqrt(pi)mu) * exp(-(mu*r12)^2)
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
double precision, intent(in) :: x, mu
|
||||
|
||||
eff_pot_gauss = mu/dsqrt(dacos(-1.d0)) * dexp(-mu*mu*x*x) - 0.25d0 * (1.d0 - derf(mu*x))**2.d0
|
||||
|
||||
end
|
||||
|
||||
! -------------------------------------------------------------------------------------------------
|
||||
! ---
|
||||
|
||||
double precision function eff_pot_fit_gauss(x)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! V(mu,r12) = -0.25 * (1 - erf(mu*r12))^2 + 1/(\sqrt(pi)mu) * exp(-(mu*r12)^2)
|
||||
!
|
||||
! but fitted with gaussians
|
||||
END_DOC
|
||||
double precision, intent(in) :: x
|
||||
integer :: i
|
||||
double precision :: alpha
|
||||
eff_pot_fit_gauss = derf(mu_erf*x)/x
|
||||
do i = 1, n_gauss_eff_pot
|
||||
alpha = expo_gauss_eff_pot(i)
|
||||
eff_pot_fit_gauss += coef_gauss_eff_pot(i) * dexp(-alpha*x*x)
|
||||
enddo
|
||||
end
|
||||
|
||||
BEGIN_PROVIDER [integer, n_fit_1_erf_x]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
!
|
||||
END_DOC
|
||||
n_fit_1_erf_x = 2
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, expos_slat_gauss_1_erf_x, (n_fit_1_erf_x)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! 1 - erf(mu*x) is fitted with a Slater and gaussian as in Eq.A15 of JCP 154, 084119 (2021)
|
||||
!
|
||||
! 1 - erf(mu*x) = e^{-expos_slat_gauss_1_erf_x(1) * mu *x} * e^{-expos_slat_gauss_1_erf_x(2) * mu^2 * x^2}
|
||||
END_DOC
|
||||
expos_slat_gauss_1_erf_x(1) = 1.09529d0
|
||||
expos_slat_gauss_1_erf_x(2) = 0.756023d0
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [double precision, expo_gauss_1_erf_x, (n_max_fit_slat)]
|
||||
&BEGIN_PROVIDER [double precision, coef_gauss_1_erf_x, (n_max_fit_slat)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! (1 - erf(mu*x)) = \sum_i coef_gauss_1_erf_x(i) * exp(-expo_gauss_1_erf_x(i) * x^2)
|
||||
!
|
||||
! This is based on a fit of (1 - erf(mu*x)) by exp(-alpha * x) exp(-beta*mu^2x^2)
|
||||
!
|
||||
! and 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 :: expos(n_max_fit_slat), alpha, beta
|
||||
|
||||
alpha = expos_slat_gauss_1_erf_x(1) * mu_erf
|
||||
call expo_fit_slater_gam(alpha, expos)
|
||||
beta = expos_slat_gauss_1_erf_x(2) * mu_erf * mu_erf
|
||||
|
||||
do i = 1, n_max_fit_slat
|
||||
expo_gauss_1_erf_x(i) = expos(i) + beta
|
||||
coef_gauss_1_erf_x(i) = coef_fit_slat_gauss(i)
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
double precision function fit_1_erf_x(x)
|
||||
|
||||
BEGIN_DOC
|
||||
! fit_1_erf_x(x) = \sum_i c_i exp (-alpha_i x^2) \approx (1 - erf(mu*x))
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i
|
||||
double precision, intent(in) :: x
|
||||
|
||||
fit_1_erf_x = 0.d0
|
||||
do i = 1, n_max_fit_slat
|
||||
fit_1_erf_x += dexp(-expo_gauss_1_erf_x(i) *x*x) * coef_gauss_1_erf_x(i)
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [double precision, expo_gauss_1_erf_x_2, (ng_fit_jast)]
|
||||
&BEGIN_PROVIDER [double precision, coef_gauss_1_erf_x_2, (ng_fit_jast)]
|
||||
|
||||
BEGIN_DOC
|
||||
! (1 - erf(mu*x))^2 = \sum_i coef_gauss_1_erf_x_2(i) * exp(-expo_gauss_1_erf_x_2(i) * x^2)
|
||||
!
|
||||
! This is based on a fit of (1 - erf(mu*x)) by exp(-alpha * x) exp(-beta*mu^2x^2)
|
||||
!
|
||||
! and the slater function exp(-alpha * x) is fitted with n_max_fit_slat gaussians
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i
|
||||
double precision :: expos(ng_fit_jast), alpha, beta, tmp
|
||||
|
||||
if(ng_fit_jast .eq. 1) then
|
||||
|
||||
coef_gauss_1_erf_x_2 = (/ 0.85345277d0 /)
|
||||
expo_gauss_1_erf_x_2 = (/ 6.23519457d0 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_1_erf_x_2(i) = tmp * expo_gauss_1_erf_x_2(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 2) then
|
||||
|
||||
coef_gauss_1_erf_x_2 = (/ 0.31030624d0 , 0.64364964d0 /)
|
||||
expo_gauss_1_erf_x_2 = (/ 55.39184787d0, 3.92151407d0 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_1_erf_x_2(i) = tmp * expo_gauss_1_erf_x_2(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 3) then
|
||||
|
||||
coef_gauss_1_erf_x_2 = (/ 0.33206082d0 , 0.52347449d0, 0.12605012d0 /)
|
||||
expo_gauss_1_erf_x_2 = (/ 19.90272209d0, 3.2671671d0 , 336.47320445d0 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_1_erf_x_2(i) = tmp * expo_gauss_1_erf_x_2(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 5) then
|
||||
|
||||
coef_gauss_1_erf_x_2 = (/ 0.02956716d0, 0.17025555d0, 0.32774114d0, 0.39034764d0, 0.07822781d0 /)
|
||||
expo_gauss_1_erf_x_2 = (/ 6467.28126d0, 46.9071990d0, 9.09617721d0, 2.76883328d0, 360.367093d0 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_1_erf_x_2(i) = tmp * expo_gauss_1_erf_x_2(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 6) then
|
||||
|
||||
coef_gauss_1_erf_x_2 = (/ 0.18331042d0 , 0.10971118d0 , 0.29949169d0 , 0.34853132d0 , 0.0394275d0 , 0.01874444d0 /)
|
||||
expo_gauss_1_erf_x_2 = (/ 2.54293498d+01, 1.40317872d+02, 7.14630801d+00, 2.65517675d+00, 1.45142619d+03, 1.00000000d+04 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_1_erf_x_2(i) = tmp * expo_gauss_1_erf_x_2(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 7) then
|
||||
|
||||
coef_gauss_1_erf_x_2 = (/ 0.0213619d0 , 0.03221511d0 , 0.29966689d0 , 0.19178934d0 , 0.06154732d0 , 0.28214555d0 , 0.11125985d0 /)
|
||||
expo_gauss_1_erf_x_2 = (/ 1.34727067d+04, 1.27166613d+03, 5.52584567d+00, 1.67753218d+01, 2.46145691d+02, 2.47971820d+00, 5.95141293d+01 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_1_erf_x_2(i) = tmp * expo_gauss_1_erf_x_2(i)
|
||||
enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 8) then
|
||||
|
||||
coef_gauss_1_erf_x_2 = (/ 0.28189124d0 , 0.19518669d0 , 0.12161735d0 , 0.24257438d0 , 0.07309656d0 , 0.042435d0 , 0.01926109d0 , 0.02393415d0 /)
|
||||
expo_gauss_1_erf_x_2 = (/ 4.69795903d+00, 1.21379451d+01, 3.55527053d+01, 2.39227172d+00, 1.14827721d+02, 4.16320213d+02, 1.52813587d+04, 1.78516557d+03 /)
|
||||
|
||||
tmp = mu_erf * mu_erf
|
||||
do i = 1, ng_fit_jast
|
||||
expo_gauss_1_erf_x_2(i) = tmp * expo_gauss_1_erf_x_2(i)
|
||||
enddo
|
||||
|
||||
!elseif(ng_fit_jast .eq. 9) then
|
||||
|
||||
! coef_gauss_1_erf_x_2 = (/ /)
|
||||
! expo_gauss_1_erf_x_2 = (/ /)
|
||||
|
||||
! tmp = mu_erf * mu_erf
|
||||
! do i = 1, ng_fit_jast
|
||||
! expo_gauss_1_erf_x_2(i) = tmp * expo_gauss_1_erf_x_2(i)
|
||||
! enddo
|
||||
|
||||
elseif(ng_fit_jast .eq. 20) then
|
||||
|
||||
ASSERT(n_max_fit_slat == 20)
|
||||
|
||||
alpha = 2.d0 * expos_slat_gauss_1_erf_x(1) * mu_erf
|
||||
call expo_fit_slater_gam(alpha, expos)
|
||||
beta = 2.d0 * expos_slat_gauss_1_erf_x(2) * mu_erf * mu_erf
|
||||
do i = 1, n_max_fit_slat
|
||||
expo_gauss_1_erf_x_2(i) = expos(i) + beta
|
||||
coef_gauss_1_erf_x_2(i) = coef_fit_slat_gauss(i)
|
||||
enddo
|
||||
|
||||
else
|
||||
|
||||
print *, ' not implemented yet'
|
||||
stop
|
||||
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
double precision function fit_1_erf_x_2(x)
|
||||
implicit none
|
||||
double precision, intent(in) :: x
|
||||
BEGIN_DOC
|
||||
! fit_1_erf_x_2(x) = \sum_i c_i exp (-alpha_i x^2) \approx (1 - erf(mu*x))^2
|
||||
END_DOC
|
||||
integer :: i
|
||||
fit_1_erf_x_2 = 0.d0
|
||||
do i = 1, n_max_fit_slat
|
||||
fit_1_erf_x_2 += dexp(-expo_gauss_1_erf_x_2(i) *x*x) * coef_gauss_1_erf_x_2(i)
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
subroutine inv_r_times_poly(r, dist_r, dist_vec, poly)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! returns
|
||||
!
|
||||
! poly(1) = x / sqrt(x^2+y^2+z^2), poly(2) = y / sqrt(x^2+y^2+z^2), poly(3) = z / sqrt(x^2+y^2+z^2)
|
||||
!
|
||||
! with the arguments
|
||||
!
|
||||
! r(1) = x, r(2) = y, r(3) = z, dist_r = sqrt(x^2+y^2+z^2)
|
||||
!
|
||||
! dist_vec(1) = sqrt(y^2+z^2), dist_vec(2) = sqrt(x^2+z^2), dist_vec(3) = sqrt(x^2+y^2)
|
||||
END_DOC
|
||||
double precision, intent(in) :: r(3), dist_r, dist_vec(3)
|
||||
double precision, intent(out):: poly(3)
|
||||
double precision :: inv_dist
|
||||
integer :: i
|
||||
if (dist_r.gt. 1.d-8)then
|
||||
inv_dist = 1.d0/dist_r
|
||||
do i = 1, 3
|
||||
poly(i) = r(i) * inv_dist
|
||||
enddo
|
||||
else
|
||||
do i = 1, 3
|
||||
if(dabs(r(i)).lt.dist_vec(i))then
|
||||
inv_dist = 1.d0/dist_r
|
||||
poly(i) = r(i) * inv_dist
|
||||
else !if(dabs(r(i)))then
|
||||
poly(i) = 1.d0
|
||||
! poly(i) = 0.d0
|
||||
endif
|
||||
enddo
|
||||
endif
|
||||
end
|
@ -174,7 +174,7 @@ double precision function general_primitive_integral_coul_shifted( dim
|
||||
general_primitive_integral_coul_shifted = fact_p * fact_q * accu * pi_5_2 * p_inv * q_inv / dsqrt(p_plus_q)
|
||||
|
||||
return
|
||||
end function general_primitive_integral_coul_shifted
|
||||
end
|
||||
!______________________________________________________________________________________________________________________
|
||||
!______________________________________________________________________________________________________________________
|
||||
|
||||
@ -354,7 +354,7 @@ double precision function general_primitive_integral_erf_shifted( dim
|
||||
general_primitive_integral_erf_shifted = fact_p * fact_q * accu * pi_5_2 * p_inv * q_inv / dsqrt(p_plus_q)
|
||||
|
||||
return
|
||||
end function general_primitive_integral_erf_shifted
|
||||
end
|
||||
!______________________________________________________________________________________________________________________
|
||||
!______________________________________________________________________________________________________________________
|
||||
|
||||
@ -362,3 +362,48 @@ end function general_primitive_integral_erf_shifted
|
||||
|
||||
|
||||
|
||||
|
||||
! ---
|
||||
|
||||
subroutine inv_r_times_poly(r, dist_r, dist_vec, poly)
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! returns
|
||||
!
|
||||
! poly(1) = x / sqrt(x^2+y^2+z^2), poly(2) = y / sqrt(x^2+y^2+z^2), poly(3) = z / sqrt(x^2+y^2+z^2)
|
||||
!
|
||||
! with the arguments
|
||||
!
|
||||
! r(1) = x, r(2) = y, r(3) = z, dist_r = sqrt(x^2+y^2+z^2)
|
||||
!
|
||||
! dist_vec(1) = sqrt(y^2+z^2), dist_vec(2) = sqrt(x^2+z^2), dist_vec(3) = sqrt(x^2+y^2)
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
double precision, intent(in) :: r(3), dist_r, dist_vec(3)
|
||||
double precision, intent(out) :: poly(3)
|
||||
integer :: i
|
||||
double precision :: inv_dist
|
||||
|
||||
if (dist_r .gt. 1.d-8)then
|
||||
inv_dist = 1.d0/dist_r
|
||||
do i = 1, 3
|
||||
poly(i) = r(i) * inv_dist
|
||||
enddo
|
||||
else
|
||||
do i = 1, 3
|
||||
if(dabs(r(i)).lt.dist_vec(i)) then
|
||||
inv_dist = 1.d0/dist_r
|
||||
poly(i) = r(i) * inv_dist
|
||||
else
|
||||
poly(i) = 1.d0
|
||||
endif
|
||||
enddo
|
||||
endif
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
|
@ -161,7 +161,7 @@ double precision function env_nucl(r)
|
||||
|
||||
else
|
||||
|
||||
print *, ' Error in grad1_env_nucl: Unknown env_type = ', env_type
|
||||
print *, ' Error in env_nucl: Unknown env_type = ', env_type
|
||||
stop
|
||||
|
||||
endif
|
||||
@ -230,7 +230,7 @@ double precision function env_nucl_square(r)
|
||||
|
||||
else
|
||||
|
||||
print *, ' Error in grad1_env_nucl: Unknown env_type = ', env_type
|
||||
print *, ' Error in env_nucl_square: Unknown env_type = ', env_type
|
||||
stop
|
||||
|
||||
endif
|
||||
|
@ -7,8 +7,7 @@ subroutine get_grad1_u12_withsq_r1_seq(r1, n_grid2, resx, resy, resz, res)
|
||||
!
|
||||
! grad_1 u(r1,r2)
|
||||
!
|
||||
! this will be integrated numerically over r2:
|
||||
! we use grid for r1 and extra_grid for r2
|
||||
! we use grid for r1 and extra_grid for r2
|
||||
!
|
||||
END_DOC
|
||||
|
||||
@ -29,13 +28,11 @@ subroutine get_grad1_u12_withsq_r1_seq(r1, n_grid2, resx, resy, resz, res)
|
||||
PROVIDE final_grid_points_extra
|
||||
|
||||
if( ((j2e_type .eq. "rs-dft") .and. (env_type .eq. "none")) .or. &
|
||||
(j2e_type .eq. "rs-dft-murho") ) then
|
||||
(j2e_type .eq. "rs-dft-murho") ) then
|
||||
|
||||
call grad1_j12_mu_r1_seq(r1, n_grid2, resx, resy, resz)
|
||||
do jpoint = 1, n_points_extra_final_grid
|
||||
res(jpoint) = resx(jpoint) * resx(jpoint) &
|
||||
+ resy(jpoint) * resy(jpoint) &
|
||||
+ resz(jpoint) * resz(jpoint)
|
||||
res(jpoint) = resx(jpoint) * resx(jpoint) + resy(jpoint) * resy(jpoint) + resz(jpoint) * resz(jpoint)
|
||||
enddo
|
||||
|
||||
elseif((j2e_type .eq. "rs-dft") .and. (env_type .ne. "none")) then
|
||||
@ -60,9 +57,7 @@ subroutine get_grad1_u12_withsq_r1_seq(r1, n_grid2, resx, resy, resz, res)
|
||||
resx(jpoint) = (gradx1_u2b(jpoint) * env_r1 + u2b_r12(jpoint) * grad1_env(1)) * env_r2(jpoint)
|
||||
resy(jpoint) = (grady1_u2b(jpoint) * env_r1 + u2b_r12(jpoint) * grad1_env(2)) * env_r2(jpoint)
|
||||
resz(jpoint) = (gradz1_u2b(jpoint) * env_r1 + u2b_r12(jpoint) * grad1_env(3)) * env_r2(jpoint)
|
||||
res (jpoint) = resx(jpoint) * resx(jpoint) &
|
||||
+ resy(jpoint) * resy(jpoint) &
|
||||
+ resz(jpoint) * resz(jpoint)
|
||||
res (jpoint) = resx(jpoint) * resx(jpoint) + resy(jpoint) * resy(jpoint) + resz(jpoint) * resz(jpoint)
|
||||
enddo
|
||||
|
||||
deallocate(env_r2, u2b_r12, gradx1_u2b, grady1_u2b, gradz1_u2b)
|
||||
|
@ -1,4 +1,6 @@
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [double precision, int2_grad1_u12_ao_num , (ao_num,ao_num,n_points_final_grid,3)]
|
||||
&BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao_num, (ao_num,ao_num,n_points_final_grid) ]
|
||||
|
||||
|
@ -10,6 +10,11 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
|
||||
! ao_two_e_tc_tot(k,i,l,j) = (ki|V^TC(r_12)|lj)
|
||||
! = <lk| V^TC(r_12) |ji> where V^TC(r_12) is the total TC operator
|
||||
! = tc_grad_and_lapl_ao(k,i,l,j) + tc_grad_square_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j)
|
||||
! AND IF(var_tc):
|
||||
!
|
||||
! ao_two_e_tot(k,i,l,j) = (ki|V^TC(r_12) + [(V^TC)(r_12)]^\dagger|lj) / 2.0
|
||||
! = tc_grad_square_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j)
|
||||
!
|
||||
!
|
||||
! where:
|
||||
!
|
||||
@ -25,7 +30,6 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
|
||||
|
||||
implicit none
|
||||
integer :: i, j, k, l, m, ipoint
|
||||
double precision :: wall1, wall0
|
||||
double precision :: weight1, ao_k_r, ao_i_r
|
||||
double precision :: der_envsq_x, der_envsq_y, der_envsq_z, lap_envsq
|
||||
double precision :: time0, time1
|
||||
@ -36,7 +40,7 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
|
||||
PROVIDE j2e_type
|
||||
PROVIDE j1e_type
|
||||
|
||||
call wall_time(wall0)
|
||||
call wall_time(time0)
|
||||
|
||||
print *, ' providing ao_two_e_tc_tot ...'
|
||||
print*, ' j2e_type: ', j2e_type
|
||||
@ -58,44 +62,6 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
|
||||
|
||||
! ---
|
||||
|
||||
PROVIDE int2_grad1_u12_ao
|
||||
|
||||
allocate(b_mat(n_points_final_grid,ao_num,ao_num,3))
|
||||
|
||||
b_mat = 0.d0
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i, k, ipoint, weight1, ao_i_r, ao_k_r) &
|
||||
!$OMP SHARED (aos_in_r_array_transp, aos_grad_in_r_array_transp_bis, b_mat, &
|
||||
!$OMP ao_num, n_points_final_grid, final_weight_at_r_vector)
|
||||
!$OMP DO SCHEDULE (static)
|
||||
do i = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
do ipoint = 1, n_points_final_grid
|
||||
|
||||
weight1 = 0.5d0 * final_weight_at_r_vector(ipoint)
|
||||
ao_i_r = aos_in_r_array_transp(ipoint,i)
|
||||
ao_k_r = aos_in_r_array_transp(ipoint,k)
|
||||
|
||||
b_mat(ipoint,k,i,1) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,1) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1))
|
||||
b_mat(ipoint,k,i,2) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,2) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2))
|
||||
b_mat(ipoint,k,i,3) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,3) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3))
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
ao_two_e_tc_tot = 0.d0
|
||||
do m = 1, 3
|
||||
call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
|
||||
, int2_grad1_u12_ao(1,1,1,m), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid &
|
||||
, 1.d0, ao_two_e_tc_tot, ao_num*ao_num)
|
||||
enddo
|
||||
deallocate(b_mat)
|
||||
|
||||
! ---
|
||||
|
||||
PROVIDE int2_grad1_u12_square_ao
|
||||
|
||||
allocate(c_mat(n_points_final_grid,ao_num,ao_num))
|
||||
@ -122,12 +88,11 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
|
||||
|
||||
FREE int2_grad1_u12_square_ao
|
||||
|
||||
if( (j2e_type .eq. "rs-dft") .and. &
|
||||
if( (tc_integ_type .eq. "semi-analytic") .and. &
|
||||
(j2e_type .eq. "rs-dft") .and. &
|
||||
((env_type .eq. "prod_gauss") .or. (env_type .eq. "sum-gauss")) .and. &
|
||||
use_ipp ) then
|
||||
|
||||
print*, " going through Manu's IPP"
|
||||
|
||||
! an additional term is added here directly instead of
|
||||
! being added in int2_grad1_u12_square_ao for performance
|
||||
|
||||
@ -170,6 +135,47 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
|
||||
|
||||
! ---
|
||||
|
||||
if(.not. var_tc) then
|
||||
|
||||
PROVIDE int2_grad1_u12_ao
|
||||
|
||||
allocate(b_mat(n_points_final_grid,ao_num,ao_num,3))
|
||||
|
||||
b_mat = 0.d0
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i, k, ipoint, weight1, ao_i_r, ao_k_r) &
|
||||
!$OMP SHARED (aos_in_r_array_transp, aos_grad_in_r_array_transp_bis, b_mat, &
|
||||
!$OMP ao_num, n_points_final_grid, final_weight_at_r_vector)
|
||||
!$OMP DO SCHEDULE (static)
|
||||
do i = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
do ipoint = 1, n_points_final_grid
|
||||
|
||||
weight1 = 0.5d0 * final_weight_at_r_vector(ipoint)
|
||||
ao_i_r = aos_in_r_array_transp(ipoint,i)
|
||||
ao_k_r = aos_in_r_array_transp(ipoint,k)
|
||||
|
||||
b_mat(ipoint,k,i,1) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,1) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1))
|
||||
b_mat(ipoint,k,i,2) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,2) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2))
|
||||
b_mat(ipoint,k,i,3) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,3) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3))
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
do m = 1, 3
|
||||
call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
|
||||
, int2_grad1_u12_ao(1,1,1,m), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid &
|
||||
, 1.d0, ao_two_e_tc_tot, ao_num*ao_num)
|
||||
enddo
|
||||
deallocate(b_mat)
|
||||
|
||||
endif ! var_tc
|
||||
|
||||
! ---
|
||||
|
||||
call sum_A_At(ao_two_e_tc_tot(1,1,1,1), ao_num*ao_num)
|
||||
|
||||
PROVIDE ao_integrals_map
|
||||
@ -191,7 +197,7 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
if(tc_integ_type .ge. "numeric") then
|
||||
if(tc_integ_type .eq. "numeric") then
|
||||
FREE int2_grad1_u12_ao_num int2_grad1_u12_square_ao_num
|
||||
endif
|
||||
|
||||
@ -214,172 +220,3 @@ END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [double precision, ao_two_e_vartc_tot, (ao_num, ao_num, ao_num, ao_num)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! CHEMIST NOTATION IS USED
|
||||
!
|
||||
! ao_two_e_vartc_tot(k,i,l,j) = (ki|V^TC(r_12)|lj)
|
||||
! = <lk| V^TC(r_12) |ji> where V^TC(r_12) is the total TC operator
|
||||
! = tc_grad_square_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j)
|
||||
!
|
||||
! where:
|
||||
!
|
||||
! tc_grad_square_ao(k,i,l,j) = -1/2 <kl | |\grad_1 u(r1,r2)|^2 + |\grad_2 u(r1,r2)|^2 | ij>
|
||||
!
|
||||
! ao_two_e_coul(k,i,l,j) = < l k | 1/r12 | j i > = ( k i | 1/r12 | l j )
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, k, l, ipoint
|
||||
double precision :: wall1, wall0
|
||||
double precision :: weight1, ao_k_r, ao_i_r
|
||||
double precision :: der_envsq_x, der_envsq_y, der_envsq_z, lap_envsq
|
||||
double precision :: time0, time1
|
||||
double precision, allocatable :: c_mat(:,:,:)
|
||||
double precision, external :: get_ao_two_e_integral
|
||||
|
||||
PROVIDE env_type
|
||||
PROVIDE j2e_type
|
||||
PROVIDE j1e_type
|
||||
|
||||
call wall_time(wall0)
|
||||
|
||||
print *, ' providing ao_two_e_vartc_tot ...'
|
||||
print*, ' j2e_type: ', j2e_type
|
||||
print*, ' j1e_type: ', j1e_type
|
||||
print*, ' env_type: ', env_type
|
||||
|
||||
if(read_tc_integ) then
|
||||
|
||||
print*, ' Reading ao_two_e_vartc_tot from ', trim(ezfio_filename) // '/work/ao_two_e_vartc_tot'
|
||||
|
||||
open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/ao_two_e_vartc_tot', action="read")
|
||||
read(11) ao_two_e_vartc_tot
|
||||
close(11)
|
||||
|
||||
else
|
||||
|
||||
PROVIDE tc_integ_type
|
||||
print*, ' approach for integrals: ', tc_integ_type
|
||||
|
||||
PROVIDE int2_grad1_u12_square_ao
|
||||
|
||||
allocate(c_mat(n_points_final_grid,ao_num,ao_num))
|
||||
|
||||
c_mat = 0.d0
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i, k, ipoint) &
|
||||
!$OMP SHARED (aos_in_r_array_transp, c_mat, ao_num, n_points_final_grid, final_weight_at_r_vector)
|
||||
!$OMP DO SCHEDULE (static)
|
||||
do i = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
do ipoint = 1, n_points_final_grid
|
||||
c_mat(ipoint,k,i) = final_weight_at_r_vector(ipoint) * aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,k)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
|
||||
, int2_grad1_u12_square_ao(1,1,1), ao_num*ao_num, c_mat(1,1,1), n_points_final_grid &
|
||||
, 0.d0, ao_two_e_vartc_tot, ao_num*ao_num)
|
||||
|
||||
FREE int2_grad1_u12_square_ao
|
||||
|
||||
if( (j2e_type .eq. "rs-dft") .and. &
|
||||
((env_type .eq. "prod_gauss") .or. (env_type .eq. "sum-gauss")) .and. &
|
||||
use_ipp ) then
|
||||
|
||||
print*, " going through Manu's IPP"
|
||||
|
||||
! an additional term is added here directly instead of
|
||||
! being added in int2_grad1_u12_square_ao for performance
|
||||
|
||||
PROVIDE int2_u2_env2
|
||||
|
||||
c_mat = 0.d0
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i, k, ipoint, weight1, ao_i_r, ao_k_r) &
|
||||
!$OMP SHARED (aos_in_r_array_transp, c_mat, ao_num, n_points_final_grid, final_weight_at_r_vector, &
|
||||
!$OMP env_square_grad, env_square_lapl, aos_grad_in_r_array_transp_bis)
|
||||
!$OMP DO SCHEDULE (static)
|
||||
do i = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
do ipoint = 1, n_points_final_grid
|
||||
|
||||
weight1 = 0.25d0 * final_weight_at_r_vector(ipoint)
|
||||
|
||||
ao_i_r = aos_in_r_array_transp(ipoint,i)
|
||||
ao_k_r = aos_in_r_array_transp(ipoint,k)
|
||||
|
||||
c_mat(ipoint,k,i) = weight1 * ( ao_k_r * ao_i_r * env_square_lapl(ipoint) &
|
||||
+ (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,1) + ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1)) * env_square_grad(ipoint,1) &
|
||||
+ (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,2) + ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2)) * env_square_grad(ipoint,2) &
|
||||
+ (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,3) + ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3)) * env_square_grad(ipoint,3) )
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
|
||||
, int2_u2_env2(1,1,1), ao_num*ao_num, c_mat(1,1,1), n_points_final_grid &
|
||||
, 1.d0, ao_two_e_vartc_tot, ao_num*ao_num)
|
||||
|
||||
FREE int2_u2_env2
|
||||
endif ! use_ipp
|
||||
|
||||
deallocate(c_mat)
|
||||
|
||||
! ---
|
||||
|
||||
call sum_A_At(ao_two_e_vartc_tot(1,1,1,1), ao_num*ao_num)
|
||||
|
||||
PROVIDE ao_integrals_map
|
||||
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP SHARED(ao_num, ao_two_e_vartc_tot, ao_integrals_map) &
|
||||
!$OMP PRIVATE(i, j, k, l)
|
||||
!$OMP DO
|
||||
do j = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
! < 1:i, 2:j | 1:k, 2:l >
|
||||
ao_two_e_vartc_tot(k,i,l,j) = ao_two_e_vartc_tot(k,i,l,j) + get_ao_two_e_integral(i, j, k, l, ao_integrals_map)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
if(tc_integ_type .ge. "numeric") then
|
||||
FREE int2_grad1_u12_ao_num int2_grad1_u12_square_ao_num
|
||||
endif
|
||||
|
||||
endif ! read_tc_integ
|
||||
|
||||
if(write_tc_integ .and. mpi_master) then
|
||||
print*, ' Saving ao_two_e_vartc_tot in ', trim(ezfio_filename) // '/work/ao_two_e_vartc_tot'
|
||||
open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/ao_two_e_vartc_tot', action="write")
|
||||
call ezfio_set_work_empty(.False.)
|
||||
write(11) ao_two_e_vartc_tot
|
||||
close(11)
|
||||
call ezfio_set_tc_keywords_io_tc_integ('Read')
|
||||
endif
|
||||
|
||||
call wall_time(time1)
|
||||
print*, ' Wall time for ao_two_e_vartc_tot (min) = ', (time1 - time0) / 60.d0
|
||||
call print_memory_usage()
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
|
@ -160,12 +160,6 @@ doc: If |true|, maximize the overlap between orthogonalized left- and right eige
|
||||
interface: ezfio,provider,ocaml
|
||||
default: False
|
||||
|
||||
[ng_fit_jast]
|
||||
type: integer
|
||||
doc: nb of Gaussians used to fit Jastrow fcts
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 20
|
||||
|
||||
[max_dim_diis_tcscf]
|
||||
type: integer
|
||||
doc: Maximum size of the DIIS extrapolation procedure
|
||||
@ -258,7 +252,7 @@ default: True
|
||||
|
||||
[tc_grid1_a]
|
||||
type: integer
|
||||
doc: size of angular grid over r1
|
||||
doc: size of angular grid over r1: [ 6 | 14 | 26 | 38 | 50 | 74 | 86 | 110 | 146 | 170 | 194 | 230 | 266 | 302 | 350 | 434 | 590 | 770 | 974 | 1202 | 1454 | 1730 | 2030 | 2354 | 2702 | 3074 | 3470 | 3890 | 4334 | 4802 | 5294 | 5810 ]
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 50
|
||||
|
||||
@ -270,19 +264,19 @@ default: 30
|
||||
|
||||
[tc_grid2_a]
|
||||
type: integer
|
||||
doc: size of angular grid over r2
|
||||
doc: size of angular grid over r2: [ 6 | 14 | 26 | 38 | 50 | 74 | 86 | 110 | 146 | 170 | 194 | 230 | 266 | 302 | 350 | 434 | 590 | 770 | 974 | 1202 | 1454 | 1730 | 2030 | 2354 | 2702 | 3074 | 3470 | 3890 | 4334 | 4802 | 5294 | 5810 ]
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 194
|
||||
default: 266
|
||||
|
||||
[tc_grid2_r]
|
||||
type: integer
|
||||
doc: size of radial grid over r2
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 50
|
||||
default: 70
|
||||
|
||||
[tc_integ_type]
|
||||
type: character*(32)
|
||||
doc: approach used to evaluate TC integrals [analytic | numeric | semi-analytic]
|
||||
doc: approach used to evaluate TC integrals [ analytic | numeric | semi-analytic ]
|
||||
interface: ezfio,ocaml,provider
|
||||
default: semi-analytic
|
||||
|
||||
|
@ -13,9 +13,9 @@
|
||||
two_e_vartc_integral_alpha = 0.d0
|
||||
two_e_vartc_integral_beta = 0.d0
|
||||
|
||||
!$OMP PARALLEL DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i, j, k, l, density_a, density_b, density, tmp_a, tmp_b, I_coul, I_kjli) &
|
||||
!$OMP SHARED (ao_num, TCSCF_density_matrix_ao_alpha, TCSCF_density_matrix_ao_beta, ao_two_e_vartc_tot, &
|
||||
!$OMP PARALLEL DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i, j, k, l, density_a, density_b, density, tmp_a, tmp_b, I_coul, I_kjli) &
|
||||
!$OMP SHARED (ao_num, TCSCF_density_matrix_ao_alpha, TCSCF_density_matrix_ao_beta, ao_two_e_tc_tot, &
|
||||
!$OMP two_e_vartc_integral_alpha, two_e_vartc_integral_beta)
|
||||
|
||||
allocate(tmp_a(ao_num,ao_num), tmp_b(ao_num,ao_num))
|
||||
@ -31,8 +31,8 @@
|
||||
do i = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
|
||||
I_coul = density * ao_two_e_vartc_tot(k,i,l,j)
|
||||
I_kjli = ao_two_e_vartc_tot(k,j,l,i)
|
||||
I_coul = density * ao_two_e_tc_tot(k,i,l,j)
|
||||
I_kjli = ao_two_e_tc_tot(k,j,l,i)
|
||||
|
||||
tmp_a(k,i) += I_coul - density_a * I_kjli
|
||||
tmp_b(k,i) += I_coul - density_b * I_kjli
|
||||
|
@ -45,7 +45,6 @@ program test_ints
|
||||
!!PROVIDE TC_HF_energy VARTC_HF_energy
|
||||
!!print *, ' TC_HF_energy = ', TC_HF_energy
|
||||
!!print *, ' VARTC_HF_energy = ', VARTC_HF_energy
|
||||
! call test_old_ints
|
||||
|
||||
call test_fock_3e_uhf_mo_cs()
|
||||
call test_fock_3e_uhf_mo_a()
|
||||
@ -796,41 +795,6 @@ end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine test_old_ints
|
||||
implicit none
|
||||
integer :: i,j,k,l
|
||||
double precision :: old, new, contrib, get_ao_tc_sym_two_e_pot
|
||||
double precision :: integral_sym , integral_nsym,accu
|
||||
PROVIDE ao_tc_sym_two_e_pot_in_map
|
||||
accu = 0.d0
|
||||
do j = 1, ao_num
|
||||
do l= 1, ao_num
|
||||
do i = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
! integral_sym = get_ao_tc_sym_two_e_pot(i, j, k, l, ao_tc_sym_two_e_pot_map)
|
||||
! ao_non_hermit_term_chemist(k,i,l,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the AO basis
|
||||
! integral_nsym = ao_non_hermit_term_chemist(k,i,l,j)
|
||||
! old = integral_sym + integral_nsym
|
||||
new = ao_tc_int_chemist_test(k,i,l,j)
|
||||
old = ao_tc_int_chemist_no_cycle(k,i,l,j)
|
||||
contrib = dabs(old - new)
|
||||
if(contrib.gt.1.d-6)then
|
||||
print*,'problem !!'
|
||||
print*,i,j,k,l
|
||||
print*,old, new, contrib
|
||||
endif
|
||||
accu += contrib
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
print*,'******'
|
||||
print*,'******'
|
||||
print*,'in test_old_ints'
|
||||
print*,'accu = ',accu/dble(ao_num**4)
|
||||
|
||||
end
|
||||
|
||||
subroutine test_int2_grad1_u12_ao_test
|
||||
implicit none
|
||||
integer :: i,j,ipoint,m,k,l
|
||||
|
Loading…
Reference in New Issue
Block a user