9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2025-01-03 09:05:39 +01:00

Introduced ao_two_e_integral_general

This commit is contained in:
Anthony Scemama 2024-12-02 16:47:57 +01:00
parent 7a0194d576
commit 80e0aaec6b

View File

@ -40,8 +40,11 @@ double precision function ao_two_e_integral(i, j, k, l)
double precision, external :: ao_two_e_integral_erf double precision, external :: ao_two_e_integral_erf
double precision, external :: ao_two_e_integral_cgtos double precision, external :: ao_two_e_integral_cgtos
double precision, external :: ao_two_e_integral_schwartz_accel double precision, external :: ao_two_e_integral_schwartz_accel
double precision, external :: ao_two_e_integral_general
double precision, external :: general_primitive_integral
logical, external :: do_schwartz_accel logical, external :: do_schwartz_accel
double precision :: coef1, coef2, coef3, coef4
if(use_cgtos) then if(use_cgtos) then
@ -58,83 +61,44 @@ double precision function ao_two_e_integral(i, j, k, l)
else else
dim1 = n_pt_max_integrals
num_i = ao_nucl(i) num_i = ao_nucl(i)
num_j = ao_nucl(j) num_j = ao_nucl(j)
num_k = ao_nucl(k) num_k = ao_nucl(k)
num_l = ao_nucl(l) num_l = ao_nucl(l)
ao_two_e_integral = 0.d0 ao_two_e_integral = 0.d0
if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k) then
do p = 1, 3
I_power(p) = ao_power(i,p)
J_power(p) = ao_power(j,p)
K_power(p) = ao_power(k,p)
L_power(p) = ao_power(l,p)
I_center(p) = nucl_coord(num_i,p)
J_center(p) = nucl_coord(num_j,p)
K_center(p) = nucl_coord(num_k,p)
L_center(p) = nucl_coord(num_l,p)
enddo
double precision :: coef1, coef2, coef3, coef4 ao_two_e_integral = ao_two_e_integral_general(i,j,k,l,general_primitive_integral)
double precision :: p_inv,q_inv
double precision :: general_primitive_integral
do p = 1, ao_prim_num(i) else
coef1 = ao_coef_normalized_ordered_transp(p,i)
do q = 1, ao_prim_num(j)
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), &
I_power,J_power,I_center,J_center,dim1)
p_inv = 1.d0/pp
do r = 1, ao_prim_num(k)
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
do s = 1, ao_prim_num(l)
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), &
K_power,L_power,K_center,L_center,dim1)
q_inv = 1.d0/qq
integral = general_primitive_integral(dim1, &
P_new,P_center,fact_p,pp,p_inv,iorder_p, &
Q_new,Q_center,fact_q,qq,q_inv,iorder_q)
ao_two_e_integral = ao_two_e_integral + coef4 * integral
enddo ! s
enddo ! r
enddo ! q
enddo ! p
else do p = 1, 3
I_power(p) = ao_power(i,p)
J_power(p) = ao_power(j,p)
K_power(p) = ao_power(k,p)
L_power(p) = ao_power(l,p)
enddo
double precision :: ERI
do p = 1, 3 do p = 1, ao_prim_num(i)
I_power(p) = ao_power(i,p) coef1 = ao_coef_normalized_ordered_transp(p,i)
J_power(p) = ao_power(j,p) do q = 1, ao_prim_num(j)
K_power(p) = ao_power(k,p) coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
L_power(p) = ao_power(l,p) do r = 1, ao_prim_num(k)
enddo coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
double precision :: ERI do s = 1, ao_prim_num(l)
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
do p = 1, ao_prim_num(i) integral = ERI( &
coef1 = ao_coef_normalized_ordered_transp(p,i) ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),&
do q = 1, ao_prim_num(j) I_power(1),J_power(1),K_power(1),L_power(1), &
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) I_power(2),J_power(2),K_power(2),L_power(2), &
do r = 1, ao_prim_num(k) I_power(3),J_power(3),K_power(3),L_power(3))
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) ao_two_e_integral = ao_two_e_integral + coef4 * integral
do s = 1, ao_prim_num(l) enddo ! s
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) enddo ! r
integral = ERI( & enddo ! q
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),& enddo ! p
I_power(1),J_power(1),K_power(1),L_power(1), &
I_power(2),J_power(2),K_power(2),L_power(2), &
I_power(3),J_power(3),K_power(3),L_power(3))
ao_two_e_integral = ao_two_e_integral + coef4 * integral
enddo ! s
enddo ! r
enddo ! q
enddo ! p
endif endif
@ -144,6 +108,76 @@ end
! --- ! ---
double precision function ao_two_e_integral_general(i, j, k, l, op)
BEGIN_DOC
! integral of the AO basis <ik|jl> or (ij|kl)
! i(r1) j(r1) 1/r12 k(r2) l(r2)
END_DOC
implicit none
include 'utils/constants.include.F'
integer, intent(in) :: i, j, k, l
double precision, external :: op ! Operator function
integer :: p, q, r, s
integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3)
integer :: iorder_p(3), iorder_q(3)
double precision :: I_center(3), J_center(3), K_center(3), L_center(3)
double precision :: integral
double precision :: P_new(0:max_dim,3),P_center(3),fact_p,pp
double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq
dim1 = n_pt_max_integrals
num_i = ao_nucl(i)
num_j = ao_nucl(j)
num_k = ao_nucl(k)
num_l = ao_nucl(l)
ao_two_e_integral_general = 0.d0
do p = 1, 3
I_power(p) = ao_power(i,p)
J_power(p) = ao_power(j,p)
K_power(p) = ao_power(k,p)
L_power(p) = ao_power(l,p)
I_center(p) = nucl_coord(num_i,p)
J_center(p) = nucl_coord(num_j,p)
K_center(p) = nucl_coord(num_k,p)
L_center(p) = nucl_coord(num_l,p)
enddo
double precision :: coef1, coef2, coef3, coef4
double precision :: p_inv,q_inv
do p = 1, ao_prim_num(i)
coef1 = ao_coef_normalized_ordered_transp(p,i)
do q = 1, ao_prim_num(j)
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), &
I_power,J_power,I_center,J_center,dim1)
p_inv = 1.d0/pp
do r = 1, ao_prim_num(k)
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
do s = 1, ao_prim_num(l)
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), &
K_power,L_power,K_center,L_center,dim1)
q_inv = 1.d0/qq
integral = op(dim1, &
P_new,P_center,fact_p,pp,p_inv,iorder_p, &
Q_new,Q_center,fact_q,qq,q_inv,iorder_q)
ao_two_e_integral_general = ao_two_e_integral_general + coef4 * integral
enddo ! s
enddo ! r
enddo ! q
enddo ! p
end
double precision function ao_two_e_integral_schwartz_accel(i,j,k,l) double precision function ao_two_e_integral_schwartz_accel(i,j,k,l)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
@ -512,7 +546,7 @@ double precision function general_primitive_integral(dim, &
double precision :: a,b,c,d,e,f,accu,pq,const double precision :: a,b,c,d,e,f,accu,pq,const
double precision :: pq_inv, p10_1, p10_2, p01_1, p01_2,pq_inv_2 double precision :: pq_inv, p10_1, p10_2, p01_1, p01_2,pq_inv_2
integer :: n_pt_tmp,n_pt_out, iorder integer :: n_pt_tmp,n_pt_out, iorder
double precision :: d1(0:max_dim),d_poly(0:max_dim),rint,d1_screened(0:max_dim) double precision :: d1(0:max_dim),d_poly(0:max_dim),d1_screened(0:max_dim)
general_primitive_integral = 0.d0 general_primitive_integral = 0.d0