diff --git a/src/BiInts/ao_bi_integrals.irp.f b/src/BiInts/ao_bi_integrals.irp.f index 990e79a6..ea7c21a4 100644 --- a/src/BiInts/ao_bi_integrals.irp.f +++ b/src/BiInts/ao_bi_integrals.irp.f @@ -136,7 +136,7 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l) double precision :: thresh 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 @@ -151,17 +151,23 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l) L_center(p) = nucl_coord(num_l,p) enddo + schwartz_kl(0,0) = 0.d0 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) + 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,& 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, & + 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) - schwartz_kl(r,s) = schwartz_kl(r,s) + Q_new,Q_center,fact_q,qq,q_inv,iorder_q) & + * coef2 + schwartz_kl(0,r) = max(schwartz_kl(0,r),schwartz_kl(s,r)) enddo + schwartz_kl(0,0) = max(schwartz_kl(0,r),schwartz_kl(0,0)) enddo do p = 1, ao_prim_num(i) @@ -175,19 +181,25 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l) 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 + 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) * & + coef2*coef2 + if (schwartz_kl(0,0)*schwartz_ij < thresh) then + cycle + endif do r = 1, ao_prim_num(k) + if (schwartz_kl(0,r)*schwartz_ij < thresh) then + cycle + endif 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 + if (schwartz_kl(s,r)*schwartz_ij < thresh) then cycle endif + coef4 = coef3*ao_coef_transp(s,l) 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), & @@ -212,14 +224,21 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l) enddo double precision :: ERI + schwartz_kl(0,0) = 0.d0 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) - 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),& 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)) + 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 + schwartz_kl(0,0) = max(schwartz_kl(0,r),schwartz_kl(0,0)) enddo 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),& 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)) + I_power(3),J_power(3),I_power(3),J_power(3))*coef2*coef2 + if (schwartz_kl(0,0)*schwartz_ij < thresh) then + cycle + 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) - coef4 = coef3*ao_coef_transp(s,l) - if (schwartz_kl(r,s)*schwartz_ij*coef4*coef4 < thresh) then + if (schwartz_kl(s,r)*schwartz_ij < thresh) then cycle endif + coef4 = coef3*ao_coef_transp(s,l) 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), & @@ -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 q = delta + gama 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)) + !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 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) implicit none 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 : ! 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 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