diff --git a/src/ao_basis/EZFIO.cfg b/src/ao_basis/EZFIO.cfg index 7f3eb7ec..ad36ca6b 100644 --- a/src/ao_basis/EZFIO.cfg +++ b/src/ao_basis/EZFIO.cfg @@ -82,12 +82,12 @@ interface: ezfio, provider [ao_expo_pw] type: double precision doc: plane wave part for each primitive GTOs |AO| -size: (3,ao_basis.ao_num,ao_basis.ao_prim_num_max) +size: (4,ao_basis.ao_num,ao_basis.ao_prim_num_max) interface: ezfio, provider [ao_expo_phase] type: double precision doc: phase shift for each primitive GTOs |AO| -size: (3,ao_basis.ao_num,ao_basis.ao_prim_num_max) +size: (4,ao_basis.ao_num,ao_basis.ao_prim_num_max) interface: ezfio, provider diff --git a/src/ao_one_e_ints/aos_cgtos.irp.f b/src/ao_one_e_ints/aos_cgtos.irp.f index c9713b7d..2792d938 100644 --- a/src/ao_one_e_ints/aos_cgtos.irp.f +++ b/src/ao_one_e_ints/aos_cgtos.irp.f @@ -30,9 +30,9 @@ END_PROVIDER ao_expo_pw_ord_transp(m,i,j) = ao_expo_pw_ord(m,j,i) ao_expo_phase_ord_transp(m,i,j) = ao_expo_phase_ord(m,j,i) enddo - ao_expo_pw_ord_transp(4,i,j) = ao_expo_pw_ord_transp(1,i,j) & - + ao_expo_pw_ord_transp(2,i,j) & - + ao_expo_pw_ord_transp(3,i,j) + ao_expo_pw_ord_transp(4,i,j) = ao_expo_pw_ord_transp(1,i,j) * ao_expo_pw_ord_transp(1,i,j) & + + ao_expo_pw_ord_transp(2,i,j) * ao_expo_pw_ord_transp(2,i,j) & + + ao_expo_pw_ord_transp(3,i,j) * ao_expo_pw_ord_transp(3,i,j) ao_expo_phase_ord_transp(4,i,j) = ao_expo_phase_ord_transp(1,j,i) & + ao_expo_phase_ord_transp(2,j,i) & + ao_expo_phase_ord_transp(3,j,i) @@ -47,10 +47,12 @@ BEGIN_PROVIDER [double precision, ao_coef_norm_cgtos, (ao_num, ao_prim_num_max)] implicit none - integer :: i, j, powA(3), nz + integer :: i, j, ii, m, powA(3), nz double precision :: norm - complex*16 :: overlap_x, overlap_y, overlap_z, C_A(3) - complex*16 :: integ1, integ2, expo + double precision :: kA2, phiA + complex*16 :: expo, expo_inv, C_A(3) + complex*16 :: overlap_x, overlap_y, overlap_z + complex*16 :: integ1, integ2, C1, C2 nz = 100 @@ -62,22 +64,31 @@ BEGIN_PROVIDER [double precision, ao_coef_norm_cgtos, (ao_num, ao_prim_num_max)] do i = 1, ao_num + ii = ao_nucl(i) powA(1) = ao_power(i,1) powA(2) = ao_power(i,2) powA(3) = ao_power(i,3) - ! TODO ! Normalization of the primitives if(primitives_normalized) then do j = 1, ao_prim_num(i) expo = ao_expo(i,j) + (0.d0, 1.d0) * ao_expo_im_cgtos(i,j) + expo_inv = (1.d0, 0.d0) / expo + do m = 1, 3 + C_A(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo_inv * ao_expo_pw(m,i,j) + enddo + phiA = ao_expo_phase(4,i,j) + KA2 = ao_expo_pw(4,i,j) - call overlap_cgaussian_xyz(C_A, C_A, expo, expo, powA, powA, overlap_x, overlap_y, overlap_z, integ1, nz) - call overlap_cgaussian_xyz(C_A, C_A, conjg(expo), expo, powA, powA, overlap_x, overlap_y, overlap_z, integ2, nz) + C1 = zexp(-(0.d0, 2.d0) * phiA - 0.5d0 * expo_inv * KA2) + C2 = zexp(-(0.5d0, 0.d0) * real(expo_inv) * KA2) - norm = 2.d0 * real(integ1 + integ2) + call overlap_cgaussian_xyz(C_A, C_A, expo, expo, powA, powA, overlap_x, overlap_y, overlap_z, integ1, nz) + call overlap_cgaussian_xyz(conjg(C_A), C_A, conjg(expo), expo, powA, powA, overlap_x, overlap_y, overlap_z, integ2, nz) + + norm = 2.d0 * real(C1 * integ1 + C2 * integ2) ao_coef_norm_cgtos(i,j) = ao_coef(i,j) / dsqrt(norm) enddo @@ -98,14 +109,14 @@ END_PROVIDER BEGIN_PROVIDER [double precision, ao_coef_norm_cgtos_ord, (ao_num, ao_prim_num_max)] &BEGIN_PROVIDER [complex*16 , ao_expo_cgtos_ord, (ao_num, ao_prim_num_max)] -&BEGIN_PROVIDER [double precision, ao_expo_pw_ord, (3, ao_num, ao_prim_num_max)] -&BEGIN_PROVIDER [double precision, ao_expo_phase_ord, (3, ao_num, ao_prim_num_max)] +&BEGIN_PROVIDER [double precision, ao_expo_pw_ord, (4, ao_num, ao_prim_num_max)] +&BEGIN_PROVIDER [double precision, ao_expo_phase_ord, (4, ao_num, ao_prim_num_max)] implicit none - integer :: i, j + integer :: i, j, m integer :: iorder(ao_prim_num_max) - double precision :: d(ao_prim_num_max,9) + double precision :: d(ao_prim_num_max,11) d = 0.d0 @@ -116,28 +127,26 @@ END_PROVIDER d(j,1) = ao_expo(i,j) d(j,2) = ao_coef_norm_cgtos(i,j) d(j,3) = ao_expo_im_cgtos(i,j) - d(j,4) = ao_expo_pw(1,i,j) - d(j,5) = ao_expo_pw(2,i,j) - d(j,6) = ao_expo_pw(3,i,j) - d(j,7) = ao_expo_phase(1,i,j) - d(j,8) = ao_expo_phase(2,i,j) - d(j,9) = ao_expo_phase(3,i,j) + + do m = 1, 4 + d(j,3+m) = ao_expo_pw(m,i,j) + d(j,7+m) = ao_expo_phase(m,i,j) + enddo enddo call dsort(d(1,1), iorder, ao_prim_num(i)) - do j = 2, 9 + do j = 2, 11 call dset_order(d(1,j), iorder, ao_prim_num(i)) enddo do j = 1, ao_prim_num(i) ao_expo_cgtos_ord (i,j) = d(j,1) + (0.d0, 1.d0) * d(j,3) ao_coef_norm_cgtos_ord(i,j) = d(j,2) - ao_expo_pw_ord(i,j,1) = d(j,4) - ao_expo_pw_ord(i,j,2) = d(j,5) - ao_expo_pw_ord(i,j,3) = d(j,6) - ao_expo_phase_ord(i,j,1) = d(j,7) - ao_expo_phase_ord(i,j,2) = d(j,8) - ao_expo_phase_ord(i,j,3) = d(j,9) + + do m = 1, 4 + ao_expo_pw_ord(m,i,j) = d(j,3+m) + ao_expo_phase_ord(m,i,j) = d(j,7+m) + enddo enddo enddo @@ -154,8 +163,10 @@ END_PROVIDER integer :: i, j, m, n, l, ii, jj, dim1, power_A(3), power_B(3) double precision :: c, overlap, overlap_x, overlap_y, overlap_z - complex*16 :: alpha, alpha_inv, A_center(3), KA2(3), phiA(3) - complex*16 :: beta, beta_inv, B_center(3), KB2(3), phiB(3) + double precision :: KA2(3), phiA(3) + double precision :: KB2(3), phiB(3) + complex*16 :: alpha, alpha_inv, A_center(3) + complex*16 :: beta, beta_inv, B_center(3) complex*16 :: C1(1:4), C2(1:4) complex*16 :: overlap1, overlap_x1, overlap_y1, overlap_z1 complex*16 :: overlap2, overlap_x2, overlap_y2, overlap_z2 @@ -199,7 +210,6 @@ END_PROVIDER alpha = ao_expo_cgtos_ord_transp(n,j) alpha_inv = (1.d0, 0.d0) / alpha - do m = 1, 3 phiA(m) = ao_expo_phase_ord_transp(m,n,j) A_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_expo_pw_ord_transp(m,n,j) @@ -210,7 +220,6 @@ END_PROVIDER beta = ao_expo_cgtos_ord_transp(l,i) beta_inv = (1.d0, 0.d0) / beta - do m = 1, 3 phiB(m) = ao_expo_phase_ord_transp(m,l,i) B_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_expo_pw_ord_transp(m,l,i) @@ -232,7 +241,7 @@ END_PROVIDER call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & overlap_x1, overlap_y1, overlap_z1, overlap1, dim1) - call overlap_cgaussian_xyz(A_center, B_center, conjg(alpha), beta, power_A, power_B, & + call overlap_cgaussian_xyz(conjg(A_center), B_center, conjg(alpha), beta, power_A, power_B, & overlap_x2, overlap_y2, overlap_z2, overlap2, dim1) overlap_x = 2.d0 * real(C1(1) * overlap_x1 + C2(1) * overlap_x2) diff --git a/src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f b/src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f index b7a9c6fe..568d4d8f 100644 --- a/src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f +++ b/src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f @@ -15,8 +15,10 @@ BEGIN_PROVIDER [double precision, ao_integrals_n_e_cgtos, (ao_num, ao_num)] integer :: power_A(3), power_B(3) integer :: i, j, k, l, m, n, ii, jj double precision :: c, Z, C_center(3) - complex*16 :: alpha, alpha_inv, A_center(3), phiA, KA2 - complex*16 :: beta, beta_inv, B_center(3), phiB, KB2 + double precision :: phiA, KA2 + double precision :: phiB, KB2 + complex*16 :: alpha, alpha_inv, A_center(3) + complex*16 :: beta, beta_inv, B_center(3) complex*16 :: C1, C2, I1, I2 complex*16 :: NAI_pol_mult_cgtos @@ -54,7 +56,7 @@ BEGIN_PROVIDER [double precision, ao_integrals_n_e_cgtos, (ao_num, ao_num)] A_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_expo_pw_ord_transp(m,n,j) enddo phiA = ao_expo_phase_ord_transp(4,n,j) - KA2 = ao_expo_pw_ord_transp(4,n,j) * ao_expo_pw_ord_transp(4,n,j) + KA2 = ao_expo_pw_ord_transp(4,n,j) do l = 1, ao_prim_num(i) @@ -65,7 +67,7 @@ BEGIN_PROVIDER [double precision, ao_integrals_n_e_cgtos, (ao_num, ao_num)] B_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_expo_pw_ord_transp(m,l,i) enddo phiB = ao_expo_phase_ord_transp(4,l,i) - KB2 = ao_expo_pw_ord_transp(4,l,i) * ao_expo_pw_ord_transp(4,l,i) + KB2 = ao_expo_pw_ord_transp(4,l,i) 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 * (conjg(alpha_inv) * KA2 + beta_inv * KB2)) @@ -79,7 +81,7 @@ BEGIN_PROVIDER [double precision, ao_integrals_n_e_cgtos, (ao_num, ao_num)] I1 = NAI_pol_mult_cgtos(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_max_integrals) - I2 = NAI_pol_mult_cgtos(A_center, B_center, power_A, power_B, conjg(alpha), beta, C_center, n_pt_max_integrals) + I2 = NAI_pol_mult_cgtos(conjg(A_center), B_center, power_A, power_B, conjg(alpha), beta, C_center, n_pt_max_integrals) c = c - Z * 2.d0 * real(C1 * I1 + C2 * I2) enddo 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 729f0a82..7f1f62cb 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 @@ -8,8 +8,10 @@ implicit none integer :: i, j, m, n, l, ii, jj, dim1, power_A(3), power_B(3) double precision :: c, deriv_tmp - complex*16 :: alpha, alpha_inv, A_center(3), KA2, phiA, C1 - complex*16 :: beta, beta_inv, B_center(3), KB2, phiB, C2 + double precision :: KA2, phiA + double precision :: KB2, phiB + complex*16 :: alpha, alpha_inv, A_center(3), C1 + complex*16 :: beta, beta_inv, B_center(3), C2 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 @@ -70,33 +72,31 @@ alpha = ao_expo_cgtos_ord_transp(n,j) alpha_inv = (1.d0, 0.d0) / alpha - do m = 1, 3 A_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_expo_pw_ord_transp(m,n,j) enddo phiA = ao_expo_phase_ord_transp(4,n,j) - KA2 = ao_expo_pw_ord_transp(4,n,j) * ao_expo_pw_ord_transp(4,n,j) + KA2 = ao_expo_pw_ord_transp(4,n,j) do l = 1, ao_prim_num(i) beta = ao_expo_cgtos_ord_transp(l,i) beta_inv = (1.d0, 0.d0) / beta - do m = 1, 3 B_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_expo_pw_ord_transp(m,l,i) enddo phiB = ao_expo_phase_ord_transp(4,l,i) - KB2 = ao_expo_pw_ord_transp(4,l,i) * ao_expo_pw_ord_transp(4,l,i) + KB2 = ao_expo_pw_ord_transp(4,l,i) c = ao_coef_cgtos_norm_ord_transp(n,j) * ao_coef_cgtos_norm_ord_transp(l,i) - 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 * (conjg(alpha_inv) * KA2 + beta_inv * KB2)) + 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(A_center, B_center, alpha, beta, power_A, power_B, & overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1) - call overlap_cgaussian_xyz(A_center, B_center, alpha, conjg(beta), power_A, power_B, & + call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & overlap_x0_2, overlap_y0_2, overlap_z0_2, overlap, dim1) ! --- @@ -106,7 +106,7 @@ call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & overlap_m2_1, overlap_y, overlap_z, overlap, dim1) - call overlap_cgaussian_xyz(A_center, B_center, alpha, conjg(beta), power_A, power_B, & + call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & overlap_m2_2, overlap_y, overlap_z, overlap, dim1) else overlap_m2_1 = (0.d0, 0.d0) @@ -117,7 +117,7 @@ call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & overlap_p2_1, overlap_y, overlap_z, overlap, dim1) - call overlap_cgaussian_xyz(A_center, B_center, alpha, conjg(beta), power_A, power_B, & + call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & overlap_p2_2, overlap_y, overlap_z, overlap, dim1) power_A(1) = power_A(1) - 2 @@ -141,7 +141,7 @@ call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & overlap_x, overlap_m2_1, overlap_y, overlap, dim1) - call overlap_cgaussian_xyz(A_center, B_center, alpha, conjg(beta), power_A, power_B, & + call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & overlap_x, overlap_m2_2, overlap_y, overlap, dim1) else overlap_m2_1 = (0.d0, 0.d0) @@ -152,7 +152,7 @@ call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & overlap_x, overlap_p2_1, overlap_y, overlap, dim1) - call overlap_cgaussian_xyz(A_center, B_center, alpha, conjg(beta), power_A, power_B, & + call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & overlap_x, overlap_p2_2, overlap_y, overlap, dim1) power_A(2) = power_A(2) - 2 @@ -176,7 +176,7 @@ call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & overlap_x, overlap_y, overlap_m2_1, overlap, dim1) - call overlap_cgaussian_xyz(A_center, B_center, alpha, conjg(beta), power_A, power_B, & + call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & overlap_x, overlap_y, overlap_m2_2, overlap, dim1) else overlap_m2_1 = (0.d0, 0.d0) @@ -187,7 +187,7 @@ call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & overlap_x, overlap_y, overlap_p2_1, overlap, dim1) - call overlap_cgaussian_xyz(A_center, B_center, alpha, conjg(beta), power_A, power_B, & + call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & overlap_x, overlap_y, overlap_p2_2, overlap, dim1) power_A(3) = power_A(3) - 2 @@ -227,11 +227,12 @@ BEGIN_PROVIDER [double precision, ao_kinetic_integrals_cgtos, (ao_num, ao_num)] END_DOC implicit none + integer :: i, j - !$OMP PARALLEL DO DEFAULT(NONE) & - !$OMP PRIVATE(i, j) & - !$OMP SHARED(ao_num, ao_kinetic_integrals_cgtos, ao_deriv2_cgtos_x, ao_deriv2_cgtos_y, ao_deriv2_cgtos_z) + !$OMP PARALLEL DO DEFAULT(NONE) & + !$OMP PRIVATE(i, j) & + !$OMP SHARED(ao_num, ao_kinetic_integrals_cgtos, ao_deriv2_cgtos_x, ao_deriv2_cgtos_y, ao_deriv2_cgtos_z) do j = 1, ao_num do i = 1, ao_num ao_kinetic_integrals_cgtos(i,j) = -0.5d0 * (ao_deriv2_cgtos_x(i,j) + & @@ -239,8 +240,9 @@ BEGIN_PROVIDER [double precision, ao_kinetic_integrals_cgtos, (ao_num, ao_num)] ao_deriv2_cgtos_z(i,j)) enddo enddo - !$OMP END PARALLEL DO + !$OMP END PARALLEL DO END_PROVIDER ! --- + diff --git a/src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f b/src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f index 0d95ced4..3ac6f1a1 100644 --- a/src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f +++ b/src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f @@ -17,10 +17,14 @@ double precision function ao_two_e_integral_cgtos(i, j, k, l) integer :: ii, jj, kk, ll, dim1, I_power(3), J_power(3), K_power(3), L_power(3) integer :: iorder_p1(3), iorder_p2(3), iorder_q1(3), iorder_q2(3) double precision :: coef1, coef2, coef3, coef4 - complex*16 :: expo1, expo1_inv, I_center(3), KI2, phiI - complex*16 :: expo2, expo2_inv, J_center(3), KJ2, phiJ - complex*16 :: expo3, expo3_inv, K_center(3), KK2, phiK - complex*16 :: expo4, expo4_inv, L_center(3), KL2, phiL + double precision :: KI2, phiI + double precision :: KJ2, phiJ + double precision :: KK2, phiK + double precision :: KL2, phiL + complex*16 :: expo1, expo1_inv, I_center(3) + complex*16 :: expo2, expo2_inv, J_center(3) + complex*16 :: expo3, expo3_inv, K_center(3) + complex*16 :: expo4, expo4_inv, L_center(3) complex*16 :: P1_new(0:max_dim,3), P1_center(3), fact_p1, pp1, p1_inv complex*16 :: P2_new(0:max_dim,3), P2_center(3), fact_p2, pp2, p2_inv complex*16 :: Q1_new(0:max_dim,3), Q1_center(3), fact_q1, qq1, q1_inv @@ -70,7 +74,7 @@ double precision function ao_two_e_integral_cgtos(i, j, k, l) I_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,p,i) enddo phiI = ao_expo_phase_ord_transp(4,p,i) - KI2 = ao_expo_pw_ord_transp(4,p,i) * ao_expo_pw_ord_transp(4,p,i) + KI2 = ao_expo_pw_ord_transp(4,p,i) do q = 1, ao_prim_num(j) @@ -81,7 +85,7 @@ double precision function ao_two_e_integral_cgtos(i, j, k, l) J_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,q,j) enddo phiJ = ao_expo_phase_ord_transp(4,q,j) - KJ2 = ao_expo_pw_ord_transp(4,q,j) * ao_expo_pw_ord_transp(4,q,j) + KJ2 = ao_expo_pw_ord_transp(4,q,j) call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, & expo1, expo2, I_power, J_power, I_center, J_center, dim1) @@ -100,7 +104,7 @@ double precision function ao_two_e_integral_cgtos(i, j, k, l) K_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo3_inv * ao_expo_pw_ord_transp(m,r,k) enddo phiK = ao_expo_phase_ord_transp(4,r,k) - KK2 = ao_expo_pw_ord_transp(4,r,k) * ao_expo_pw_ord_transp(4,r,k) + KK2 = ao_expo_pw_ord_transp(4,r,k) do s = 1, ao_prim_num(l) @@ -111,7 +115,7 @@ double precision function ao_two_e_integral_cgtos(i, j, k, l) L_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo4_inv * ao_expo_pw_ord_transp(m,s,l) enddo phiL = ao_expo_phase_ord_transp(4,s,l) - KL2 = ao_expo_pw_ord_transp(4,s,l) * ao_expo_pw_ord_transp(4,s,l) + KL2 = ao_expo_pw_ord_transp(4,s,l) call give_explicit_cpoly_and_cgaussian(Q1_new, Q1_center, qq1, fact_q1, iorder_q1, & expo3, expo4, K_power, L_power, K_center, L_center, dim1) @@ -189,7 +193,7 @@ double precision function ao_two_e_integral_cgtos(i, j, k, l) I_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,p,i) enddo phiI = ao_expo_phase_ord_transp(4,p,i) - KI2 = ao_expo_pw_ord_transp(4,p,i) * ao_expo_pw_ord_transp(4,p,i) + KI2 = ao_expo_pw_ord_transp(4,p,i) do q = 1, ao_prim_num(j) @@ -200,7 +204,7 @@ double precision function ao_two_e_integral_cgtos(i, j, k, l) J_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,q,j) enddo phiJ = ao_expo_phase_ord_transp(4,q,j) - KJ2 = ao_expo_pw_ord_transp(4,q,j) * ao_expo_pw_ord_transp(4,q,j) + KJ2 = ao_expo_pw_ord_transp(4,q,j) do r = 1, ao_prim_num(k) @@ -211,7 +215,7 @@ double precision function ao_two_e_integral_cgtos(i, j, k, l) K_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo3_inv * ao_expo_pw_ord_transp(m,r,k) enddo phiK = ao_expo_phase_ord_transp(4,r,k) - KK2 = ao_expo_pw_ord_transp(4,r,k) * ao_expo_pw_ord_transp(4,r,k) + KK2 = ao_expo_pw_ord_transp(4,r,k) do s = 1, ao_prim_num(l) @@ -222,7 +226,7 @@ double precision function ao_two_e_integral_cgtos(i, j, k, l) L_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo4_inv * ao_expo_pw_ord_transp(m,s,l) enddo phiL = ao_expo_phase_ord_transp(4,s,l) - KL2 = ao_expo_pw_ord_transp(4,s,l) * ao_expo_pw_ord_transp(4,s,l) + KL2 = ao_expo_pw_ord_transp(4,s,l) C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL) & - 0.25d0 * (expo1_inv * KI2 + expo2_inv * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) @@ -313,10 +317,14 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) integer :: ii, jj, kk, ll, dim1, I_power(3), J_power(3), K_power(3), L_power(3) integer :: iorder_p1(3), iorder_p2(3), iorder_q1(3), iorder_q2(3) double precision :: coef1, coef2, coef3, coef4 - complex*16 :: expo1, expo1_inv, I_center(3), KI2, phiI - complex*16 :: expo2, expo2_inv, J_center(3), KJ2, phiJ - complex*16 :: expo3, expo3_inv, K_center(3), KK2, phiK - complex*16 :: expo4, expo4_inv, L_center(3), KL2, phiL + double precision :: KI2, phiI + double precision :: KJ2, phiJ + double precision :: KK2, phiK + double precision :: KL2, phiL + complex*16 :: expo1, expo1_inv, I_center(3) + complex*16 :: expo2, expo2_inv, J_center(3) + complex*16 :: expo3, expo3_inv, K_center(3) + complex*16 :: expo4, expo4_inv, L_center(3) complex*16 :: P1_new(0:max_dim,3), P1_center(3), fact_p1, pp1, p1_inv complex*16 :: P2_new(0:max_dim,3), P2_center(3), fact_p2, pp2, p2_inv complex*16 :: Q1_new(0:max_dim,3), Q1_center(3), fact_q1, qq1, q1_inv @@ -366,7 +374,7 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) K_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,r,k) enddo phiK = ao_expo_phase_ord_transp(4,r,k) - KK2 = ao_expo_pw_ord_transp(4,r,k) * ao_expo_pw_ord_transp(4,r,k) + KK2 = ao_expo_pw_ord_transp(4,r,k) schwartz_kl(0,r) = 0.d0 do s = 1, ao_prim_num(l) @@ -378,7 +386,7 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) L_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,s,l) enddo phiL = ao_expo_phase_ord_transp(4,s,l) - KL2 = ao_expo_pw_ord_transp(4,s,l) * ao_expo_pw_ord_transp(4,s,l) + KL2 = ao_expo_pw_ord_transp(4,s,l) call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, & expo1, expo2, K_power, L_power, K_center, L_center, dim1) @@ -393,7 +401,7 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) !C3 = C2 C4 = zexp((0.d0, 2.d0) * (phiK - phiL) - 0.5d0 * (conjg(expo1_inv) * KK2 + expo2_inv * KL2)) C5 = zexp(-(0.d0, 2.d0) * phiK - 0.5d0 * (expo1_inv * KK2 + real(expo2_inv) * KL2)) - C6 = zexp(-0.5d0 * (real(expo1_inv) * KK2 + real(expo2_inv) * KL2)) + C6 = zexp(-(0.5d0, 0.d0) * (real(expo1_inv) * KK2 + real(expo2_inv) * KL2)) !C7 = C6 !C8 = conjg(C5) @@ -452,7 +460,7 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) I_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,p,i) enddo phiI = ao_expo_phase_ord_transp(4,p,i) - KI2 = ao_expo_pw_ord_transp(4,p,i) * ao_expo_pw_ord_transp(4,p,i) + KI2 = ao_expo_pw_ord_transp(4,p,i) do q = 1, ao_prim_num(j) @@ -463,7 +471,7 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) J_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,q,j) enddo phiJ = ao_expo_phase_ord_transp(4,q,j) - KJ2 = ao_expo_pw_ord_transp(4,q,j) * ao_expo_pw_ord_transp(4,q,j) + KJ2 = ao_expo_pw_ord_transp(4,q,j) call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, & expo1, expo2, I_power, J_power, I_center, J_center, dim1) @@ -478,7 +486,7 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) !C3 = C2 C4 = zexp((0.d0, 2.d0) * (phiI - phiJ) - 0.5d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2)) C5 = zexp(-(0.d0, 2.d0) * phiI - 0.5d0 * (expo1_inv * KI2 + real(expo2_inv) * KJ2)) - C6 = zexp(-0.5d0 * (real(expo1_inv) * KI2 + real(expo2_inv) * KJ2)) + C6 = zexp(-(0.5d0, 0.d0) * (real(expo1_inv) * KI2 + real(expo2_inv) * KJ2)) !C7 = C6 !C8 = conjg(C5) @@ -533,7 +541,7 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) K_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo3_inv * ao_expo_pw_ord_transp(m,r,k) enddo phiK = ao_expo_phase_ord_transp(4,r,k) - KK2 = ao_expo_pw_ord_transp(4,r,k) * ao_expo_pw_ord_transp(4,r,k) + KK2 = ao_expo_pw_ord_transp(4,r,k) do s = 1, ao_prim_num(l) if(schwartz_kl(s,r)*schwartz_ij < thr) cycle @@ -545,7 +553,7 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) L_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo4_inv * ao_expo_pw_ord_transp(m,s,l) enddo phiL = ao_expo_phase_ord_transp(4,s,l) - KL2 = ao_expo_pw_ord_transp(4,s,l) * ao_expo_pw_ord_transp(4,s,l) + KL2 = ao_expo_pw_ord_transp(4,s,l) call give_explicit_cpoly_and_cgaussian(Q1_new, Q1_center, qq1, fact_q1, iorder_q1, & expo3, expo4, K_power, L_power, K_center, L_center, dim1) @@ -624,7 +632,7 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) K_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,r,k) enddo phiK = ao_expo_phase_ord_transp(4,r,k) - KK2 = ao_expo_pw_ord_transp(4,r,k) * ao_expo_pw_ord_transp(4,r,k) + KK2 = ao_expo_pw_ord_transp(4,r,k) schwartz_kl(0,r) = 0.d0 do s = 1, ao_prim_num(l) @@ -636,14 +644,14 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) L_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,s,l) enddo phiL = ao_expo_phase_ord_transp(4,s,l) - KL2 = ao_expo_pw_ord_transp(4,s,l) * ao_expo_pw_ord_transp(4,s,l) + KL2 = ao_expo_pw_ord_transp(4,s,l) C1 = zexp(-(0.d0, 2.d0) * (phiK + phiL) - 0.5d0 * (expo1_inv * KK2 + expo2_inv * KL2)) C2 = zexp(-(0.d0, 2.d0) * phiL - 0.5d0 * (real(expo1_inv) * KK2 + expo2_inv * KL2)) !C3 = C2 C4 = zexp((0.d0, 2.d0) * (phiK - phiL) - 0.5d0 * (conjg(expo1_inv) * KK2 + expo2_inv * KL2)) C5 = zexp(-(0.d0, 2.d0) * phiK - 0.5d0 * (expo1_inv * KK2 + real(expo2_inv) * KL2)) - C6 = zexp(-0.5d0 * (real(expo1_inv) * KK2 + real(expo2_inv) * KL2)) + C6 = zexp(-(0.5d0, 0.d0) * (real(expo1_inv) * KK2 + real(expo2_inv) * KL2)) !C7 = C6 !C8 = conjg(C5) @@ -708,7 +716,7 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) I_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,p,i) enddo phiI = ao_expo_phase_ord_transp(4,p,i) - KI2 = ao_expo_pw_ord_transp(4,p,i) * ao_expo_pw_ord_transp(4,p,i) + KI2 = ao_expo_pw_ord_transp(4,p,i) do q = 1, ao_prim_num(j) @@ -719,14 +727,14 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) J_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,q,j) enddo phiJ = ao_expo_phase_ord_transp(4,q,j) - KJ2 = ao_expo_pw_ord_transp(4,q,j) * ao_expo_pw_ord_transp(4,q,j) + KJ2 = ao_expo_pw_ord_transp(4,q,j) C1 = zexp(-(0.d0, 2.d0) * (phiI + phiJ) - 0.5d0 * (expo1_inv * KI2 + expo2_inv * KJ2)) C2 = zexp(-(0.d0, 2.d0) * phiJ - 0.5d0 * (real(expo1_inv) * KI2 + expo2_inv * KJ2)) !C3 = C2 C4 = zexp((0.d0, 2.d0) * (phiI - phiJ) - 0.5d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2)) C5 = zexp(-(0.d0, 2.d0) * phiI - 0.5d0 * (expo1_inv * KI2 + real(expo2_inv) * KJ2)) - C6 = zexp(-0.5d0 * (real(expo1_inv) * KI2 + real(expo2_inv) * KJ2)) + C6 = zexp(-(0.5d0, 0.d0) * (real(expo1_inv) * KI2 + real(expo2_inv) * KJ2)) !C7 = C6 !C8 = conjg(C5) @@ -788,7 +796,7 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) K_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo3_inv * ao_expo_pw_ord_transp(m,r,k) enddo phiK = ao_expo_phase_ord_transp(4,r,k) - KK2 = ao_expo_pw_ord_transp(4,r,k) * ao_expo_pw_ord_transp(4,r,k) + KK2 = ao_expo_pw_ord_transp(4,r,k) do s = 1, ao_prim_num(l) if(schwartz_kl(s,r)*schwartz_ij < thr) cycle @@ -800,7 +808,7 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) L_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo4_inv * ao_expo_pw_ord_transp(m,s,l) enddo phiL = ao_expo_phase_ord_transp(4,s,l) - KL2 = ao_expo_pw_ord_transp(4,s,l) * ao_expo_pw_ord_transp(4,s,l) + KL2 = ao_expo_pw_ord_transp(4,s,l) C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL) & - 0.25d0 * (expo1_inv * KI2 + expo2_inv * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) diff --git a/src/utils/cgtos_one_e.irp.f b/src/utils/cgtos_one_e.irp.f index 43ca8224..80b82eaf 100644 --- a/src/utils/cgtos_one_e.irp.f +++ b/src/utils/cgtos_one_e.irp.f @@ -42,12 +42,12 @@ complex*16 function overlap_cgaussian_x(A_center, B_center, alpha, beta, power_A overlap_cgaussian_x *= fact_p -end function overlap_cgaussian_x +end ! --- -subroutine overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B & - , overlap_x, overlap_y, overlap_z, overlap, dim ) +subroutine overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & + overlap_x, overlap_y, overlap_z, overlap, dim ) BEGIN_DOC ! @@ -113,7 +113,7 @@ subroutine overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, powe overlap = overlap_x * overlap_y * overlap_z -end subroutine overlap_cgaussian_xyz +end ! ---