PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-31 00:25:53 +01:00
Code Releases Activity

improv cpx erf

Browse Source
This commit is contained in:
AbdAmmar 2024-10-13 01:33:24 +02:00
parent 23a242052d
commit cd3c0cb1a9
3 changed files with 184 additions and 111 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
src/hartree_fock/deb_ao_2e_int.irp.f
View File

@ -3,10 +3,10 @@ program deb_ao_2e_int
implicit none
call check_ao_two_e_integral_cosgtos()
!call check_ao_two_e_integral_cosgtos()
!call check_crint1()
!call check_crint2()
!call check_crint3()
call check_crint3()
end
@ -26,31 +26,35 @@ subroutine check_ao_two_e_integral_cosgtos()
acc = 0.d0
nrm = 0.d0
i = 11
j = 100
k = 74
l = 104
! do i = 1, ao_num
! do k = 1, ao_num
! do j = 1, ao_num
! do l = 1, ao_num
!i = 11
!j = 100
!k = 74
!l = 104
do i = 1, ao_num
do k = 1, ao_num
j = i
l = k
!do j = 1, ao_num
! do l = 1, ao_num
tmp1 = ao_two_e_integral (i, j, k, l)
tmp2 = ao_two_e_integral_cosgtos(i, j, k, l)
dif = abs(tmp1 - tmp2)
!if(dif .gt. 1d-10) then
if(tmp1 .lt. 0.d0) then
print*, ' error on:', i, j, k, l
print*, tmp1, tmp2, dif
!stop
endif
!endif
acc += dif
nrm += abs(tmp1)
! enddo
! enddo
! enddo
! enddo
enddo
enddo
print *, ' acc (%) = ', dif * 100.d0 / nrm
@ -218,7 +222,7 @@ subroutine check_crint3()
t_int1 = 0.d0
t_int2 = 0.d0
n_test = 5
n_test = 1
acc_re = 0.d0
nrm_re = 0.d0
@ -249,9 +253,14 @@ subroutine check_crint3()
n = int(x)
!if(n.eq.0) cycle
n = 0
!rho = (-6.83897018210218d0, -7.24479852507338d0)
rho = (-9.83206247355480d0, 0.445269582329036d0)
print*, " n = ", n
print*, " rho = ", real(rho), aimag(rho)
call wall_time(t1)
int1_old = crint_2(n, rho)
!n_quad = 10000000

196
src/utils/cpx_boys.irp.f
View File

