LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-23 12:56:14 +01:00
Code Issues Releases Wiki Activity

Gained another 21% on AO integrals

Browse Source
This commit is contained in:
Anthony Scemama 2014-09-24 23:23:17 +02:00
parent 56df7257af
commit 1efb1c3687
1 changed files with 44 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

64
src/BiInts/ao_bi_integrals.irp.f
View File

@ -136,7 +136,7 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l)
double precision :: thresh double precision :: thresh
thresh = ao_integrals_threshold*ao_integrals_threshold thresh = ao_integrals_threshold*ao_integrals_threshold
allocate(schwartz_kl(ao_prim_num(k),ao_prim_num(l))) allocate(schwartz_kl(0:ao_prim_num(l),0:ao_prim_num(k)))
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
@ -151,17 +151,23 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l)
L_center(p) = nucl_coord(num_l,p) L_center(p) = nucl_coord(num_l,p)
enddo enddo
schwartz_kl(0,0) = 0.d0
do r = 1, ao_prim_num(k) do r = 1, ao_prim_num(k)
coef1 = ao_coef_transp(r,k)*ao_coef_transp(r,k)
schwartz_kl(0,r) = 0.d0
do s = 1, ao_prim_num(l) do s = 1, ao_prim_num(l)
coef2 = coef1 * ao_coef_transp(s,l) * ao_coef_transp(s,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)
q_inv = 1.d0/qq q_inv = 1.d0/qq
schwartz_kl(r,s) = general_primitive_integral(dim1, & schwartz_kl(s,r) = 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, &
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) * coef2
schwartz_kl(0,r) = max(schwartz_kl(0,r),schwartz_kl(s,r))
enddo enddo
schwartz_kl(0,0) = max(schwartz_kl(0,r),schwartz_kl(0,0))
enddo enddo
do p = 1, ao_prim_num(i) do p = 1, ao_prim_num(i)
@ -177,17 +183,23 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l)
p_inv = 1.d0/pp p_inv = 1.d0/pp
schwartz_ij = general_primitive_integral(dim1, & 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, &
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 coef2*coef2
if (schwartz_kl(0,0)*schwartz_ij < thresh) then
cycle
endif
do r = 1, ao_prim_num(k) do r = 1, ao_prim_num(k)
if (schwartz_kl(0,r)*schwartz_ij < thresh) then
cycle
endif
double precision :: coef3 double precision :: coef3
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
coef4 = coef3*ao_coef_transp(s,l) if (schwartz_kl(s,r)*schwartz_ij < thresh) then
if (schwartz_kl(r,s)*schwartz_ij*coef4*coef4 < thresh) then
cycle cycle
endif endif
coef4 = coef3*ao_coef_transp(s,l)
double precision :: general_primitive_integral double precision :: general_primitive_integral
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), &
@ -212,14 +224,21 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l)
enddo enddo
double precision :: ERI double precision :: ERI
schwartz_kl(0,0) = 0.d0
do r = 1, ao_prim_num(k) do r = 1, ao_prim_num(k)
coef1 = ao_coef_transp(r,k)*ao_coef_transp(r,k)
schwartz_kl(0,r) = 0.d0
do s = 1, ao_prim_num(l) do s = 1, ao_prim_num(l)
schwartz_kl(r,s) = ERI( & coef2 = coef1*ao_coef_transp(s,l)*ao_coef_transp(s,l)
schwartz_kl(s,r) = ERI( &
ao_expo_transp(r,k),ao_expo_transp(s,l),ao_expo_transp(r,k),ao_expo_transp(s,l),& 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(1),L_power(1),K_power(1),L_power(1), &
K_power(2),L_power(2),K_power(2),L_power(2), & K_power(2),L_power(2),K_power(2),L_power(2), &
K_power(3),L_power(3),K_power(3),L_power(3)) K_power(3),L_power(3),K_power(3),L_power(3)) * &
coef2
schwartz_kl(0,r) = max(schwartz_kl(0,r),schwartz_kl(s,r))
enddo enddo
schwartz_kl(0,0) = max(schwartz_kl(0,r),schwartz_kl(0,0))
enddo enddo
do p = 1, ao_prim_num(i) do p = 1, ao_prim_num(i)
@ -230,14 +249,20 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l)
ao_expo_transp(p,i),ao_expo_transp(q,j),ao_expo_transp(p,i),ao_expo_transp(q,j),& 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(1),J_power(1),I_power(1),J_power(1), &
I_power(2),J_power(2),I_power(2),J_power(2), & I_power(2),J_power(2),I_power(2),J_power(2), &
I_power(3),J_power(3),I_power(3),J_power(3)) I_power(3),J_power(3),I_power(3),J_power(3))*coef2*coef2
do r = 1, ao_prim_num(k) if (schwartz_kl(0,0)*schwartz_ij < thresh) then
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 cycle
endif endif
do r = 1, ao_prim_num(k)
if (schwartz_kl(0,r)*schwartz_ij < thresh) then
cycle
endif
coef3 = coef2*ao_coef_transp(r,k)
do s = 1, ao_prim_num(l)
if (schwartz_kl(s,r)*schwartz_ij < thresh) then
cycle
endif
coef4 = coef3*ao_coef_transp(s,l)
integral = ERI( & integral = ERI( &
ao_expo_transp(p,i),ao_expo_transp(q,j),ao_expo_transp(r,k),ao_expo_transp(s,l),& 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(1),J_power(1),K_power(1),L_power(1), &
@ -663,8 +688,9 @@ double precision function ERI(alpha,beta,delta,gama,a_x,b_x,c_x,d_x,a_y,b_y,c_y,
p = alpha + beta p = alpha + beta
q = delta + gama q = delta + gama
ASSERT (p+q >= 0.d0) ASSERT (p+q >= 0.d0)
n_pt = n_pt_sup(a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z)
coeff = pi_5_2 / (p * q * dsqrt(p+q)) coeff = pi_5_2 / (p * q * dsqrt(p+q))
!DIR$ FORCEINLINE
n_pt = n_pt_sup(a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z)
if (n_pt == 0) then if (n_pt == 0) then
ERI = coeff ERI = coeff
@ -768,12 +794,12 @@ end
integer function n_pt_sup(a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z) integer function n_pt_sup(a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Returns the upper boundary of the degree of the polynom involved in the ! Returns the upper boundary of the degree of the polynomial involved in the
! bielctronic integral : ! bielctronic integral :
! Ix(a_x,b_x,c_x,d_x) * Iy(a_y,b_y,c_y,d_y) * Iz(a_z,b_z,c_z,d_z) ! Ix(a_x,b_x,c_x,d_x) * Iy(a_y,b_y,c_y,d_y) * Iz(a_z,b_z,c_z,d_z)
END_DOC END_DOC
integer :: a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z integer :: a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z
n_pt_sup = 2 * ( (a_x+b_x+c_x+d_x) + (a_y+b_y+c_y+d_y) + (a_z+b_z+c_z+d_z) ) n_pt_sup = ishft( a_x+b_x+c_x+d_x + a_y+b_y+c_y+d_y + a_z+b_z+c_z+d_z,1 )
end end