mirror of
https://github.com/LCPQ/quantum_package
synced 2024-08-25 05:51:46 +02:00
35% acceleration of AO integrals
This commit is contained in:
parent
f99c4c16cf
commit
56df7257af
@ -13,7 +13,12 @@ double precision function ao_bielec_integral(i,j,k,l)
|
|||||||
double precision :: P_new(0:max_dim,3),P_center(3),fact_p,pp
|
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
|
double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq
|
||||||
integer :: iorder_p(3), iorder_q(3)
|
integer :: iorder_p(3), iorder_q(3)
|
||||||
logical :: do_p, do_q, do_r
|
double precision :: ao_bielec_integral_schwartz_accel
|
||||||
|
|
||||||
|
if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then
|
||||||
|
ao_bielec_integral = ao_bielec_integral_schwartz_accel(i,j,k,l)
|
||||||
|
return
|
||||||
|
endif
|
||||||
|
|
||||||
PROVIDE all_utils
|
PROVIDE all_utils
|
||||||
|
|
||||||
@ -24,8 +29,6 @@ double precision function ao_bielec_integral(i,j,k,l)
|
|||||||
num_k = ao_nucl(k)
|
num_k = ao_nucl(k)
|
||||||
num_l = ao_nucl(l)
|
num_l = ao_nucl(l)
|
||||||
ao_bielec_integral = 0.d0
|
ao_bielec_integral = 0.d0
|
||||||
double precision :: thresh
|
|
||||||
! thresh = ao_integrals_threshold
|
|
||||||
|
|
||||||
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
|
do p = 1, 3
|
||||||
@ -39,31 +42,19 @@ double precision function ao_bielec_integral(i,j,k,l)
|
|||||||
L_center(p) = nucl_coord(num_l,p)
|
L_center(p) = nucl_coord(num_l,p)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
do_p = .True.
|
|
||||||
do p = 1, ao_prim_num(i)
|
do p = 1, ao_prim_num(i)
|
||||||
double precision :: coef1
|
double precision :: coef1
|
||||||
if (.not.do_p) exit
|
|
||||||
do_p = .False.
|
|
||||||
do_q = .True.
|
|
||||||
coef1 = ao_coef_transp(p,i)
|
coef1 = ao_coef_transp(p,i)
|
||||||
do q = 1, ao_prim_num(j)
|
do q = 1, ao_prim_num(j)
|
||||||
double precision :: coef2
|
double precision :: coef2
|
||||||
if (.not.do_q) exit
|
|
||||||
do_q = .False.
|
|
||||||
do_r = .True.
|
|
||||||
coef2 = coef1*ao_coef_transp(q,j)
|
coef2 = coef1*ao_coef_transp(q,j)
|
||||||
double precision :: p_inv,q_inv
|
double precision :: p_inv,q_inv
|
||||||
call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
|
call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
|
||||||
ao_expo_transp(p,i),ao_expo_transp(q,j), &
|
ao_expo_transp(p,i),ao_expo_transp(q,j), &
|
||||||
I_power,J_power,I_center,J_center,dim1)
|
I_power,J_power,I_center,J_center,dim1)
|
||||||
if (abs(fact_p) < 1.d-8) then
|
|
||||||
exit
|
|
||||||
endif
|
|
||||||
p_inv = 1.d0/pp
|
p_inv = 1.d0/pp
|
||||||
do r = 1, ao_prim_num(k)
|
do r = 1, ao_prim_num(k)
|
||||||
double precision :: coef3
|
double precision :: coef3
|
||||||
if (.not.do_r) exit
|
|
||||||
do_r = .False.
|
|
||||||
coef3 = coef2*ao_coef_transp(r,k)
|
coef3 = coef2*ao_coef_transp(r,k)
|
||||||
do s = 1, ao_prim_num(l)
|
do s = 1, ao_prim_num(l)
|
||||||
double precision :: coef4
|
double precision :: coef4
|
||||||
@ -72,20 +63,11 @@ double precision function ao_bielec_integral(i,j,k,l)
|
|||||||
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
|
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
|
||||||
ao_expo_transp(r,k),ao_expo_transp(s,l), &
|
ao_expo_transp(r,k),ao_expo_transp(s,l), &
|
||||||
K_power,L_power,K_center,L_center,dim1)
|
K_power,L_power,K_center,L_center,dim1)
|
||||||
if (abs(fact_p*fact_q) < 1.d-8) then
|
|
||||||
exit
|
|
||||||
endif
|
|
||||||
q_inv = 1.d0/qq
|
q_inv = 1.d0/qq
|
||||||
integral = general_primitive_integral(dim1, &
|
integral = general_primitive_integral(dim1, &
|
||||||
P_new,P_center,fact_p,pp,p_inv,iorder_p, &
|
P_new,P_center,fact_p,pp,p_inv,iorder_p, &
|
||||||
Q_new,Q_center,fact_q,qq,q_inv,iorder_q)
|
Q_new,Q_center,fact_q,qq,q_inv,iorder_q)
|
||||||
ao_bielec_integral += coef4 * integral
|
ao_bielec_integral += coef4 * integral
|
||||||
! if (abs(integral) < thresh) then
|
|
||||||
! exit
|
|
||||||
! endif
|
|
||||||
do_r = .True.
|
|
||||||
do_q = .True.
|
|
||||||
do_p = .True.
|
|
||||||
enddo ! s
|
enddo ! s
|
||||||
enddo ! r
|
enddo ! r
|
||||||
enddo ! q
|
enddo ! q
|
||||||
@ -100,20 +82,12 @@ double precision function ao_bielec_integral(i,j,k,l)
|
|||||||
L_power(p) = ao_power(l,p)
|
L_power(p) = ao_power(l,p)
|
||||||
enddo
|
enddo
|
||||||
double precision :: ERI
|
double precision :: ERI
|
||||||
do_p = .True.
|
|
||||||
do p = 1, ao_prim_num(i)
|
do p = 1, ao_prim_num(i)
|
||||||
if (.not.do_p) exit
|
|
||||||
do_p = .False.
|
|
||||||
do_q = .True.
|
|
||||||
coef1 = ao_coef_transp(p,i)
|
coef1 = ao_coef_transp(p,i)
|
||||||
do q = 1, ao_prim_num(j)
|
do q = 1, ao_prim_num(j)
|
||||||
if (.not.do_q) exit
|
|
||||||
do_q = .False.
|
|
||||||
do_r = .True.
|
|
||||||
coef2 = coef1*ao_coef_transp(q,j)
|
coef2 = coef1*ao_coef_transp(q,j)
|
||||||
do r = 1, ao_prim_num(k)
|
do r = 1, ao_prim_num(k)
|
||||||
if (.not.do_r) exit
|
|
||||||
do_r = .False.
|
|
||||||
coef3 = coef2*ao_coef_transp(r,k)
|
coef3 = coef2*ao_coef_transp(r,k)
|
||||||
do s = 1, ao_prim_num(l)
|
do s = 1, ao_prim_num(l)
|
||||||
coef4 = coef3*ao_coef_transp(s,l)
|
coef4 = coef3*ao_coef_transp(s,l)
|
||||||
@ -123,12 +97,6 @@ double precision function ao_bielec_integral(i,j,k,l)
|
|||||||
I_power(2),J_power(2),K_power(2),L_power(2), &
|
I_power(2),J_power(2),K_power(2),L_power(2), &
|
||||||
I_power(3),J_power(3),K_power(3),L_power(3))
|
I_power(3),J_power(3),K_power(3),L_power(3))
|
||||||
ao_bielec_integral += coef4 * integral
|
ao_bielec_integral += coef4 * integral
|
||||||
! if (abs(integral) < thresh) then
|
|
||||||
! exit
|
|
||||||
! endif
|
|
||||||
do_r = .True.
|
|
||||||
do_q = .True.
|
|
||||||
do_p = .True.
|
|
||||||
enddo ! s
|
enddo ! s
|
||||||
enddo ! r
|
enddo ! r
|
||||||
enddo ! q
|
enddo ! q
|
||||||
@ -138,6 +106,154 @@ double precision function ao_bielec_integral(i,j,k,l)
|
|||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
|
double precision function ao_bielec_integral_schwartz_accel(i,j,k,l)
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! integral of the AO basis <ik|jl> or (ij|kl)
|
||||||
|
! i(r1) j(r1) 1/r12 k(r2) l(r2)
|
||||||
|
END_DOC
|
||||||
|
integer,intent(in) :: i,j,k,l
|
||||||
|
integer :: p,q,r,s
|
||||||
|
double precision :: I_center(3),J_center(3),K_center(3),L_center(3)
|
||||||
|
integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3)
|
||||||
|
double precision :: integral
|
||||||
|
include 'include/constants.F'
|
||||||
|
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
|
||||||
|
integer :: iorder_p(3), iorder_q(3)
|
||||||
|
double precision, allocatable :: schwartz_kl(:,:)
|
||||||
|
double precision :: schwartz_ij
|
||||||
|
|
||||||
|
PROVIDE all_utils
|
||||||
|
|
||||||
|
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_bielec_integral_schwartz_accel = 0.d0
|
||||||
|
double precision :: thresh
|
||||||
|
thresh = ao_integrals_threshold*ao_integrals_threshold
|
||||||
|
|
||||||
|
allocate(schwartz_kl(ao_prim_num(k),ao_prim_num(l)))
|
||||||
|
|
||||||
|
|
||||||
|
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
|
||||||
|
|
||||||
|
do r = 1, ao_prim_num(k)
|
||||||
|
do s = 1, ao_prim_num(l)
|
||||||
|
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
|
||||||
|
ao_expo_transp(r,k),ao_expo_transp(s,l), &
|
||||||
|
K_power,L_power,K_center,L_center,dim1)
|
||||||
|
q_inv = 1.d0/qq
|
||||||
|
schwartz_kl(r,s) = general_primitive_integral(dim1, &
|
||||||
|
Q_new,Q_center,fact_q,qq,q_inv,iorder_q, &
|
||||||
|
Q_new,Q_center,fact_q,qq,q_inv,iorder_q)
|
||||||
|
schwartz_kl(r,s) = schwartz_kl(r,s)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
do p = 1, ao_prim_num(i)
|
||||||
|
double precision :: coef1
|
||||||
|
coef1 = ao_coef_transp(p,i)
|
||||||
|
do q = 1, ao_prim_num(j)
|
||||||
|
double precision :: coef2
|
||||||
|
coef2 = coef1*ao_coef_transp(q,j)
|
||||||
|
double precision :: p_inv,q_inv
|
||||||
|
call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
|
||||||
|
ao_expo_transp(p,i),ao_expo_transp(q,j), &
|
||||||
|
I_power,J_power,I_center,J_center,dim1)
|
||||||
|
p_inv = 1.d0/pp
|
||||||
|
schwartz_ij = general_primitive_integral(dim1, &
|
||||||
|
P_new,P_center,fact_p,pp,p_inv,iorder_p, &
|
||||||
|
P_new,P_center,fact_p,pp,p_inv,iorder_p)
|
||||||
|
schwartz_ij = schwartz_ij
|
||||||
|
do r = 1, ao_prim_num(k)
|
||||||
|
double precision :: coef3
|
||||||
|
coef3 = coef2*ao_coef_transp(r,k)
|
||||||
|
do s = 1, ao_prim_num(l)
|
||||||
|
double precision :: coef4
|
||||||
|
coef4 = coef3*ao_coef_transp(s,l)
|
||||||
|
if (schwartz_kl(r,s)*schwartz_ij*coef4*coef4 < thresh) then
|
||||||
|
cycle
|
||||||
|
endif
|
||||||
|
double precision :: general_primitive_integral
|
||||||
|
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
|
||||||
|
ao_expo_transp(r,k),ao_expo_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_bielec_integral_schwartz_accel += 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 r = 1, ao_prim_num(k)
|
||||||
|
do s = 1, ao_prim_num(l)
|
||||||
|
schwartz_kl(r,s) = ERI( &
|
||||||
|
ao_expo_transp(r,k),ao_expo_transp(s,l),ao_expo_transp(r,k),ao_expo_transp(s,l),&
|
||||||
|
K_power(1),L_power(1),K_power(1),L_power(1), &
|
||||||
|
K_power(2),L_power(2),K_power(2),L_power(2), &
|
||||||
|
K_power(3),L_power(3),K_power(3),L_power(3))
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
do p = 1, ao_prim_num(i)
|
||||||
|
coef1 = ao_coef_transp(p,i)
|
||||||
|
do q = 1, ao_prim_num(j)
|
||||||
|
coef2 = coef1*ao_coef_transp(q,j)
|
||||||
|
schwartz_ij = ERI( &
|
||||||
|
ao_expo_transp(p,i),ao_expo_transp(q,j),ao_expo_transp(p,i),ao_expo_transp(q,j),&
|
||||||
|
I_power(1),J_power(1),I_power(1),J_power(1), &
|
||||||
|
I_power(2),J_power(2),I_power(2),J_power(2), &
|
||||||
|
I_power(3),J_power(3),I_power(3),J_power(3))
|
||||||
|
do r = 1, ao_prim_num(k)
|
||||||
|
coef3 = coef2*ao_coef_transp(r,k)
|
||||||
|
do s = 1, ao_prim_num(l)
|
||||||
|
coef4 = coef3*ao_coef_transp(s,l)
|
||||||
|
if (schwartz_kl(r,s)*schwartz_ij*coef4*coef4 < thresh) then
|
||||||
|
cycle
|
||||||
|
endif
|
||||||
|
integral = ERI( &
|
||||||
|
ao_expo_transp(p,i),ao_expo_transp(q,j),ao_expo_transp(r,k),ao_expo_transp(s,l),&
|
||||||
|
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_bielec_integral_schwartz_accel += coef4 * integral
|
||||||
|
enddo ! s
|
||||||
|
enddo ! r
|
||||||
|
enddo ! q
|
||||||
|
enddo ! p
|
||||||
|
|
||||||
|
endif
|
||||||
|
deallocate (schwartz_kl)
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
integer function ao_l4(i,j,k,l)
|
integer function ao_l4(i,j,k,l)
|
||||||
implicit none
|
implicit none
|
||||||
@ -501,6 +617,9 @@ double precision function general_primitive_integral(dim, &
|
|||||||
|
|
||||||
!DEC$ FORCEINLINE
|
!DEC$ FORCEINLINE
|
||||||
call multiply_poly(Ix_pol,n_Ix,Iy_pol,n_Iy,d_poly,n_pt_tmp)
|
call multiply_poly(Ix_pol,n_Ix,Iy_pol,n_Iy,d_poly,n_pt_tmp)
|
||||||
|
if (n_pt_tmp == -1) then
|
||||||
|
return
|
||||||
|
endif
|
||||||
n_pt_out = n_pt_tmp+n_Iz
|
n_pt_out = n_pt_tmp+n_Iz
|
||||||
do i=0,n_pt_out
|
do i=0,n_pt_out
|
||||||
d1(i)=0.d0
|
d1(i)=0.d0
|
||||||
|
Loading…
Reference in New Issue
Block a user