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

fixed pw-centering-related bug in cGTOs kinetic integrals

Browse Source
This commit is contained in:
AbdAmmar 2024-10-18 13:42:51 +02:00
parent 4da6a2aea0
commit 45481ac08e
3 changed files with 157 additions and 91 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -246,11 +246,11 @@ END_PROVIDER
C2(3) = zexp((0.d0, 1.d0) * (phiA(3) - phiB(3)) - 0.25d0 * (conjg(alpha_inv) * KA2(3) + beta_inv * KB2(3)))
C2(4) = C2(1) * C2(2) * C2(3)
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, &
overlap_x1, overlap_y1, overlap_z1, overlap1, dim1)
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x1, overlap_y1, overlap_z1, overlap1, dim1)
call overlap_cgaussian_xyz(conjg(Ae_center), Be_center, conjg(alpha), beta, power_A, power_B, conjg(Ap_center), Bp_center, &
overlap_x2, overlap_y2, overlap_z2, overlap2, dim1)
call overlap_cgaussian_xyz(conjg(Ae_center), Be_center, conjg(alpha), beta, power_A, power_B, &
conjg(Ap_center), Bp_center, overlap_x2, overlap_y2, overlap_z2, overlap2, dim1)
overlap_x = 2.d0 * real(C1(1) * overlap_x1 + C2(1) * overlap_x2)
overlap_y = 2.d0 * real(C1(2) * overlap_y1 + C2(2) * overlap_y2)

232
src/ao_one_e_ints/one_e_kin_integrals_cgtos.irp.f
View File

