From b37256e6c616861475e45c0c85217a6b192fada0 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Sat, 24 Aug 2024 20:44:47 +0200 Subject: [PATCH] added & tested torus overlap --- src/ao_one_e_ints/ao_overlap.irp.f | 1 + src/ao_one_e_ints/aos_torus.irp.f | 96 ++++--- src/utils/torus_integration.irp.f | 365 ++++++++++++++++++++++++ src/utils/torus_one_e_integration.irp.f | 78 +++++ 4 files changed, 495 insertions(+), 45 deletions(-) create mode 100644 src/utils/torus_integration.irp.f create mode 100644 src/utils/torus_one_e_integration.irp.f diff --git a/src/ao_one_e_ints/ao_overlap.irp.f b/src/ao_one_e_ints/ao_overlap.irp.f index 0622efd6..a0c50521 100644 --- a/src/ao_one_e_ints/ao_overlap.irp.f +++ b/src/ao_one_e_ints/ao_overlap.irp.f @@ -43,6 +43,7 @@ enddo elseif(do_torus) then + !print*, ' do_torus for ao_overlap ?', do_torus do j = 1, ao_num do i = 1, ao_num diff --git a/src/ao_one_e_ints/aos_torus.irp.f b/src/ao_one_e_ints/aos_torus.irp.f index 96918154..409a5863 100644 --- a/src/ao_one_e_ints/aos_torus.irp.f +++ b/src/ao_one_e_ints/aos_torus.irp.f @@ -30,51 +30,57 @@ ao_overlap_torus_y = 0.d0 ao_overlap_torus_z = 0.d0 -! dim1=100 -! !$OMP PARALLEL DO SCHEDULE(GUIDED) & -! !$OMP DEFAULT(NONE) & -! !$OMP PRIVATE(A_center,B_center,power_A,power_B,& -! !$OMP overlap_x,overlap_y, overlap_z, overlap, & -! !$OMP alpha, beta,i,j,c) & -! !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & -! !$OMP ao_overlap_x,ao_overlap_y,ao_overlap_z,ao_overlap,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, & -! !$OMP ao_expo_ordered_transp,dim1) -! do j=1,ao_num -! A_center(1) = nucl_coord( ao_nucl(j), 1 ) -! A_center(2) = nucl_coord( ao_nucl(j), 2 ) -! A_center(3) = nucl_coord( ao_nucl(j), 3 ) -! 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 ) -! B_center(2) = nucl_coord( ao_nucl(i), 2 ) -! B_center(3) = nucl_coord( ao_nucl(i), 3 ) -! power_B(1) = ao_power( i, 1 ) -! power_B(2) = ao_power( i, 2 ) -! power_B(3) = ao_power( i, 3 ) -! do n = 1,ao_prim_num(j) -! alpha = ao_expo_ordered_transp(n,j) -! do l = 1, ao_prim_num(i) -! beta = ao_expo_ordered_transp(l,i) -! call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) -! c = ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i) -! ao_overlap(i,j) += c * overlap -! if(isnan(ao_overlap(i,j)))then -! print*,'i,j',i,j -! print*,'l,n',l,n -! print*,'c,overlap',c,overlap -! print*,overlap_x,overlap_y,overlap_z -! stop -! endif -! ao_overlap_x(i,j) += c * overlap_x -! ao_overlap_y(i,j) += c * overlap_y -! ao_overlap_z(i,j) += c * overlap_z -! enddo -! enddo -! enddo -! enddo -! !$OMP END PARALLEL DO + dim1 = 100 + + !$OMP PARALLEL & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(i, j, n, l, A_center, B_center, power_A, power_B, & + !$OMP alpha, beta, c, overlap_x, overlap_y, overlap_z, overlap) & + !$OMP SHARED(nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1, & + !$OMP ao_coef_normalized_ordered_transp, ao_expo_ordered_transp, torus_length, & + !$OMP ao_overlap_torus_x, ao_overlap_torus_y, ao_overlap_torus_z, ao_overlap_torus) + !$OMP DO SCHEDULE(GUIDED) + do j = 1, ao_num + + A_center(1) = nucl_coord(ao_nucl(j),1) + A_center(2) = nucl_coord(ao_nucl(j),2) + A_center(3) = nucl_coord(ao_nucl(j),3) + + 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) + B_center(2) = nucl_coord(ao_nucl(i),2) + B_center(3) = nucl_coord(ao_nucl(i),3) + + power_B(1) = ao_power(i,1) + power_B(2) = ao_power(i,2) + power_B(3) = ao_power(i,3) + + do n = 1, ao_prim_num(j) + alpha = ao_expo_ordered_transp(n,j) + + do l = 1, ao_prim_num(i) + beta = ao_expo_ordered_transp(l,i) + + call overlap_gaussian_torus_xyz(A_center, B_center, alpha, beta, power_A, power_B, torus_length, & + overlap_x, overlap_y, overlap_z, overlap, dim1) + + c = ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i) + + ao_overlap_torus (i,j) += c * overlap + ao_overlap_torus_x(i,j) += c * overlap_x + ao_overlap_torus_y(i,j) += c * overlap_y + ao_overlap_torus_z(i,j) += c * overlap_z + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL END_PROVIDER diff --git a/src/utils/torus_integration.irp.f b/src/utils/torus_integration.irp.f new file mode 100644 index 00000000..5c01ee1b --- /dev/null +++ b/src/utils/torus_integration.irp.f @@ -0,0 +1,365 @@ + +! --- + +subroutine gaussian_product_torus_x(a, xa, b, xb, Lx, k, p, xp) + + BEGIN_DOC + ! + ! TODO + ! + ! Gaussian product in 1D. + ! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2} + ! + END_DOC + + implicit none + + double precision, intent(in) :: a, b ! Exponents + double precision, intent(in) :: xa, xb ! Centers + double precision, intent(in) :: Lx + double precision, intent(out) :: p ! New exponent + double precision, intent(out) :: xp ! New center + double precision, intent(out) :: k ! Constant + + double precision :: p_inv + double precision :: xab + + p = a + b + + p_inv = 1.d0 / p + + xab = xa - xb + if (dabs(xab) > 0.5d0*Lx) then + xab = Lx - dabs(xab) + end if + + k = a * b * p_inv * xab * xab + if (k > 400.d0) then + k = 0.d0 + xp = 0.d0 + return + endif + k = dexp(-k) + + ! xp = (a*xa + b*xb) * p_inv + + xab = xa - xb + if (dabs(xab) > 0.5d0*Lx) then + + if (xa > xb) then + xp = ( a * xa + b * (xb + Lx) ) * p_inv + elseif (xa < xb) then + xp = ( a * (xa + Lx) + b * xb ) * p_inv + else + xp = ( a * xa + b * xb ) * p_inv + end if + + else + + xp = ( a * xa + b * xb ) * p_inv + + end if + +end + +! --- + +subroutine gaussian_product_torus(a, xa, b, xb, torus_L, k, p, xp) + + BEGIN_DOC + ! + ! TODO + ! + ! Gaussian product in 1D. + ! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2} + ! + END_DOC + + implicit none + + double precision, intent(in) :: a, b ! Exponents + double precision, intent(in) :: xa(3), xb(3) ! Centers + double precision, intent(in) :: torus_L(3) + double precision, intent(out) :: p ! New exponent + double precision, intent(out) :: xp(3) ! New center + double precision, intent(out) :: k ! Constant + + double precision :: p_inv + double precision :: xab(3) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xab + + p = a + b + + p_inv = 1.d0 / p + + xab(1) = xa(1) - xb(1) + if (dabs(xab(1)) > 0.5d0*torus_L(1)) then + xab(1) = torus_L(1) - dabs(xab(1)) + end if + + xab(2) = xa(2) - xb(2) + if (dabs(xab(2)) > 0.5d0*torus_L(2)) then + xab(2) = torus_L(2) - dabs(xab(2)) + end if + + xab(3) = xa(3) - xb(3) + if (dabs(xab(3)) > 0.5d0*torus_L(3)) then + xab(3) = torus_L(3) - dabs(xab(3)) + end if + + k = a * b * p_inv * (xab(1)*xab(1) + xab(2)*xab(2) + xab(3)*xab(3)) + if (k > 40.d0) then + k = 0.d0 + xp(1) = 0.d0 + xp(2) = 0.d0 + xp(3) = 0.d0 + return + endif + k = dexp(-k) + +! xp(1) = (a * xa(1) + b * xb(1)) * p_inv + + xab(1) = xa(1) - xb(1) + if (dabs(xab(1)) > 0.5d0*torus_L(1)) then + + if (xa(1) > xb(1)) then + xp(1) = ( a * xa(1) + b * (xb(1) + torus_L(1)) ) * p_inv + elseif (xa(1) < xb(1)) then + xp(1) = ( a * (xa(1) + torus_L(1)) + b * xb(1) ) * p_inv + else + xp(1) = ( a * xa(1) + b * xb(1) ) * p_inv + end if + + else + + xp(1) = ( a * xa(1) + b * xb(1) ) * p_inv + + end if + + +! xp(2) = (a * xa(2) + b * xb(2)) * p_inv + + xab(2) = xa(2) - xb(2) + if (dabs(xab(2)) > 0.5d0*torus_L(2)) then + + if (xa(2) > xb(2)) then + xp(2) = ( a * xa(2) + b * (xb(2) + torus_L(2)) ) * p_inv + elseif (xa(2) < xb(2)) then + xp(2) = ( a * (xa(2) + torus_L(2)) + b * xb(2) ) * p_inv + else + xp(2) = ( a * xa(2) + b * xb(2) ) * p_inv + end if + + else + + xp(2) = ( a * xa(2) + b * xb(2) ) * p_inv + + end if + + +! xp(3) = (a * xa(3) + b * xb(3)) * p_inv + + xab(3) = xa(3) - xb(3) + if (dabs(xab(3)) > 0.5d0*torus_L(3)) then + + if (xa(3) > xb(3)) then + xp(3) = ( a * xa(3) + b * (xb(3) + torus_L(3)) ) * p_inv + elseif (xa(3) < xb(3)) then + xp(3) = ( a * (xa(3) + torus_L(3)) + b * xb(3) ) * p_inv + else + xp(3) = ( a * xa(3) + b * xb(3) ) * p_inv + end if + + else + + xp(3) = ( a * xa(3) + b * xb(3) ) * p_inv + + end if + + +end + +! --- + +subroutine give_explicit_poly_and_gaussian_torus(P_new, P_center, p, fact_k, iorder, alpha, beta, a, b, A_center, B_center, torus_L, dim) + + BEGIN_DOC + ! + ! TODO + ! + ! Transforms the product of + ! (x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) + ! into + ! fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 ) + ! * [ sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y ] exp (- p (y-P_center(2))^2 ) + ! * [ sum (l_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z ] exp (- p (z-P_center(3))^2 ) + ! + ! WARNING ::: IF fact_k is too smal then: + ! returns a "s" function centered in zero + ! with an inifinite exponent and a zero polynom coef + ! + END_DOC + + implicit none + + include 'constants.include.F' + + integer, intent(in) :: dim + integer, intent(in) :: a(3), b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1) + double precision, intent(in) :: alpha, beta ! exponents + double precision, intent(in) :: A_center(3) ! A center + double precision, intent(in) :: B_center (3) ! B center + double precision, intent(in) :: torus_L(3) + integer, intent(out) :: iorder(3) ! i_order(i) = order of the polynomials + double precision, intent(out) :: P_center(3) ! new center + double precision, intent(out) :: p ! new exponent + double precision, intent(out) :: fact_k ! constant factor + double precision, intent(out) :: P_new(0:max_dim,3) ! polynomial + + integer :: n_new, i + double precision :: P_a(0:max_dim,3), P_b(0:max_dim,3) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: P_a, P_b + + iorder(1) = 0 + iorder(2) = 0 + iorder(3) = 0 + + P_new(0,1) = 0.d0 + P_new(0,2) = 0.d0 + P_new(0,3) = 0.d0 + + !DIR$ FORCEINLINE + call gaussian_product_torus(alpha, A_center, beta, B_center, torus_L, fact_k, p, P_center) + if (fact_k < thresh) then + ! IF fact_k is too smal then: + ! returns a "s" function centered in zero + ! with an inifinite exponent and a zero polynom coef + P_center = 0.d0 + p = 1.d+15 + P_new = 0.d0 + iorder = 0 + fact_k = 0.d0 + return + endif + + !DIR$ FORCEINLINE + call recentered_poly2_torus(P_a(0,1), A_center(1), P_center(1), a(1), P_b(0,1), B_center(1), P_center(1), b(1), torus_L(1)) + iorder(1) = a(1) + b(1) + do i = 0, iorder(1) + P_new(i,1) = 0.d0 + enddo + n_new = 0 + !DIR$ FORCEINLINE + call multiply_poly(P_a(0,1), a(1), P_b(0,1), b(1), P_new(0,1), n_new) + + !DIR$ FORCEINLINE + call recentered_poly2_torus(P_a(0,2), A_center(2), P_center(2), a(2), P_b(0,2), B_center(2), P_center(2), b(2), torus_L(2)) + iorder(2) = a(2) + b(2) + do i = 0, iorder(2) + P_new(i,2) = 0.d0 + enddo + n_new = 0 + !DIR$ FORCEINLINE + call multiply_poly(P_a(0,2), a(2), P_b(0,2), b(2), P_new(0,2), n_new) + + !DIR$ FORCEINLINE + call recentered_poly2_torus(P_a(0,3), A_center(3), P_center(3), a(3), P_b(0,3), B_center(3), P_center(3), b(3), torus_L(3)) + iorder(3) = a(3) + b(3) + do i = 0, iorder(3) + P_new(i,3) = 0.d0 + enddo + n_new = 0 + !DIR$ FORCEINLINE + call multiply_poly(P_a(0,3), a(3), P_b(0,3), b(3), P_new(0,3), n_new) + +end + +! --- + +subroutine recentered_poly2_torus(P_new, x_A, x_P, a, P_new2, x_B, x_Q, b, Lx) + + BEGIN_DOC + ! + ! TODO + ! + ! Recenter two polynomials + ! + END_DOC + + implicit none + + integer, intent(in) :: a, b + double precision, intent(in) :: x_A, x_P, x_B, x_Q + double precision, intent(in) :: Lx + double precision, intent(out) :: P_new(0:a), P_new2(0:b) + + integer :: i, minab, maxab + double precision :: sign + double precision :: pows_a(-2:a+b+4), pows_b(-2:a+b+4) + double precision, external :: binom_func + + PROVIDE binom_transp + + if((a<0) .or. (b<0)) then + P_new = 0.d0 + P_new2 = 0.d0 + return + endif + + maxab = max(a, b) + minab = max(min(a, b), 0) + + pows_a(0) = 1.d0 + pows_b(0) = 1.d0 + + pows_a(1) = (x_P - x_A) + if(pows_a(1) >= 0.d0) then + sign = 1.d0 + else + sign = -1.d0 + endif + if (dabs(pows_a(1)) > 0.5d0*Lx) then + pows_a(1) = -1.d0 * sign * (Lx - dabs(pows_a(1))) + end if + + pows_b(1) = (x_Q - x_B) + if(pows_b(1) >= 0.d0) then + sign = 1.d0 + else + sign = -1.d0 + endif + if (dabs(pows_b(1)) > 0.5d0*Lx) then + pows_b(1) = -1.d0 * sign * (Lx - dabs(pows_b(1))) + end if + + do i = 2, maxab + pows_a(i) = pows_a(i-1) * pows_a(1) + pows_b(i) = pows_b(i-1) * pows_b(1) + enddo + + P_new (0) = pows_a(a) + P_new2(0) = pows_b(b) + + do i = 1, min(minab, 20) + P_new (i) = binom_transp(a-i,a) * pows_a(a-i) + P_new2(i) = binom_transp(b-i,b) * pows_b(b-i) + enddo + + do i = minab+1, min(a, 20) + P_new (i) = binom_transp(a-i,a) * pows_a(a-i) + enddo + do i = minab+1, min(b, 20) + P_new2(i) = binom_transp(b-i,b) * pows_b(b-i) + enddo + + do i = 101, a + P_new(i) = binom_func(a,a-i) * pows_a(a-i) + enddo + do i = 101, b + P_new2(i) = binom_func(b,b-i) * pows_b(b-i) + enddo + +end + +! --- + diff --git a/src/utils/torus_one_e_integration.irp.f b/src/utils/torus_one_e_integration.irp.f new file mode 100644 index 00000000..f147b0f8 --- /dev/null +++ b/src/utils/torus_one_e_integration.irp.f @@ -0,0 +1,78 @@ + +! --- + +subroutine overlap_gaussian_torus_xyz(A_center, B_center, alpha, beta, power_A, power_B, torus_L, overlap_x, overlap_y, overlap_z, overlap, dim) + + BEGIN_DOC + ! + ! TODO + ! + ! 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' + + implicit none + + integer, intent(in) :: dim ! dimension maximum for the arrays representing the polynomials + integer, intent(in) :: power_A(3), power_B(3) ! power of the x1 functions + double precision, intent(in) :: A_center(3), B_center(3) ! center of the x1 functions + double precision, intent(in) :: alpha, beta + double precision, intent(in) :: torus_L(3) + double precision, intent(out) :: overlap_x, overlap_y, overlap_z, overlap + + integer :: i + integer :: iorder_p(3) + integer :: nmax + double precision :: P_new(0:max_dim,3), P_center(3), fact_p, p + double precision :: fact_p_x, fact_p_y, fact_p_z + double precision :: F_integral_tab(0:max_dim) + double precision, external :: F_integral + + + ! fact_p = fact_p_x x fact_p_y x fact_p_z + call give_explicit_poly_and_gaussian_torus(P_new, P_center, p, fact_p, iorder_p, alpha, beta, power_A, power_B, A_center, B_center, torus_L, 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 + + 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) + + do i = 1, iorder_p(1) + overlap_x = overlap_x + P_new(i,1) * F_integral_tab(i) + enddo + call gaussian_product_torus_x(alpha, A_center(1), beta, B_center(1), torus_L(1), fact_p_x, p, P_center(1)) + overlap_x *= fact_p_x + + do i = 1,iorder_p(2) + overlap_y = overlap_y + P_new(i,2) * F_integral_tab(i) + enddo + call gaussian_product_torus_x(alpha, A_center(2), beta, B_center(2), torus_L(2), fact_p_y, p, P_center(2)) + overlap_y *= fact_p_y + + do i = 1,iorder_p(3) + overlap_z = overlap_z + P_new(i,3) * F_integral_tab(i) + enddo + call gaussian_product_torus_x(alpha, A_center(3), beta, B_center(3), torus_L(3), fact_p_z, p, P_center(3)) + overlap_z *= fact_p_z + + overlap = overlap_x * overlap_y * overlap_z + +end + +! --- +