From 45481ac08e3b8377161650ff57b7972108f0a7d3 Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Fri, 18 Oct 2024 13:42:51 +0200 Subject: [PATCH] fixed pw-centering-related bug in cGTOs kinetic integrals --- src/ao_one_e_ints/aos_cgtos.irp.f | 8 +- .../one_e_kin_integrals_cgtos.irp.f | 232 +++++++++++------- src/utils/cgtos_one_e.irp.f | 8 +- 3 files changed, 157 insertions(+), 91 deletions(-) diff --git a/src/ao_one_e_ints/aos_cgtos.irp.f b/src/ao_one_e_ints/aos_cgtos.irp.f index f89dfd09..f926c846 100644 --- a/src/ao_one_e_ints/aos_cgtos.irp.f +++ b/src/ao_one_e_ints/aos_cgtos.irp.f @@ -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) diff --git a/src/ao_one_e_ints/one_e_kin_integrals_cgtos.irp.f b/src/ao_one_e_ints/one_e_kin_integrals_cgtos.irp.f index b55c37f8..f2557fc3 100644 --- a/src/ao_one_e_ints/one_e_kin_integrals_cgtos.irp.f +++ b/src/ao_one_e_ints/one_e_kin_integrals_cgtos.irp.f @@ -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 diff --git a/src/utils/cgtos_one_e.irp.f b/src/utils/cgtos_one_e.irp.f index dffb4e47..83616969 100644 --- a/src/utils/cgtos_one_e.irp.f +++ b/src/utils/cgtos_one_e.irp.f @@ -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