@ -16,7 +16,6 @@ complex*16 function crint_1(n, rho)
complex*16 :: sq_rho, rho_inv, rho_exp
complex*16 :: crint_smallz, cpx_erf_1
complex*16 :: cpx_erf_2
rho_re = real (rho)
rho_im = aimag(rho)
@ -48,10 +47,6 @@ complex*16 function crint_1(n, rho)
rho_exp = 0.5d0 * zexp(-rho)
rho_inv = (1.d0, 0.d0) / rho
!print*, sq_rho_re, sq_rho_im
!print*, cpx_erf_1(sq_rho_re, sq_rho_im)
!print*, cpx_erf_2(sq_rho_re, sq_rho_im)
crint_1 = 0.5d0 * sqpi * cpx_erf_1(sq_rho_re, sq_rho_im) / sq_rho
mmax = n
if(mmax .gt. 0) then
@ -185,7 +180,7 @@ complex*16 function crint_smallz(n, rho)
complex*16, intent(in) :: rho
integer, parameter :: kmax = 40
double precision, parameter :: eps = 1.d-13
double precision, parameter :: eps = 1.d-10
integer :: k
double precision :: delta_mod
@ -207,7 +202,8 @@ complex*16 function crint_smallz(n, rho)
if(delta_mod > eps) then
write(*,*) ' pb in crint_smallz !'
write(*,*) ' n, rho = ', n, rho
write(*,*) ' n, rho = ', n, rho
write(*,*) ' value = ', crint_smallz
write(*,*) ' delta_mod = ', delta_mod
!stop 1
endif
@ -229,29 +225,93 @@ complex*16 function crint_2(n, rho)
complex*16, external :: crint_smallz
if(abs(rho) < 10.d0) then
if(abs(rho) < 3.5d0) then
if(abs(rho) .lt. 0.35d0) then
select case(n)
case(0)
crint_2 = (((((((((1.3122532963802805073d-08 * rho &
- 1.450385222315046877d-07) * rho &
+ 1.458916900093370682d-06) * rho &
- 0.132275132275132275d-04) * rho &
+ 0.106837606837606838d-03) * rho &
- 0.757575757575757576d-03) * rho &
+ 0.462962962962962963d-02) * rho &
- 0.238095238095238095d-01) * rho &
+ 0.10000000000000000000d0) * rho &
- 0.33333333333333333333d0) * rho &
+ 1.0d0
case(1)
crint_2 = (((((((((1.198144314086343d-08 * rho &
- 1.312253296380281d-07) * rho &
+ 1.305346700083542d-06) * rho &
- 1.167133520074696d-05) * rho &
+ 9.259259259259259d-05) * rho &
- 6.410256410256410d-04) * rho &
+ 3.787878787878788d-03) * rho &
- 1.851851851851852d-02) * rho &
+ 7.142857142857142d-02) * rho &
- 2.000000000000000d-01) * rho &
+ 3.333333333333333d-01
case(2)
crint_2 = (((((((((1.102292768959436d-08 * rho &
- 1.198144314086343d-07) * rho &
+ 1.181027966742252d-06) * rho &
- 1.044277360066834d-05) * rho &
+ 8.169934640522875d-05) * rho &
- 5.555555555555556d-04) * rho &
+ 3.205128205128205d-03) * rho &
- 1.515151515151515d-02) * rho &
+ 5.555555555555555d-02) * rho &
- 1.428571428571428d-01) * rho &
+ 2.000000000000000d-01
case(3)
crint_2 = (((((((((1.020641452740218d-08 * rho &
- 1.102292768959436d-07) * rho &
+ 1.078329882677709d-06) * rho &
- 9.448223733938020d-06) * rho &
+ 7.309941520467836d-05) * rho &
- 4.901960784313725d-04) * rho &
+ 2.777777777777778d-03) * rho &
- 1.282051282051282d-02) * rho &
+ 4.545454545454546d-02) * rho &
- 1.111111111111111d-01) * rho &
+ 1.428571428571428d-01
case default
tmp = dble(n + n + 1)
crint_2 = (((((((((2.755731922398589d-07 * rho / (tmp + 20.d0) &
- 2.755731922398589d-06 / (tmp + 18.d0)) * rho &
+ 2.480158730158730d-05 / (tmp + 16.d0)) * rho &
- 1.984126984126984d-04 / (tmp + 14.d0)) * rho &
+ 1.388888888888889d-03 / (tmp + 12.d0)) * rho &
- 8.333333333333333d-03 / (tmp + 10.d0)) * rho &
+ 4.166666666666666d-02 / (tmp + 8.d0)) * rho &
- 1.666666666666667d-01 / (tmp + 6.d0)) * rho &
+ 5.000000000000000d-01 / (tmp + 4.d0)) * rho &
- 1.000000000000000d+00 / (tmp + 2.d0)) * rho &
+ 1.0d0 / tmp
end select
if(abs(rho) .lt. 1d-6) then
tmp = dble(n + n)
rho2 = rho * rho
crint_2 = - 0.16666666666666666d0 * rho * rho2 / (tmp + 7.d0) &
+ 0.5d0 * rho2 / (tmp + 5.d0) &
- rho / (tmp + 3.d0) &
+ 1.d0 / (tmp + 1.d0)
else
crint_2 = crint_smallz(n, rho)
endif
else
if(real(rho) .ge. 0.d0) then
call zboysfun(n, rho, vals)
crint_2 = vals(n)
else
call zboysfunnrp(n, rho, vals)
crint_2 = vals(n) * zexp(-rho)
endif
endif
endif
return
@ -380,29 +440,82 @@ subroutine crint_2_vec(n_max, rho, vals)
abs_rho = abs(rho)
if(abs_rho < 10.d0) then
if(abs_rho < 3.5d0) then
if(abs_rho .lt. 1d-6) then
if(abs_rho .lt. 0.35d0) then
! use finite expansion for very small rho
vals(0) = (((((((((1.3122532963802805073d-08 * rho &
- 1.450385222315046877d-07) * rho &
+ 1.458916900093370682d-06) * rho &
- 0.132275132275132275d-04) * rho &
+ 0.106837606837606838d-03) * rho &
- 0.757575757575757576d-03) * rho &
+ 0.462962962962962963d-02) * rho &
- 0.238095238095238095d-01) * rho &
+ 0.10000000000000000000d0) * rho &
- 0.33333333333333333333d0) * rho &
+ 1.0d0
!! rho^2 / 2
!rho2 = 0.5d0 * rho * rho
!! rho^3 / 6
!rho3 = 0.3333333333333333d0 * rho * rho2
!vals(0) = 1.d0 - 0.3333333333333333d0 * rho + 0.2d0 * rho2 - 0.14285714285714285d0 * rho3
!do n = 1, n_max
! tmp = 2.d0 * dble(n)
! vals(n) = 1.d0 / (tmp + 1.d0) - rho / (tmp + 3.d0) &
! + rho2 / (tmp + 5.d0) - rho3 / (tmp + 7.d0)
!enddo
if(n > 0) then
! TODO (last term)
vals(0) = 1.d0 + rho * (-0.3333333333333333d0 + rho*(0.1d0 - 0.047619047619047616d0*rho))
do n = 1, n_max
tmp = 2.d0 * dble(n)
vals(n) = 1.d0 / (tmp + 1.d0) + rho * (-1.d0/(tmp+3.d0) + rho*(0.5d0/(tmp+5.d0) - 0.3333333333333333d0*rho/(tmp+7.d0)))
enddo
vals(1) = (((((((((1.198144314086343d-08 * rho &
- 1.312253296380281d-07) * rho &
+ 1.305346700083542d-06) * rho &
- 1.167133520074696d-05) * rho &
+ 9.259259259259259d-05) * rho &
- 6.410256410256410d-04) * rho &
+ 3.787878787878788d-03) * rho &
- 1.851851851851852d-02) * rho &
+ 7.142857142857142d-02) * rho &
- 2.000000000000000d-01) * rho &
+ 3.333333333333333d-01
if(n > 1) then
vals(2) = (((((((((1.102292768959436d-08 * rho &
- 1.198144314086343d-07) * rho &
+ 1.181027966742252d-06) * rho &
- 1.044277360066834d-05) * rho &
+ 8.169934640522875d-05) * rho &
- 5.555555555555556d-04) * rho &
+ 3.205128205128205d-03) * rho &
- 1.515151515151515d-02) * rho &
+ 5.555555555555555d-02) * rho &
- 1.428571428571428d-01) * rho &
+ 2.000000000000000d-01
if(n > 2) then
vals(3) = (((((((((1.020641452740218d-08 * rho &
- 1.102292768959436d-07) * rho &
+ 1.078329882677709d-06) * rho &
- 9.448223733938020d-06) * rho &
+ 7.309941520467836d-05) * rho &
- 4.901960784313725d-04) * rho &
+ 2.777777777777778d-03) * rho &
- 1.282051282051282d-02) * rho &
+ 4.545454545454546d-02) * rho &
- 1.111111111111111d-01) * rho &
+ 1.428571428571428d-01
do n = 4, n_max
tmp = dble(n + n + 1)
vals(n) = (((((((((2.755731922398589d-07 * rho / (tmp + 20.d0) &
- 2.755731922398589d-06 / (tmp + 18.d0)) * rho &
+ 2.480158730158730d-05 / (tmp + 16.d0)) * rho &
- 1.984126984126984d-04 / (tmp + 14.d0)) * rho &
+ 1.388888888888889d-03 / (tmp + 12.d0)) * rho &
- 8.333333333333333d-03 / (tmp + 10.d0)) * rho &
+ 4.166666666666666d-02 / (tmp + 8.d0)) * rho &
- 1.666666666666667d-01 / (tmp + 6.d0)) * rho &
+ 5.000000000000000d-01 / (tmp + 4.d0)) * rho &
- 1.000000000000000d+00 / (tmp + 2.d0)) * rho &
+ 1.0d0 / tmp
enddo
endif ! n > 2
endif ! n > 1
endif ! n > 0
else
@ -445,7 +558,7 @@ subroutine crint_smallz_vec(n_max, rho, vals)
complex*16, intent(out) :: vals(0:n_max)
integer, parameter :: kmax = 40
double precision, parameter :: eps = 1.d-13
double precision, parameter :: eps = 1.d-10
integer :: k, n
complex*16 :: ct, delta_k
@ -475,11 +588,12 @@ subroutine crint_smallz_vec(n_max, rho, vals)
endif
enddo
!if(abs(delta_k) > eps) then
! write(*,*) ' pb in crint_smallz_vec !'
! write(*,*) ' n, rho = ', n, rho
! write(*,*) ' |delta_k| = ', abs(delta_k)
!endif
if(abs(delta_k) > eps) then
write(*,*) ' pb in crint_smallz_vec !'
write(*,*) ' n, rho = ', n, rho
write(*,*) ' value = ', vals(n)
write(*,*) ' |delta_k| = ', abs(delta_k)
endif
enddo
deallocate(rho_k)

