From 9407a6c7a98a8efbb78503b9f9d5404e1b7e67ab Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Mon, 14 Oct 2024 16:36:23 +0200 Subject: [PATCH] added Kinetic integ for general cGTOs --- src/ao_one_e_ints/kin_ao_ints.irp.f | 13 +- .../one_e_kin_integrals_cgtos.irp.f | 169 ++++++++++-------- 2 files changed, 102 insertions(+), 80 deletions(-) diff --git a/src/ao_one_e_ints/kin_ao_ints.irp.f b/src/ao_one_e_ints/kin_ao_ints.irp.f index ea5c2670..49eb56ad 100644 --- a/src/ao_one_e_ints/kin_ao_ints.irp.f +++ b/src/ao_one_e_ints/kin_ao_ints.irp.f @@ -24,7 +24,6 @@ double precision :: d_a_2,d_2 if(use_cgtos) then - !print*, 'use_cgtos for ao_kinetic_integrals ?', use_cgtos do j = 1, ao_num do i = 1, ao_num @@ -92,8 +91,8 @@ power_A(1) = power_A(1)-2 double precision :: deriv_tmp - deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(1) +1.d0) * overlap_x0 & - +power_A(1) * (power_A(1)-1.d0) * d_a_2 & + deriv_tmp = (-2.d0 * alpha * (2.d0 * dble(power_A(1)) +1.d0) * overlap_x0 & + +dble(power_A(1)) * (dble(power_A(1))-1.d0) * d_a_2 & +4.d0 * alpha * alpha * d_2 )*overlap_y0*overlap_z0 ao_deriv2_x(i,j) += c*deriv_tmp @@ -107,8 +106,8 @@ call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,d_2,overlap_z,overlap,dim1) power_A(2) = power_A(2)-2 - deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(2) +1.d0 ) * overlap_y0 & - +power_A(2) * (power_A(2)-1.d0) * d_a_2 & + deriv_tmp = (-2.d0 * alpha * (2.d0 * dble(power_A(2)) +1.d0 ) * overlap_y0 & + +dble(power_A(2)) * (dble(power_A(2))-1.d0) * d_a_2 & +4.d0 * alpha * alpha * d_2 )*overlap_x0*overlap_z0 ao_deriv2_y(i,j) += c*deriv_tmp @@ -122,8 +121,8 @@ call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,overlap_z,d_2,overlap,dim1) power_A(3) = power_A(3)-2 - deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(3) +1.d0 ) * overlap_z0 & - +power_A(3) * (power_A(3)-1.d0) * d_a_2 & + deriv_tmp = (-2.d0 * alpha * (2.d0 * dble(power_A(3)) +1.d0 ) * overlap_z0 & + +dble(power_A(3)) * (dble(power_A(3))-1.d0) * d_a_2 & +4.d0 * alpha * alpha * d_2 )*overlap_x0*overlap_y0 ao_deriv2_z(i,j) += c*deriv_tmp 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 2d4d0d17..96ac308f 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 @@ -1,14 +1,15 @@ ! --- - BEGIN_PROVIDER [ double precision, ao_deriv2_cgtos_x, (ao_num, ao_num) ] -&BEGIN_PROVIDER [ double precision, ao_deriv2_cgtos_y, (ao_num, ao_num) ] -&BEGIN_PROVIDER [ double precision, ao_deriv2_cgtos_z, (ao_num, ao_num) ] + BEGIN_PROVIDER [double precision, ao_deriv2_cgtos_x, (ao_num, ao_num)] +&BEGIN_PROVIDER [double precision, ao_deriv2_cgtos_y, (ao_num, ao_num)] +&BEGIN_PROVIDER [double precision, ao_deriv2_cgtos_z, (ao_num, ao_num)] implicit none - integer :: i, j, n, l, dim1, power_A(3), power_B(3) + integer :: i, j, m, n, l, ii, jj, dim1, power_A(3), power_B(3) double precision :: c, deriv_tmp - complex*16 :: alpha, beta, A_center(3), B_center(3) + complex*16 :: alpha, alpha_inv, A_center(3), KA2, phiA, C1 + complex*16 :: beta, beta_inv, B_center(3), KB2, phiB, 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 @@ -27,87 +28,109 @@ beta = (0.1d0, 0.d0) power_A = 1 power_B = 0 - 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, beta, power_A, power_B, & + overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1) ! --- - !$OMP PARALLEL DO SCHEDULE(GUIDED) & - !$OMP DEFAULT(NONE) & - !$OMP PRIVATE( A_center, B_center, power_A, power_B, alpha, beta, i, j, l, n, c & - !$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, 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_deriv2_cgtos_x, ao_deriv2_cgtos_y, ao_deriv2_cgtos_z ) + !$OMP PARALLEL DO SCHEDULE(GUIDED) & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(i, j, m, n, l, ii, jj, c, C1, C2, & + !$OMP A_center, power_A, alpha, alpha_inv, KA2, phiA, & + !$OMP B_center, power_B, beta, beta_inv, KB2, phiB, & + !$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) do j = 1, ao_num - A_center(1) = nucl_coord(ao_nucl(j),1) * (1.d0, 0.d0) - A_center(2) = nucl_coord(ao_nucl(j),2) * (1.d0, 0.d0) - A_center(3) = nucl_coord(ao_nucl(j),3) * (1.d0, 0.d0) - power_A(1) = ao_power(j,1) - power_A(2) = ao_power(j,2) - power_A(3) = ao_power(j,3) + + jj = ao_nucl(j) + power_A(1) = ao_power(j,1) + power_A(2) = ao_power(j,2) + power_A(3) = ao_power(j,3) do i = 1, ao_num - B_center(1) = nucl_coord(ao_nucl(i),1) * (1.d0, 0.d0) - B_center(2) = nucl_coord(ao_nucl(i),2) * (1.d0, 0.d0) - B_center(3) = nucl_coord(ao_nucl(i),3) * (1.d0, 0.d0) - power_B(1) = ao_power(i,1) - power_B(2) = ao_power(i,2) - power_B(3) = ao_power(i,3) + + ii = ao_nucl(i) + power_B(1) = ao_power(i,1) + power_B(2) = ao_power(i,2) + power_B(3) = ao_power(i,3) ao_deriv2_cgtos_x(i,j) = 0.d0 ao_deriv2_cgtos_y(i,j) = 0.d0 ao_deriv2_cgtos_z(i,j) = 0.d0 do n = 1, ao_prim_num(j) + 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) do l = 1, ao_prim_num(i) - c = ao_coef_cgtos_norm_ord_transp(n,j) * ao_coef_cgtos_norm_ord_transp(l,i) + beta = ao_expo_cgtos_ord_transp(l,i) + beta_inv = (1.d0, 0.d0) / beta - 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 ) + 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) - call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B & - , overlap_x0_2, overlap_y0_2, overlap_z0_2, overlap, dim1 ) + 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)) + + 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, & + overlap_x0_2, overlap_y0_2, overlap_z0_2, overlap, dim1) ! --- power_A(1) = power_A(1) - 2 if(power_A(1) > -1) then - 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, 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 & - , overlap_m2_2, overlap_y, overlap_z, overlap, dim1 ) + call overlap_cgaussian_xyz(A_center, 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) overlap_m2_2 = (0.d0, 0.d0) endif power_A(1) = power_A(1) + 4 - 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, 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 & - , overlap_p2_2, overlap_y, overlap_z, overlap, dim1 ) + call overlap_cgaussian_xyz(A_center, 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 - deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * power_A(1) + 1.d0) * overlap_x0_1 & - + power_A(1) * (power_A(1) - 1.d0) * overlap_m2_1 & + 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 - deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * power_A(1) + 1.d0) * overlap_x0_2 & - + power_A(1) * (power_A(1) - 1.d0) * overlap_m2_2 & + 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 = 2.d0 * real(deriv_tmp_1 + deriv_tmp_2) + deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2) ao_deriv2_cgtos_x(i,j) += c * deriv_tmp @@ -115,34 +138,34 @@ power_A(2) = power_A(2) - 2 if(power_A(2) > -1) then - 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, 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 & - , overlap_x, overlap_m2_2, overlap_y, overlap, dim1 ) + call overlap_cgaussian_xyz(A_center, 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) overlap_m2_2 = (0.d0, 0.d0) endif power_A(2) = power_A(2) + 4 - 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, 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 & - , overlap_x, overlap_p2_2, overlap_y, overlap, dim1 ) + call overlap_cgaussian_xyz(A_center, 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 - deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * power_A(2) + 1.d0) * overlap_y0_1 & - + power_A(2) * (power_A(2) - 1.d0) * overlap_m2_1 & + 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 - deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * power_A(2) + 1.d0) * overlap_y0_2 & - + power_A(2) * (power_A(2) - 1.d0) * overlap_m2_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 = 2.d0 * real(deriv_tmp_1 + deriv_tmp_2) + deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2) ao_deriv2_cgtos_y(i,j) += c * deriv_tmp @@ -150,34 +173,34 @@ power_A(3) = power_A(3) - 2 if(power_A(3) > -1) then - 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, 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 & - , overlap_x, overlap_y, overlap_m2_2, overlap, dim1 ) + call overlap_cgaussian_xyz(A_center, 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) overlap_m2_2 = (0.d0, 0.d0) endif power_A(3) = power_A(3) + 4 - 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, 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 & - , overlap_x, overlap_y, overlap_p2_2, overlap, dim1 ) + call overlap_cgaussian_xyz(A_center, 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 - deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * power_A(3) + 1.d0) * overlap_z0_1 & - + power_A(3) * (power_A(3) - 1.d0) * overlap_m2_1 & + 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 * power_A(3) + 1.d0) * overlap_z0_2 & - + power_A(3) * (power_A(3) - 1.d0) * overlap_m2_2 & + 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 - deriv_tmp = 2.d0 * real(deriv_tmp_1 + deriv_tmp_2) + deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2) ao_deriv2_cgtos_z(i,j) += c * deriv_tmp