9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-10-06 16:15:57 +02:00
qp2/plugins/local/ao_many_one_e_ints/ao_gaus_gauss.irp.f
2024-01-15 12:02:38 +01:00

506 lines
15 KiB
Fortran

! ---
subroutine overlap_gauss_xyz_r12_ao(D_center,delta,i,j,gauss_ints)
implicit none
BEGIN_DOC
! gauss_ints(m) = \int dr AO_i(r) AO_j(r) x/y/z e^{-delta |r-D_center|^2}
!
! with m == 1 ==> x, m == 2 ==> y, m == 3 ==> z
END_DOC
integer, intent(in) :: i,j
double precision, intent(in) :: D_center(3), delta
double precision, intent(out) :: gauss_ints(3)
integer :: num_a,num_b,power_A(3), power_B(3),l,k,m
double precision :: A_center(3), B_center(3),overlap_gauss_r12,alpha,beta,gauss_ints_tmp(3)
gauss_ints = 0.d0
if(ao_overlap_abs(j,i).lt.1.d-12)then
return
endif
num_A = ao_nucl(i)
power_A(1:3)= ao_power(i,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
num_B = ao_nucl(j)
power_B(1:3)= ao_power(j,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do l=1,ao_prim_num(i)
alpha = ao_expo_ordered_transp(l,i)
do k=1,ao_prim_num(j)
beta = ao_expo_ordered_transp(k,j)
call overlap_gauss_xyz_r12(D_center,delta,A_center,B_center,power_A,power_B,alpha,beta,gauss_ints_tmp)
do m = 1, 3
gauss_ints(m) += gauss_ints_tmp(m) * ao_coef_normalized_ordered_transp(l,i) &
* ao_coef_normalized_ordered_transp(k,j)
enddo
enddo
enddo
end
double precision function overlap_gauss_xyz_r12_ao_specific(D_center,delta,i,j,mx)
implicit none
BEGIN_DOC
! \int dr AO_i(r) AO_j(r) x/y/z e^{-delta |r-D_center|^2}
!
! with mx == 1 ==> x, mx == 2 ==> y, mx == 3 ==> z
END_DOC
integer, intent(in) :: i,j,mx
double precision, intent(in) :: D_center(3), delta
integer :: num_a,num_b,power_A(3), power_B(3),l,k
double precision :: gauss_int
double precision :: A_center(3), B_center(3),overlap_gauss_r12,alpha,beta
double precision :: overlap_gauss_xyz_r12_specific
overlap_gauss_xyz_r12_ao_specific = 0.d0
if(ao_overlap_abs(j,i).lt.1.d-12)then
return
endif
num_A = ao_nucl(i)
power_A(1:3)= ao_power(i,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
num_B = ao_nucl(j)
power_B(1:3)= ao_power(j,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do l=1,ao_prim_num(i)
alpha = ao_expo_ordered_transp(l,i)
do k=1,ao_prim_num(j)
beta = ao_expo_ordered_transp(k,j)
gauss_int = overlap_gauss_xyz_r12_specific(D_center,delta,A_center,B_center,power_A,power_B,alpha,beta,mx)
overlap_gauss_xyz_r12_ao_specific = gauss_int * ao_coef_normalized_ordered_transp(l,i) &
* ao_coef_normalized_ordered_transp(k,j)
enddo
enddo
end
subroutine overlap_gauss_r12_all_ao(D_center,delta,aos_ints)
implicit none
double precision, intent(in) :: D_center(3), delta
double precision, intent(out):: aos_ints(ao_num,ao_num)
integer :: num_a,num_b,power_A(3), power_B(3),l,k,i,j
double precision :: A_center(3), B_center(3),overlap_gauss_r12,alpha,beta,analytical_j
aos_ints = 0.d0
do i = 1, ao_num
do j = 1, ao_num
if(ao_overlap_abs(j,i).lt.1.d-12)cycle
num_A = ao_nucl(i)
power_A(1:3)= ao_power(i,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
num_B = ao_nucl(j)
power_B(1:3)= ao_power(j,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do l=1,ao_prim_num(i)
alpha = ao_expo_ordered_transp(l,i)
do k=1,ao_prim_num(j)
beta = ao_expo_ordered_transp(k,j)
analytical_j = overlap_gauss_r12(D_center,delta,A_center,B_center,power_A,power_B,alpha,beta)
aos_ints(j,i) += analytical_j * ao_coef_normalized_ordered_transp(l,i) &
* ao_coef_normalized_ordered_transp(k,j)
enddo
enddo
enddo
enddo
end
! ---
! TODO :: PUT CYCLES IN LOOPS
double precision function overlap_gauss_r12_ao(D_center, delta, i, j)
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
double precision, intent(in) :: D_center(3), delta
integer :: power_A(3), power_B(3), l, k
double precision :: A_center(3), B_center(3), alpha, beta, coef, coef1, analytical_j
double precision, external :: overlap_gauss_r12
overlap_gauss_r12_ao = 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)
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
analytical_j = overlap_gauss_r12(D_center, delta, A_center, B_center, power_A, power_B, alpha, beta)
overlap_gauss_r12_ao += coef * analytical_j
enddo
enddo
end
! --
double precision function overlap_abs_gauss_r12_ao(D_center, delta, i, j)
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
double precision, intent(in) :: D_center(3), delta
integer :: power_A(3), power_B(3), l, k
double precision :: A_center(3), B_center(3), alpha, beta, coef, coef1, analytical_j
double precision, external :: overlap_abs_gauss_r12
overlap_abs_gauss_r12_ao = 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)
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
analytical_j = overlap_abs_gauss_r12(D_center, delta, A_center, B_center, power_A, power_B, alpha, beta)
overlap_abs_gauss_r12_ao += dabs(coef * analytical_j)
enddo
enddo
end
! --
subroutine overlap_gauss_r12_ao_v(D_center, LD_D, delta, i, j, resv, LD_resv, n_points)
BEGIN_DOC
!
! \int dr AO_i(r) AO_j(r) e^{-delta |r-D_center|^2}
!
! n_points: nb of integrals <= min(LD_D, LD_resv)
!
END_DOC
implicit none
integer, intent(in) :: i, j, LD_D, LD_resv, n_points
double precision, intent(in) :: D_center(LD_D,3), delta
double precision, intent(out) :: resv(LD_resv)
integer :: ipoint
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(:)
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, LD_D, delta, A_center, B_center, power_A, power_B, alpha, beta, analytical_j, n_points, 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)
BEGIN_DOC
!
! \int dr AO_i(r) AO_j(r) e^{-beta |r-B_center^2|} e^{-delta |r-D_center|^2}
!
END_DOC
implicit none
integer, intent(in) :: i, j
double precision, intent(in) :: B_center(3), beta, D_center(3), delta
integer :: power_A1(3), power_A2(3), l, k
double precision :: A1_center(3), A2_center(3), alpha1, alpha2, coef1, coef12, analytical_j
double precision :: G_center(3), gama, fact_g, gama_inv
double precision, external :: overlap_gauss_r12, overlap_gauss_r12_ao
if(beta .lt. 1d-10) then
overlap_gauss_r12_ao_with1s = overlap_gauss_r12_ao(D_center, delta, i, j)
return
endif
overlap_gauss_r12_ao_with1s = 0.d0
if(ao_overlap_abs(j,i) .lt. 1.d-12) then
return
endif
! e^{-beta |r-B_center^2|} e^{-delta |r-D_center|^2} = fact_g e^{-gama |r - G|^2}
gama = beta + delta
gama_inv = 1.d0 / gama
G_center(1) = (beta * B_center(1) + delta * D_center(1)) * gama_inv
G_center(2) = (beta * B_center(2) + delta * D_center(2)) * gama_inv
G_center(3) = (beta * B_center(3) + delta * D_center(3)) * gama_inv
fact_g = beta * delta * gama_inv * ( (B_center(1) - D_center(1)) * (B_center(1) - D_center(1)) &
+ (B_center(2) - D_center(2)) * (B_center(2) - D_center(2)) &
+ (B_center(3) - D_center(3)) * (B_center(3) - D_center(3)) )
if(fact_g .gt. 10d0) return
fact_g = dexp(-fact_g)
! ---
power_A1(1:3) = ao_power(i,1:3)
power_A2(1:3) = ao_power(j,1:3)
A1_center(1:3) = nucl_coord(ao_nucl(i),1:3)
A2_center(1:3) = nucl_coord(ao_nucl(j),1:3)
do l = 1, ao_prim_num(i)
alpha1 = ao_expo_ordered_transp (l,i)
coef1 = fact_g * ao_coef_normalized_ordered_transp(l,i)
if(dabs(coef1) .lt. 1d-12) cycle
do k = 1, ao_prim_num(j)
alpha2 = ao_expo_ordered_transp (k,j)
coef12 = coef1 * ao_coef_normalized_ordered_transp(k,j)
if(dabs(coef12) .lt. 1d-12) cycle
analytical_j = overlap_gauss_r12(G_center, gama, A1_center, A2_center, power_A1, power_A2, alpha1, alpha2)
overlap_gauss_r12_ao_with1s += coef12 * analytical_j
enddo
enddo
end
! ---
subroutine overlap_gauss_r12_ao_with1s_v(B_center, beta, D_center, LD_D, delta, i, j, resv, LD_resv, n_points)
BEGIN_DOC
!
! \int dr AO_i(r) AO_j(r) e^{-beta |r-B_center^2|} e^{-delta |r-D_center|^2}
! using an array of D_centers.
!
END_DOC
implicit none
integer, intent(in) :: i, j, n_points, LD_D, LD_resv
double precision, intent(in) :: B_center(3), beta, D_center(LD_D,3), delta
double precision, intent(out) :: resv(LD_resv)
integer :: ipoint
integer :: power_A1(3), power_A2(3), l, k
double precision :: A1_center(3), A2_center(3), alpha1, alpha2, coef1
double precision :: coef12, coef12f
double precision :: gama, gama_inv
double precision :: bg, dg, bdg
double precision, allocatable :: fact_g(:), G_center(:,:), analytical_j(:)
if(ao_overlap_abs(j,i) .lt. 1.d-12) then
return
endif
ASSERT(beta .gt. 0.d0)
if(beta .lt. 1d-10) then
call overlap_gauss_r12_ao_v(D_center, LD_D, delta, i, j, resv, LD_resv, n_points)
return
endif
resv(:) = 0.d0
! e^{-beta |r-B_center^2|} e^{-delta |r-D_center|^2} = fact_g e^{-gama |r - G|^2}
gama = beta + delta
gama_inv = 1.d0 / gama
power_A1(1:3) = ao_power(i,1:3)
power_A2(1:3) = ao_power(j,1:3)
A1_center(1:3) = nucl_coord(ao_nucl(i),1:3)
A2_center(1:3) = nucl_coord(ao_nucl(j),1:3)
allocate(fact_g(n_points), G_center(n_points,3), analytical_j(n_points))
bg = beta * gama_inv
dg = delta * gama_inv
bdg = bg * delta
do ipoint = 1, n_points
G_center(ipoint,1) = bg * B_center(1) + dg * D_center(ipoint,1)
G_center(ipoint,2) = bg * B_center(2) + dg * D_center(ipoint,2)
G_center(ipoint,3) = bg * B_center(3) + dg * D_center(ipoint,3)
fact_g(ipoint) = bdg * ( (B_center(1) - D_center(ipoint,1)) * (B_center(1) - D_center(ipoint,1)) &
+ (B_center(2) - D_center(ipoint,2)) * (B_center(2) - D_center(ipoint,2)) &
+ (B_center(3) - D_center(ipoint,3)) * (B_center(3) - D_center(ipoint,3)) )
if(fact_g(ipoint) < 10d0) then
fact_g(ipoint) = dexp(-fact_g(ipoint))
else
fact_g(ipoint) = 0.d0
endif
enddo
do l = 1, ao_prim_num(i)
alpha1 = ao_expo_ordered_transp (l,i)
coef1 = ao_coef_normalized_ordered_transp(l,i)
do k = 1, ao_prim_num(j)
alpha2 = ao_expo_ordered_transp (k,j)
coef12 = coef1 * ao_coef_normalized_ordered_transp(k,j)
if(dabs(coef12) .lt. 1d-12) cycle
call overlap_gauss_r12_v(G_center, n_points, gama, A1_center, A2_center, power_A1, power_A2, alpha1, alpha2, analytical_j, n_points, n_points)
do ipoint = 1, n_points
coef12f = coef12 * fact_g(ipoint)
resv(ipoint) += coef12f * analytical_j(ipoint)
enddo
enddo
enddo
deallocate(fact_g, G_center, analytical_j)
end
! ---
subroutine overlap_gauss_r12_ao_012(D_center, delta, i, j, ints)
BEGIN_DOC
!
! Computes the following integrals :
!
! ints(1) = $\int_{-\infty}^{infty} dr x^0 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
!
! ints(2) = $\int_{-\infty}^{infty} dr x^1 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
! ints(3) = $\int_{-\infty}^{infty} dr y^1 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
! ints(4) = $\int_{-\infty}^{infty} dr z^1 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
!
! ints(5) = $\int_{-\infty}^{infty} dr x^2 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
! ints(6) = $\int_{-\infty}^{infty} dr y^2 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
! ints(7) = $\int_{-\infty}^{infty} dr z^2 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
!
END_DOC
include 'utils/constants.include.F'
implicit none
integer, intent(in) :: i, j
double precision, intent(in) :: delta, D_center(3)
double precision, intent(out) :: ints(7)
integer :: k, l, m
integer :: power_A(3), power_B(3), power_A1(3), power_A2(3)
double precision :: A_center(3), B_center(3), alpha, beta, coef1, coef
double precision :: integral0, integral1, integral2
double precision, external :: overlap_gauss_r12
ints = 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)
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
integral0 = overlap_gauss_r12(D_center, delta, A_center, B_center, power_A, power_B, alpha, beta)
ints(1) += coef * integral0
do m = 1, 3
power_A1 = power_A
power_A1(m) += 1
integral1 = overlap_gauss_r12(D_center, delta, A_center, B_center, power_A1, power_B, alpha, beta)
ints(1+m) += coef * (integral1 + A_center(m)*integral0)
power_A2 = power_A
power_A2(m) += 2
integral2 = overlap_gauss_r12(D_center, delta, A_center, B_center, power_A2, power_B, alpha, beta)
ints(4+m) += coef * (integral2 + A_center(m) * (2.d0*integral1 + A_center(m)*integral0))
enddo
enddo ! k
enddo ! l
return
end
! ---