@ -10,13 +10,15 @@
double precision :: c, deriv_tmp
double precision :: KA2, phiA
double precision :: KB2, phiB
double precision :: aa
complex*16 :: alpha, alpha_inv, Ae_center(3), Ap_center(3), C1
complex*16 :: beta, beta_inv, Be_center(3), Bp_center(3), C2
complex*16 :: xa
complex*16 :: overlap_x, overlap_y, overlap_z, overlap
complex*16 :: overlap_x0_1, overlap_y0_1, overlap_z0_1
complex*16 :: overlap_x0_2, overlap_y0_2, overlap_z0_2
complex*16 :: overlap_m2_1, overlap_p2_1
complex*16 :: overlap_m2_2, overlap_p2_2
complex*16 :: overlap_m2_1, overlap_m1_1, overlap_p1_1, overlap_p2_1
complex*16 :: overlap_m2_2, overlap_m1_2, overlap_p1_2, overlap_p2_2
complex*16 :: deriv_tmp_1, deriv_tmp_2
@ -32,26 +34,27 @@
beta = (0.1d0, 0.d0)
power_A = 1
power_B = 0
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, &
overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1)
! ---
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(i, j, m, n, l, ii, jj, c, C1, C2, &
!$OMP Ae_center, power_A, alpha, alpha_inv, KA2, phiA, Ap_center, &
!$OMP Be_center, power_B, beta, beta_inv, KB2, phiB, Bp_center, &
!$OMP deriv_tmp, deriv_tmp_1, deriv_tmp_2, &
!$OMP overlap_x, overlap_y, overlap_z, overlap, &
!$OMP overlap_m2_1, overlap_p2_1, overlap_m2_2, overlap_p2_2, &
!$OMP overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap_x0_2, &
!$OMP overlap_y0_2, overlap_z0_2) &
!$OMP SHARED(nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1, &
!$OMP ao_coef_cgtos_norm_ord_transp, ao_expo_cgtos_ord_transp, &
!$OMP ao_expo_pw_ord_transp, ao_expo_phase_ord_transp, &
!$OMP ao_deriv2_cgtos_x, ao_deriv2_cgtos_y, ao_deriv2_cgtos_z)
!$OMP PARALLEL &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(i, j, m, n, l, ii, jj, c, aa, xa, C1, C2, &
!$OMP Ae_center, power_A, alpha, alpha_inv, KA2, phiA, Ap_center, &
!$OMP Be_center, power_B, beta, beta_inv, KB2, phiB, Bp_center, &
!$OMP deriv_tmp, deriv_tmp_1, deriv_tmp_2, &
!$OMP overlap_x, overlap_y, overlap_z, overlap, &
!$OMP overlap_m2_1, overlap_m1_1, overlap_p1_1, overlap_p2_1, &
!$OMP overlap_m2_2, overlap_m1_2, overlap_p1_2, overlap_p2_2, &
!$OMP overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap_x0_2, &
!$OMP overlap_y0_2, overlap_z0_2) &
!$OMP SHARED(nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1, &
!$OMP ao_coef_cgtos_norm_ord_transp, ao_expo_cgtos_ord_transp, &
!$OMP ao_expo_pw_ord_transp, ao_expo_phase_ord_transp, &
!$OMP ao_deriv2_cgtos_x, ao_deriv2_cgtos_y, ao_deriv2_cgtos_z)
!$OMP DO SCHEDULE(GUIDED)
do j = 1, ao_num
jj = ao_nucl(j)
@ -97,42 +100,62 @@
C1 = zexp((0.d0, 1.d0) * (-phiA - phiB) - 0.25d0 * (alpha_inv * KA2 + beta_inv * KB2))
C2 = zexp((0.d0, 1.d0) * (-phiA + phiB) - 0.25d0 * (alpha_inv * KA2 + conjg(beta_inv) * KB2))
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, &
overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), &
overlap_x0_2, overlap_y0_2, overlap_z0_2, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_x0_2, overlap_y0_2, overlap_z0_2, overlap, dim1)
! ---
power_A(1) = power_A(1) - 2
power_A(1) = power_A(1) - 1
if(power_A(1) > -1) then
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, &
overlap_m2_1, overlap_y, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), &
overlap_m2_2, overlap_y, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_m1_1, overlap_y, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_m1_2, overlap_y, overlap_z, overlap, dim1)
power_A(1) = power_A(1) - 1
if(power_A(1) > -1) then
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_m2_1, overlap_y, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_m2_2, overlap_y, overlap_z, overlap, dim1)
else
overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0)
endif
power_A(1) = power_A(1) + 1
else
overlap_m1_1 = (0.d0, 0.d0)
overlap_m1_2 = (0.d0, 0.d0)
overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0)
endif
power_A(1) = power_A(1) + 1
power_A(1) = power_A(1) + 4
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, &
overlap_p2_1, overlap_y, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), &
overlap_p2_2, overlap_y, overlap_z, overlap, dim1)
power_A(1) = power_A(1) + 1
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_p1_1, overlap_y, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_p1_2, overlap_y, overlap_z, overlap, dim1)
power_A(1) = power_A(1) + 1
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_p2_1, overlap_y, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_p2_2, overlap_y, overlap_z, overlap, dim1)
power_A(1) = power_A(1) - 2
deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * dble(power_A(1)) + 1.d0) * overlap_x0_1 &
+ dble(power_A(1)) * (dble(power_A(1)) - 1.d0) * overlap_m2_1 &
+ 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_y0_1 * overlap_z0_1
aa = dble(power_A(1))
xa = Ap_center(1) - Ae_center(1)
deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * dble(power_A(1)) + 1.d0) * overlap_x0_2 &
+ dble(power_A(1)) * (dble(power_A(1)) - 1.d0) * overlap_m2_2 &
+ 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_y0_2 * overlap_z0_2
deriv_tmp_1 = aa * (aa - 1.d0) * overlap_m2_1 - 4.d0 * alpha * aa * xa * overlap_m1_1 &
+ 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_x0_1 &
+ 8.d0 * alpha * alpha * (xa * overlap_p1_1 + 0.5d0 * overlap_p2_1)
deriv_tmp_1 = deriv_tmp_1 * overlap_y0_1 * overlap_z0_1
deriv_tmp_2 = aa * (aa - 1.d0) * overlap_m2_2 - 4.d0 * alpha * aa * xa * overlap_m1_2 &
+ 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_x0_2 &
+ 8.d0 * alpha * alpha * (xa * overlap_p1_2 + 0.5d0 * overlap_p2_2)
deriv_tmp_2 = deriv_tmp_2 * overlap_y0_2 * overlap_z0_2
deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2)
@ -140,34 +163,55 @@
! ---
power_A(2) = power_A(2) - 2
power_A(2) = power_A(2) - 1
if(power_A(2) > -1) then
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, &
overlap_x, overlap_m2_1, overlap_y, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), &
overlap_x, overlap_m2_2, overlap_y, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x, overlap_m1_1, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_x, overlap_m1_2, overlap_z, overlap, dim1)
power_A(2) = power_A(2) - 1
if(power_A(2) > -1) then
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x, overlap_m2_1, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_x, overlap_m2_2, overlap_z, overlap, dim1)
else
overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0)
endif
power_A(2) = power_A(2) + 1
else
overlap_m1_1 = (0.d0, 0.d0)
overlap_m1_2 = (0.d0, 0.d0)
overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0)
endif
power_A(2) = power_A(2) + 1
power_A(2) = power_A(2) + 4
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, &
overlap_x, overlap_p2_1, overlap_y, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), &
overlap_x, overlap_p2_2, overlap_y, overlap, dim1)
power_A(2) = power_A(2) + 1
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x, overlap_p1_1, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_x, overlap_p1_2, overlap_z, overlap, dim1)
power_A(2) = power_A(2) + 1
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x, overlap_p2_1, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_x, overlap_p2_2, overlap_z, overlap, dim1)
power_A(2) = power_A(2) - 2
deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * dble(power_A(2)) + 1.d0) * overlap_y0_1 &
+ dble(power_A(2)) * (dble(power_A(2)) - 1.d0) * overlap_m2_1 &
+ 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_x0_1 * overlap_z0_1
aa = dble(power_A(2))
xa = Ap_center(2) - Ae_center(2)
deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * dble(power_A(2)) + 1.d0) * overlap_y0_2 &
+ dble(power_A(2)) * (dble(power_A(2)) - 1.d0) * overlap_m2_2 &
+ 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_x0_2 * overlap_z0_2
deriv_tmp_1 = aa * (aa - 1.d0) * overlap_m2_1 - 4.d0 * alpha * aa * xa * overlap_m1_1 &
+ 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_y0_1 &
+ 8.d0 * alpha * alpha * (xa * overlap_p1_1 + 0.5d0 * overlap_p2_1)
deriv_tmp_1 = deriv_tmp_1 * overlap_x0_1 * overlap_z0_1
deriv_tmp_2 = aa * (aa - 1.d0) * overlap_m2_2 - 4.d0 * alpha * aa * xa * overlap_m1_2 &
+ 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_y0_2 &
+ 8.d0 * alpha * alpha * (xa * overlap_p1_2 + 0.5d0 * overlap_p2_2)
deriv_tmp_2 = deriv_tmp_2 * overlap_x0_2 * overlap_z0_2
deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2)
@ -175,34 +219,55 @@
! ---
power_A(3) = power_A(3) - 2
power_A(3) = power_A(3) - 1
if(power_A(3) > -1) then
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, &
overlap_x, overlap_y, overlap_m2_1, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), &
overlap_x, overlap_y, overlap_m2_2, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x, overlap_y, overlap_m1_1, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_m1_2, overlap, dim1)
power_A(3) = power_A(3) - 1
if(power_A(3) > -1) then
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x, overlap_y, overlap_m2_1, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_m2_2, overlap, dim1)
else
overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0)
endif
power_A(3) = power_A(3) + 1
else
overlap_m1_1 = (0.d0, 0.d0)
overlap_m1_2 = (0.d0, 0.d0)
overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0)
endif
power_A(3) = power_A(3) + 1
power_A(3) = power_A(3) + 4
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, &
overlap_x, overlap_y, overlap_p2_1, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), &
overlap_x, overlap_y, overlap_p2_2, overlap, dim1)
power_A(3) = power_A(3) + 1
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x, overlap_y, overlap_p1_1, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_p1_2, overlap, dim1)
power_A(3) = power_A(3) + 1
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x, overlap_y, overlap_p2_1, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_p2_2, overlap, dim1)
power_A(3) = power_A(3) - 2
deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * dble(power_A(3)) + 1.d0) * overlap_z0_1 &
+ dble(power_A(3)) * (dble(power_A(3)) - 1.d0) * overlap_m2_1 &
+ 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_x0_1 * overlap_y0_1
deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * dble(power_A(3)) + 1.d0) * overlap_z0_2 &
+ dble(power_A(3)) * (dble(power_A(3)) - 1.d0) * overlap_m2_2 &
+ 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_x0_2 * overlap_y0_2
aa = dble(power_A(3))
xa = Ap_center(3) - Ae_center(3)
deriv_tmp_1 = aa * (aa - 1.d0) * overlap_m2_1 - 4.d0 * alpha * aa * xa * overlap_m1_1 &
+ 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_z0_1 &
+ 8.d0 * alpha * alpha * (xa * overlap_p1_1 + 0.5d0 * overlap_p2_1)
deriv_tmp_1 = deriv_tmp_1 * overlap_x0_1 * overlap_y0_1
deriv_tmp_2 = aa * (aa - 1.d0) * overlap_m2_2 - 4.d0 * alpha * aa * xa * overlap_m1_2 &
+ 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_z0_2 &
+ 8.d0 * alpha * alpha * (xa * overlap_p1_2 + 0.5d0 * overlap_p2_2)
deriv_tmp_2 = deriv_tmp_2 * overlap_x0_2 * overlap_y0_2
deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2)
@ -214,7 +279,8 @@
enddo
enddo
enddo
!$OMP END PARALLEL DO
!$OMP END DO
!$OMP END PARALLEL
END_PROVIDER

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

@ -79,10 +79,10 @@ subroutine overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, pow
alpha, beta, power_A, power_B, Ae_center, Be_center, Ap_center, Bp_center, dim)
if(zabs(fact_p) .lt. 1.d-14) then
overlap_x = (1.d-10, 0.d0)
overlap_y = (1.d-10, 0.d0)
overlap_z = (1.d-10, 0.d0)
overlap = (1.d-10, 0.d0)
overlap_x = (0.d0, 0.d0)
overlap_y = (0.d0, 0.d0)
overlap_z = (0.d0, 0.d0)
overlap = (0.d0, 0.d0)
return
endif