LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 19:43:32 +01:00
Code Releases Activity

fixed bug in cGTOS overlaps

Browse Source
This commit is contained in:
AbdAmmar 2024-10-18 09:50:39 +02:00
parent ff63afa797
commit 4da6a2aea0
4 changed files with 66 additions and 55 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/ao_basis/aos.irp.f
View File

@ -75,10 +75,6 @@ END_PROVIDER
enddo
endif
powA(1) = ao_power(i,1)
powA(2) = ao_power(i,2)
powA(3) = ao_power(i,3)
! Normalization of the contracted basis functions
if (ao_normalized) then
norm = 0.d0

13
src/ao_one_e_ints/aos_cgtos.irp.f
View File

@ -4,6 +4,7 @@
BEGIN_PROVIDER [double precision, ao_coef_cgtos_norm_ord_transp, (ao_prim_num_max, ao_num)]
implicit none
integer :: i, j
do j = 1, ao_num
@ -62,9 +63,9 @@ BEGIN_PROVIDER [double precision, ao_coef_norm_cgtos, (ao_num, ao_prim_num_max)]
powA(2) = ao_power(i,2)
powA(3) = ao_power(i,3)
! Normalization of the primitives
if(primitives_normalized) then
! Normalization of the primitives
do j = 1, ao_prim_num(i)
expo = ao_expo(i,j) + (0.d0, 1.d0) * ao_expo_im(i,j)
@ -81,11 +82,15 @@ BEGIN_PROVIDER [double precision, ao_coef_norm_cgtos, (ao_num, ao_prim_num_max)]
C1 = zexp(-(0.d0, 2.d0) * phiA - 0.5d0 * expo_inv * KA2)
C2 = zexp(-(0.5d0, 0.d0) * real(expo_inv) * KA2)
call overlap_cgaussian_xyz(C_Ae, C_Ae, expo, expo, powA, powA, C_Ap, C_Ap, overlap_x, overlap_y, overlap_z, integ1, nz)
call overlap_cgaussian_xyz(conjg(C_Ae), C_Ae, conjg(expo), expo, powA, powA, conjg(C_Ap), C_Ap, overlap_x, overlap_y, overlap_z, integ2, nz)
call overlap_cgaussian_xyz(C_Ae, C_Ae, expo, expo, powA, powA, &
C_Ap, C_Ap, overlap_x, overlap_y, overlap_z, integ1, nz)
call overlap_cgaussian_xyz(conjg(C_Ae), C_Ae, conjg(expo), expo, powA, powA, &
conjg(C_Ap), C_Ap, overlap_x, overlap_y, overlap_z, integ2, nz)
norm = 2.d0 * real(C1 * integ1 + C2 * integ2)
!ao_coef_norm_cgtos(i,j) = 1.d0 / dsqrt(norm)
ao_coef_norm_cgtos(i,j) = ao_coef(i,j) / dsqrt(norm)
enddo
@ -95,7 +100,7 @@ BEGIN_PROVIDER [double precision, ao_coef_norm_cgtos, (ao_num, ao_prim_num_max)]
ao_coef_norm_cgtos(i,j) = ao_coef(i,j)
enddo
endif
endif ! primitives_normalized
enddo

46
src/utils/cgtos_one_e.irp.f
View File

@ -46,12 +46,13 @@ end
! ---
subroutine overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, &
overlap_x, overlap_y, overlap_z, overlap, dim)
subroutine overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x, overlap_y, overlap_z, overlap, dim)
BEGIN_DOC
!
! S_x = \int (x - Ap_x)^{a_x} exp(-\alpha (x - Ae_x)^2) (x - Bp_x)^{b_x} exp(-beta (x - Be_x)^2) dx
! S_x = \int (x - Ap_x)^{a_x} exp(-\alpha (x - Ae_x)^2)
! (x - Bp_x)^{b_x} exp(-\beta (x - Be_x)^2) dx
!
! S = S_x S_y S_z
!
@ -70,8 +71,9 @@ subroutine overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, pow
integer :: i, nmax, iorder_p(3)
complex*16 :: P_new(0:max_dim,3), P_center(3), fact_p, p, inv_sq_p
complex*16 :: F_integral_tab(0:max_dim)
complex*16 :: ab, arg
complex*16 :: Fc_integral
complex*16, external :: Fc_integral
call give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_p, iorder_p, &
alpha, beta, power_A, power_B, Ae_center, Be_center, Ap_center, Bp_center, dim)
@ -85,36 +87,52 @@ subroutine overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, pow
endif
nmax = maxval(iorder_p)
inv_sq_p = (1.d0, 0.d0) / zsqrt(p)
do i = 0, nmax
F_integral_tab(i) = Fc_integral(i, inv_sq_p)
enddo
overlap_x = P_new(0,1) * F_integral_tab(0)
overlap_y = P_new(0,2) * F_integral_tab(0)
overlap_z = P_new(0,3) * F_integral_tab(0)
ab = alpha * beta * inv_sq_p
arg = ab * (Ae_center(1) - Be_center(1)) &
* (Ae_center(1) - Be_center(1))
if(real(arg) > 40.d0) then
overlap_x = (0.d0, 0.d0)
else
overlap_x = P_new(0,1) * F_integral_tab(0)
do i = 1, iorder_p(1)
overlap_x = overlap_x + P_new(i,1) * F_integral_tab(i)
enddo
call cgaussian_product_x(alpha, Ap_center(1), beta, Bp_center(1), fact_p, p, P_center(1))
overlap_x *= fact_p
overlap_x = overlap_x * zexp(-arg)
endif
arg = ab * (Ae_center(2) - Be_center(2)) &
* (Ae_center(2) - Be_center(2))
if(real(arg) > 40.d0) then
overlap_y = (0.d0, 0.d0)
else
overlap_y = P_new(0,2) * F_integral_tab(0)
do i = 1, iorder_p(2)
overlap_y = overlap_y + P_new(i,2) * F_integral_tab(i)
enddo
call cgaussian_product_x(alpha, Ap_center(2), beta, Bp_center(2), fact_p, p, P_center(2))
overlap_y *= fact_p
overlap_y = overlap_y * zexp(-arg)
endif
arg = ab * (Ae_center(3) - Be_center(3)) &
* (Ae_center(3) - Be_center(3))
if(real(arg) > 40.d0) then
overlap_z = (0.d0, 0.d0)
else
overlap_z = P_new(0,3) * F_integral_tab(0)
do i = 1, iorder_p(3)
overlap_z = overlap_z + P_new(i,3) * F_integral_tab(i)
enddo
call cgaussian_product_x(alpha, Ap_center(3), beta, Bp_center(3), fact_p, p, P_center(3))
overlap_z *= fact_p
overlap_z = overlap_z * zexp(-arg)
endif
overlap = overlap_x * overlap_y * overlap_z
return
end
! ---