64
src/utils/cpx_erf.irp.f
View File

@ -201,49 +201,6 @@ end
! ---
complex*16 function cpx_erf_2(x, y)
BEGIN_DOC
!
! compute erf(z) for z = x + i y
!
! Beylkin & Sharma, J. Chem. Phys. 155, 174117 (2021)
! https://doi.org/10.1063/5.0062444
!
END_DOC
implicit none
double precision, intent(in) :: x, y
double precision :: yabs
complex*16 :: z
yabs = dabs(y)
if(yabs .lt. 1.d-15) then
cpx_erf_2 = (1.d0, 0.d0) * derf(x)
return
else
z = x + (0.d0, 1.d0) * y
if(x .ge. 0.d0) then
call zboysfun00(z, cpx_erf_2)
else
call zboysfun00nrp(z, cpx_erf_2)
cpx_erf_2 = cpx_erf_2 * zexp(-z)
endif
endif
return
end
! ---
subroutine zboysfun00(z, val)
BEGIN_DOC
@ -333,20 +290,18 @@ subroutine zboysfun00(z, val)
complex*16, intent(out) :: val
integer :: k
complex*16 :: z1, zz, y
zz = zexp(-z)
complex*16 :: z1, y
if(abs(z) .ge. 100.0d0) then
! large |z|
z1 = (1.0d0, 0.d0) / zsqrt(z)
y = (1.0d0, 0.d0) / z
z1 = (1.d0, 0.d0) / zsqrt(z)
y = (1.d0, 0.d0) / z
val = asymcoef(7)
do k = 6, 1, -1
val = val * y + asymcoef(k)
enddo
val = zz * val * y + z1 * sqpio2
val = zexp(-z) * val * y + z1 * sqpio2
else if(abs(z) .le. 0.35d0) then
@ -359,12 +314,7 @@ subroutine zboysfun00(z, val)
else
! intermediate |z|
val = sqpio2 / zsqrt(z) - 0.5d0 * zz * sum(ff(1:22)/(z+pp(1:22)))
!val = (0.d0, 0.d0)
!do k = 1, 22
! val += ff(k) / (z + pp(k))
!enddo
!val = sqpio2 / zsqrt(z) - 0.5d0 * zz * val
val = sqpio2 / zsqrt(z) - 0.5d0 * zexp(-z) * sum(ff(1:22)/(z+pp(1:22)))
endif
@ -533,8 +483,8 @@ subroutine zboysfun00nrp(z, val)
if(abs(z) .ge. 100.0d0) then
! large |z|
z1 = (1.0d0, 0.d0) / zsqrt(z)
y = (1.0d0, 0.d0) / z
z1 = (1.d0, 0.d0) / zsqrt(z)
y = (1.d0, 0.d0) / z
val = asymcoef(7)
do k = 6, 1, -1
val = val * y + asymcoef(k)