double precision function ao_tc_sym_two_e_pot(i,j,k,l) implicit none BEGIN_DOC ! integral of the AO basis or (ij|kl) ! i(r1) j(r1) (tc_pot(r12,mu)) k(r2) l(r2) ! ! where (tc_pot(r12,mu)) is the scalar part of the potential EXCLUDING the term erf(mu r12)/r12. ! ! See Eq. (32) of JCP 154, 084119 (2021). 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 'utils/constants.include.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 double precision :: scw_gauss_int,general_primitive_integral_gauss 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_tc_sym_two_e_pot = 0.d0 double precision :: thr thr = ao_integrals_threshold*ao_integrals_threshold allocate(schwartz_kl(0:ao_prim_num(l),0:ao_prim_num(k))) double precision :: coef3 double precision :: coef2 double precision :: p_inv,q_inv double precision :: coef1 double precision :: coef4 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 schwartz_kl(0,0) = 0.d0 do r = 1, ao_prim_num(k) coef1 = ao_coef_normalized_ordered_transp(r,k)*ao_coef_normalized_ordered_transp(r,k) schwartz_kl(0,r) = 0.d0 do s = 1, ao_prim_num(l) coef2 = coef1 * ao_coef_normalized_ordered_transp(s,l) * 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 scw_gauss_int = general_primitive_integral_gauss(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(s,r) = dabs(scw_gauss_int * 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) 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 scw_gauss_int = general_primitive_integral_gauss(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 = dabs(scw_gauss_int * coef2*coef2) if (schwartz_kl(0,0)*schwartz_ij < thr) then cycle endif do r = 1, ao_prim_num(k) if (schwartz_kl(0,r)*schwartz_ij < thr) then cycle endif coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) do s = 1, ao_prim_num(l) if (schwartz_kl(s,r)*schwartz_ij < thr) then cycle endif 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_gauss(dim1, & P_new,P_center,fact_p,pp,p_inv,iorder_p, & Q_new,Q_center,fact_q,qq,q_inv,iorder_q) ao_tc_sym_two_e_pot = ao_tc_sym_two_e_pot + coef4 * integral enddo ! s enddo ! r enddo ! q enddo ! p deallocate (schwartz_kl) end double precision function general_primitive_integral_gauss(dim, & P_new,P_center,fact_p,p,p_inv,iorder_p, & Q_new,Q_center,fact_q,q,q_inv,iorder_q) implicit none BEGIN_DOC ! Computes the integral where p,q,r,s are Gaussian primitives END_DOC integer,intent(in) :: dim include 'utils/constants.include.F' double precision, intent(in) :: P_new(0:max_dim,3),P_center(3),fact_p,p,p_inv double precision, intent(in) :: Q_new(0:max_dim,3),Q_center(3),fact_q,q,q_inv integer, intent(in) :: iorder_p(3) integer, intent(in) :: iorder_q(3) double precision :: r_cut,gama_r_cut,rho,dist double precision :: dx(0:max_dim),Ix_pol(0:max_dim),dy(0:max_dim),Iy_pol(0:max_dim),dz(0:max_dim),Iz_pol(0:max_dim) integer :: n_Ix,n_Iy,n_Iz,nx,ny,nz double precision :: bla integer :: ix,iy,iz,jx,jy,jz,i 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 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 :: thr thr = ao_integrals_threshold general_primitive_integral_gauss = 0.d0 !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: dx,Ix_pol,dy,Iy_pol,dz,Iz_pol !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: d1, d_poly ! Gaussian Product ! ---------------- pq = p_inv*0.5d0*q_inv pq_inv = 0.5d0/(p+q) p10_1 = q*pq ! 1/(2p) p01_1 = p*pq ! 1/(2q) pq_inv_2 = pq_inv+pq_inv p10_2 = pq_inv_2 * p10_1*q !0.5d0*q/(pq + p*p) p01_2 = pq_inv_2 * p01_1*p !0.5d0*p/(q*q + pq) accu = 0.d0 iorder = iorder_p(1)+iorder_q(1)+iorder_p(1)+iorder_q(1) do ix=0,iorder Ix_pol(ix) = 0.d0 enddo n_Ix = 0 do ix = 0, iorder_p(1) if (abs(P_new(ix,1)) < thr) cycle a = P_new(ix,1) do jx = 0, iorder_q(1) d = a*Q_new(jx,1) if (abs(d) < thr) cycle !DIR$ FORCEINLINE call give_polynom_mult_center_x(P_center(1),Q_center(1),ix,jx,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dx,nx) !DIR$ FORCEINLINE call add_poly_multiply(dx,nx,d,Ix_pol,n_Ix) enddo enddo if (n_Ix == -1) then return endif iorder = iorder_p(2)+iorder_q(2)+iorder_p(2)+iorder_q(2) do ix=0, iorder Iy_pol(ix) = 0.d0 enddo n_Iy = 0 do iy = 0, iorder_p(2) if (abs(P_new(iy,2)) > thr) then b = P_new(iy,2) do jy = 0, iorder_q(2) e = b*Q_new(jy,2) if (abs(e) < thr) cycle !DIR$ FORCEINLINE call give_polynom_mult_center_x(P_center(2),Q_center(2),iy,jy,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dy,ny) !DIR$ FORCEINLINE call add_poly_multiply(dy,ny,e,Iy_pol,n_Iy) enddo endif enddo if (n_Iy == -1) then return endif iorder = iorder_p(3)+iorder_q(3)+iorder_p(3)+iorder_q(3) do ix=0,iorder Iz_pol(ix) = 0.d0 enddo n_Iz = 0 do iz = 0, iorder_p(3) if (abs(P_new(iz,3)) > thr) then c = P_new(iz,3) do jz = 0, iorder_q(3) f = c*Q_new(jz,3) if (abs(f) < thr) cycle !DIR$ FORCEINLINE call give_polynom_mult_center_x(P_center(3),Q_center(3),iz,jz,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dz,nz) !DIR$ FORCEINLINE call add_poly_multiply(dz,nz,f,Iz_pol,n_Iz) enddo endif enddo if (n_Iz == -1) then return endif rho = p*q *pq_inv_2 dist = (P_center(1) - Q_center(1))*(P_center(1) - Q_center(1)) + & (P_center(2) - Q_center(2))*(P_center(2) - Q_center(2)) + & (P_center(3) - Q_center(3))*(P_center(3) - Q_center(3)) const = dist*rho n_pt_tmp = n_Ix+n_Iy do i=0,n_pt_tmp d_poly(i)=0.d0 enddo !DIR$ FORCEINLINE 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 do i=0,n_pt_out d1(i)=0.d0 enddo !DIR$ FORCEINLINE call multiply_poly(d_poly ,n_pt_tmp ,Iz_pol,n_Iz,d1,n_pt_out) double precision :: aa,c_a,t_a,rho_old,w_a,pi_3,prefactor,inv_pq_3_2 double precision :: gauss_int integer :: m gauss_int = 0.d0 pi_3 = pi*pi*pi inv_pq_3_2 = (p_inv * q_inv)**(1.5d0) rho_old = (p*q)/(p+q) prefactor = pi_3 * inv_pq_3_2 * fact_p * fact_q do i = 1, n_gauss_eff_pot ! browse the gaussians with different expo/coef !do i = 1, n_gauss_eff_pot-1 aa = expo_gauss_eff_pot(i) c_a = coef_gauss_eff_pot(i) t_a = dsqrt( aa /(rho_old + aa) ) w_a = dexp(-t_a*t_a*rho_old*dist) accu = 0.d0 ! evaluation of the polynom Ix(t_a) * Iy(t_a) * Iz(t_a) do m = 0, n_pt_out,2 accu += d1(m) * (t_a)**(dble(m)) enddo ! equation A8 of PRA-70-062505 (2004) of Toul. Col. Sav. gauss_int = gauss_int + c_a * prefactor * (1.d0 - t_a*t_a)**(1.5d0) * w_a * accu enddo general_primitive_integral_gauss = gauss_int end subroutine compute_ao_integrals_gauss_jl(j,l,n_integrals,buffer_i,buffer_value) implicit none use map_module BEGIN_DOC ! Parallel client for AO integrals END_DOC integer, intent(in) :: j,l integer,intent(out) :: n_integrals integer(key_kind),intent(out) :: buffer_i(ao_num*ao_num) real(integral_kind),intent(out) :: buffer_value(ao_num*ao_num) integer :: i,k double precision :: cpu_1,cpu_2, wall_1, wall_2 double precision :: integral, wall_0 double precision :: thr,ao_tc_sym_two_e_pot integer :: kk, m, j1, i1 logical, external :: ao_two_e_integral_zero thr = ao_integrals_threshold n_integrals = 0 j1 = j+ishft(l*l-l,-1) do k = 1, ao_num ! r1 i1 = ishft(k*k-k,-1) if (i1 > j1) then exit endif do i = 1, k i1 += 1 if (i1 > j1) then exit endif ! if (ao_two_e_integral_zero(i,j,k,l)) then if (.False.) then cycle endif if (ao_two_e_integral_erf_schwartz(i,k)*ao_two_e_integral_erf_schwartz(j,l) < thr ) then cycle endif !DIR$ FORCEINLINE integral = ao_tc_sym_two_e_pot(i,k,j,l) ! i,k : r1 j,l : r2 if (abs(integral) < thr) then cycle endif n_integrals += 1 !DIR$ FORCEINLINE call two_e_integrals_index(i,j,k,l,buffer_i(n_integrals)) buffer_value(n_integrals) = integral enddo enddo end