24
src/utils/cgtos_utils.irp.f
View File

@ -91,8 +91,8 @@ subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder,
include 'constants.include.F'
integer, intent(in) :: dim, a(3), b(3)
complex*16, intent(in) :: alpha, Ap_center(3), Ae_center(3)
complex*16, intent(in) :: beta, Bp_center(3), Be_center(3)
complex*16, intent(in) :: alpha, Ae_center(3), Ap_center(3)
complex*16, intent(in) :: beta, Be_center(3), Bp_center(3)
integer, intent(out) :: iorder(3)
complex*16, intent(out) :: p, P_center(3), fact_k, P_new(0:max_dim,3)
@ -167,13 +167,12 @@ subroutine cgaussian_product(a, xa, b, xb, k, p, xp)
complex*16, intent(in) :: a, b, xa(3), xb(3)
complex*16, intent(out) :: p, k, xp(3)
double precision :: tmp_mod
complex*16 :: p_inv, xab(3), ab
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xab
ASSERT (REAL(a) > 0.)
ASSERT (REAL(b) > 0.)
ASSERT (real(a) > 0.)
ASSERT (real(b) > 0.)
! new exponent
p = a + b
@ -186,8 +185,7 @@ subroutine cgaussian_product(a, xa, b, xb, k, p, xp)
ab = a * b * p_inv
k = ab * (xab(1)*xab(1) + xab(2)*xab(2) + xab(3)*xab(3))
tmp_mod = dsqrt(real(k)*real(k) + aimag(k)*aimag(k))
if(tmp_mod .gt. 40.d0) then
if(real(k) .gt. 40.d0) then
k = (0.d0, 0.d0)
xp(1:3) = (0.d0, 0.d0)
return
@ -229,7 +227,7 @@ subroutine cgaussian_product_x(a, xa, b, xb, k, p, xp)
ab = a * b * p_inv
k = ab * xab * xab
if(zabs(k) > 40.d0) then
if(real(k) > 40.d0) then
k = (0.d0, 0.d0)
xp = (0.d0, 0.d0)
return
@ -300,7 +298,6 @@ subroutine add_cpoly(b, nb, c, nc, d, nd)
complex*16, intent(out) :: d(0:nb+nc)
integer :: ib
double precision :: tmp_mod
complex*16 :: tmp
nd = nb + nc
@ -309,11 +306,9 @@ subroutine add_cpoly(b, nb, c, nc, d, nd)
enddo
tmp = d(nd)
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp))
do while( (tmp_mod .lt. 1.d-15) .and. (nd >= 0) )
do while( (zabs(tmp) .lt. 1.d-15) .and. (nd >= 0) )
nd -= 1
tmp = d(nd)
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp))
if(nd < 0) exit
enddo
@ -336,7 +331,6 @@ subroutine add_cpoly_multiply(b, nb, cst, d, nd)
complex*16, intent(inout) :: d(0:max(nb, nd))
integer :: ib
double precision :: tmp_mod
complex*16 :: tmp
nd = max(nd, nb)
@ -347,12 +341,10 @@ subroutine add_cpoly_multiply(b, nb, cst, d, nd)
enddo
tmp = d(nd)
tmp_mod = dsqrt(real(tmp)*real(tmp) + aimag(tmp)*aimag(tmp))
do while(tmp_mod .lt. 1.d-15)
do while(zabs(tmp) .lt. 1.d-15)
nd -= 1
if(nd < 0) exit
tmp = d(nd)
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp))
enddo
endif