From c1bdfe7e93489d03adf2d69b2d8d89ee5e72b49e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 20 Nov 2022 14:52:27 +0100 Subject: [PATCH] Starting to transpose --- src/ao_many_one_e_ints/ao_gaus_gauss.irp.f | 60 +++++++++++++++++-- src/utils/one_e_integration.irp.f | 67 +++++++++++++++++++++- 2 files changed, 119 insertions(+), 8 deletions(-) diff --git a/src/ao_many_one_e_ints/ao_gaus_gauss.irp.f b/src/ao_many_one_e_ints/ao_gaus_gauss.irp.f index 19cfe4cd..a657d6b1 100644 --- a/src/ao_many_one_e_ints/ao_gaus_gauss.irp.f +++ b/src/ao_many_one_e_ints/ao_gaus_gauss.irp.f @@ -150,6 +150,58 @@ double precision function overlap_gauss_r12_ao(D_center, delta, i, j) end function overlap_gauss_r12_ao +! -- + +subroutine overlap_gauss_r12_ao_v(D_center, delta, i, j, resv, n_points) + + BEGIN_DOC + ! \int dr AO_i(r) AO_j(r) e^{-delta |r-D_center|^2} + END_DOC + + implicit none + integer, intent(in) :: i, j, n_points + double precision, intent(in) :: D_center(3,n_points), delta + double precision, intent(out) :: resv(n_points) + + integer :: power_A(3), power_B(3), l, k + double precision :: A_center(3), B_center(3), alpha, beta, coef, coef1 + double precision, allocatable :: analytical_j(:) + + double precision, external :: overlap_gauss_r12 + integer :: ipoint + + resv(:) = 0.d0 + if(ao_overlap_abs(j,i).lt.1.d-12) then + return + endif + + power_A(1:3) = ao_power(i,1:3) + power_B(1:3) = ao_power(j,1:3) + + A_center(1:3) = nucl_coord(ao_nucl(i),1:3) + B_center(1:3) = nucl_coord(ao_nucl(j),1:3) + + allocate(analytical_j(n_points)) + do l = 1, ao_prim_num(i) + alpha = ao_expo_ordered_transp (l,i) + coef1 = ao_coef_normalized_ordered_transp(l,i) + + do k = 1, ao_prim_num(j) + beta = ao_expo_ordered_transp(k,j) + coef = coef1 * ao_coef_normalized_ordered_transp(k,j) + + if(dabs(coef) .lt. 1d-12) cycle + + call overlap_gauss_r12_v(D_center, delta, A_center, B_center, power_A, power_B, alpha, beta, analytical_j, n_points) + do ipoint=1, n_points + resv(ipoint) = resv(ipoint) + coef*analytical_j(ipoint) + enddo + enddo + enddo + deallocate(analytical_j) + +end + ! --- double precision function overlap_gauss_r12_ao_with1s(B_center, beta, D_center, delta, i, j) @@ -241,9 +293,6 @@ subroutine overlap_gauss_r12_ao_with1s_v(B_center, beta, D_center, delta, i, j, double precision :: gama, gama_inv double precision :: bg, dg, bdg - double precision, external :: overlap_gauss_r12, overlap_gauss_r12_ao - - double precision, external :: overlap_gauss_r12_ao_with1s integer :: ipoint double precision, allocatable :: fact_g(:), G_center(:,:), analytical_j(:) @@ -255,9 +304,7 @@ subroutine overlap_gauss_r12_ao_with1s_v(B_center, beta, D_center, delta, i, j, ASSERT(beta .gt. 0.d0) if(beta .lt. 1d-10) then - do ipoint=1,n_points - resv(ipoint) = overlap_gauss_r12_ao(D_center(1,ipoint), delta, i, j) - enddo + call overlap_gauss_r12_ao_v(D_center, delta, i, j, resv, n_points) return endif @@ -317,6 +364,7 @@ subroutine overlap_gauss_r12_ao_with1s_v(B_center, beta, D_center, delta, i, j, enddo enddo + deallocate (fact_g, G_center, analytical_j ) end diff --git a/src/utils/one_e_integration.irp.f b/src/utils/one_e_integration.irp.f index 9c1d2445..e175042b 100644 --- a/src/utils/one_e_integration.irp.f +++ b/src/utils/one_e_integration.irp.f @@ -92,9 +92,9 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,& overlap = overlap_x * overlap_y * overlap_z end - + ! --- - + subroutine overlap_x_abs(A_center, B_center, alpha, beta, power_A, power_B, overlap_x, lower_exp_val, dx, nx) BEGIN_DOC @@ -151,4 +151,67 @@ subroutine overlap_x_abs(A_center, B_center, alpha, beta, power_A, power_B, over end +! --- +subroutine overlap_gaussian_xyz_v(A_center,B_center,alpha,beta,power_A,& + power_B,overlap_x,overlap_y,overlap_z,overlap,dim, n_points) + implicit none + BEGIN_DOC + !.. math:: + ! + ! S_x = \int (x-A_x)^{a_x} exp(-\alpha(x-A_x)^2) (x-B_x)^{b_x} exp(-beta(x-B_x)^2) dx \\ + ! S = S_x S_y S_z + ! + END_DOC + include 'constants.include.F' + integer,intent(in) :: dim, n_points + double precision,intent(in) :: A_center(3),B_center(3) ! center of the x1 functions + double precision, intent(in) :: alpha,beta + integer,intent(in) :: power_A(3), power_B(3) ! power of the x1 functions + double precision, intent(out) :: overlap_x(n_points),overlap_y(n_points),overlap_z(n_points),overlap(n_points) + double precision :: P_new(0:max_dim,3),P_center(3),fact_p,p + double precision :: F_integral_tab(0:max_dim) + integer :: iorder_p(3) + + call give_explicit_poly_and_gaussian(P_new,P_center,p,fact_p,iorder_p,alpha,beta,power_A,power_B,A_center,B_center,dim) + if(fact_p.lt.1d-20)then + overlap_x = 1.d-10 + overlap_y = 1.d-10 + overlap_z = 1.d-10 + overlap = 1.d-10 + return + endif + integer :: nmax + double precision :: F_integral + nmax = maxval(iorder_p) + do i = 0,nmax + F_integral_tab(i) = F_integral(i,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) + + integer :: i + do i = 1,iorder_p(1) + overlap_x = overlap_x + P_new(i,1) * F_integral_tab(i) + enddo + call gaussian_product_x(alpha,A_center(1),beta,B_center(1),fact_p,p,P_center(1)) + overlap_x *= fact_p + + do i = 1,iorder_p(2) + overlap_y = overlap_y + P_new(i,2) * F_integral_tab(i) + enddo + call gaussian_product_x(alpha,A_center(2),beta,B_center(2),fact_p,p,P_center(2)) + overlap_y *= fact_p + + do i = 1,iorder_p(3) + overlap_z = overlap_z + P_new(i,3) * F_integral_tab(i) + enddo + call gaussian_product_x(alpha,A_center(3),beta,B_center(3),fact_p,p,P_center(3)) + overlap_z *= fact_p + + overlap = overlap_x * overlap_y * overlap_z + +end + +! ---