diff --git a/src/ao_basis/aos.irp.f b/src/ao_basis/aos.irp.f index 1cbd3976..09604487 100644 --- a/src/ao_basis/aos.irp.f +++ b/src/ao_basis/aos.irp.f @@ -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 diff --git a/src/ao_one_e_ints/aos_cgtos.irp.f b/src/ao_one_e_ints/aos_cgtos.irp.f index b5dd5263..f89dfd09 100644 --- a/src/ao_one_e_ints/aos_cgtos.irp.f +++ b/src/ao_one_e_ints/aos_cgtos.irp.f @@ -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 diff --git a/src/utils/cgtos_one_e.irp.f b/src/utils/cgtos_one_e.irp.f index a8773ce0..dffb4e47 100644 --- a/src/utils/cgtos_one_e.irp.f +++ b/src/utils/cgtos_one_e.irp.f @@ -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,11 +71,12 @@ 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) + 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) @@ -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 - 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 + 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 + overlap_x = overlap_x * zexp(-arg) + endif - 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 + 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 + overlap_y = overlap_y * zexp(-arg) + endif - 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 + 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 + overlap_z = overlap_z * zexp(-arg) + endif overlap = overlap_x * overlap_y * overlap_z + return end ! --- diff --git a/src/utils/cgtos_utils.irp.f b/src/utils/cgtos_utils.irp.f index 40c8a838..075d9bbc 100644 --- a/src/utils/cgtos_utils.irp.f +++ b/src/utils/cgtos_utils.irp.f @@ -66,7 +66,7 @@ end ! --- subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder, & - alpha, beta, a, b, Ae_center, Be_center, Ap_center, Bp_center, dim) + alpha, beta, a, b, Ae_center, Be_center, Ap_center, Bp_center, dim) BEGIN_DOC ! @@ -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 @@ -185,9 +184,8 @@ subroutine cgaussian_product(a, xa, b, xb, k, p, xp) p_inv = (1.d0, 0.d0) / p 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 + k = ab * (xab(1)*xab(1) + xab(2)*xab(2) + xab(3)*xab(3)) + if(real(k) .gt. 40.d0) then k = (0.d0, 0.d0) xp(1:3) = (0.d0, 0.d0) return @@ -228,8 +226,8 @@ subroutine cgaussian_product_x(a, xa, b, xb, k, p, xp) p_inv = (1.d0, 0.d0) / p ab = a * b * p_inv - k = ab * xab*xab - if(zabs(k) > 40.d0) then + k = ab * xab * xab + 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 @@ -308,12 +305,10 @@ subroutine add_cpoly(b, nb, c, nc, d, nd) d(ib) = d(ib) + c(ib) + b(ib) 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) ) + tmp = d(nd) + 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)) + tmp = d(nd) 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) @@ -346,13 +340,11 @@ subroutine add_cpoly_multiply(b, nb, cst, d, nd) d(ib) = d(ib) + cst * b(ib) enddo - tmp = d(nd) - tmp_mod = dsqrt(real(tmp)*real(tmp) + aimag(tmp)*aimag(tmp)) - do while(tmp_mod .lt. 1.d-15) + tmp = d(nd) + 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)) + tmp = d(nd) enddo endif