diff --git a/plugins/Integrals_erf/EZFIO.cfg b/plugins/Integrals_erf/EZFIO.cfg new file mode 100644 index 00000000..fb08245f --- /dev/null +++ b/plugins/Integrals_erf/EZFIO.cfg @@ -0,0 +1,84 @@ +[do_direct_integrals] +type: logical +doc: Compute integrals on the fly +interface: ezfio,provider,ocaml +default: False +ezfio_name: direct + +[no_vvvv_integrals] +type: logical +doc: If True, computes all integrals except for the integrals having 4 virtual index +interface: ezfio,provider,ocaml +default: False +ezfio_name: no_vvvv_integrals + +[no_ivvv_integrals] +type: logical +doc: Can be switched on only if no_vvvv_integrals is True, then do not computes the integrals having 3 virtual index and 1 belonging to the core inactive active orbitals +interface: ezfio,provider,ocaml +default: False +ezfio_name: no_ivvv_integrals + +[no_vvv_integrals] +type: logical +doc: Can be switched on only if no_vvvv_integrals is True, then do not computes the integrals having 3 virtual orbitals +interface: ezfio,provider,ocaml +default: False +ezfio_name: no_vvv_integrals + +[disk_access_mo_integrals] +type: Disk_access +doc: Read/Write MO integrals from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + +[disk_access_ao_integrals_standard] +type: Disk_access +doc: Read/Write AO integrals_standard from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + + +[disk_access_mo_integrals_standard] +type: Disk_access +doc: Read/Write MO integrals_standard from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + +[disk_access_ao_integrals] +type: Disk_access +doc: Read/Write AO integrals from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + + + +[ao_integrals_threshold] +type: Threshold +doc: If || < ao_integrals_threshold then is zero +interface: ezfio,provider,ocaml +default: 1.e-15 +ezfio_name: threshold_ao + +[mo_integrals_threshold] +type: Threshold +doc: If || < ao_integrals_threshold then is zero +interface: ezfio,provider,ocaml +default: 1.e-15 +ezfio_name: threshold_mo + +[mu_erf] +type: double precision +doc: cutting of the interaction in the range separated model +interface: ezfio,provider,ocaml +default: 0.5 +ezfio_name: mu_erf + + +[long_range] +type: logical +doc: if true, compute all the integrals using the long range interaction +interface: ezfio,provider,ocaml +default: False +ezfio_name: long_range + diff --git a/plugins/Integrals_erf/NEEDED_CHILDREN_MODULES b/plugins/Integrals_erf/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..152711f3 --- /dev/null +++ b/plugins/Integrals_erf/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Pseudo Bitmask ZMQ diff --git a/plugins/Integrals_erf/ao_bi_integrals.irp.f b/plugins/Integrals_erf/ao_bi_integrals.irp.f new file mode 100644 index 00000000..b4c9adb5 --- /dev/null +++ b/plugins/Integrals_erf/ao_bi_integrals.irp.f @@ -0,0 +1,1121 @@ +double precision function ao_bielec_integral(i,j,k,l) + implicit none + BEGIN_DOC + ! integral of the AO basis 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 '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 :: 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 + + 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 = 0.d0 + + 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 + double precision :: p_inv,q_inv + double precision :: general_primitive_integral + + 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 = 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 = ao_bielec_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, 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) + 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) + integral = ERI( & + ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_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 = ao_bielec_integral + coef4 * integral + enddo ! s + enddo ! r + enddo ! q + enddo ! p + + endif + +end + +double precision function ao_bielec_integral_schwartz_accel(i,j,k,l) + implicit none + BEGIN_DOC + ! integral of the AO basis 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 '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 + + 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 :: thr + thr = ao_integrals_threshold*ao_integrals_threshold + + 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 + 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 + 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) & + * 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) + double precision :: coef1 + coef1 = ao_coef_normalized_ordered_transp(p,i) + do q = 1, ao_prim_num(j) + double precision :: coef2 + coef2 = coef1*ao_coef_normalized_ordered_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_ordered_transp(p,i),ao_expo_ordered_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) * & + 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 + double precision :: coef3 + coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) + do s = 1, ao_prim_num(l) + double precision :: coef4 + if (schwartz_kl(s,r)*schwartz_ij < thr) then + cycle + endif + coef4 = coef3*ao_coef_normalized_ordered_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_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_bielec_integral_schwartz_accel = 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 + + 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) + schwartz_kl(s,r) = ERI( & + ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),ao_expo_ordered_transp(r,k),ao_expo_ordered_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)) * & + 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) + schwartz_ij = ERI( & + ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(p,i),ao_expo_ordered_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))*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) + integral = ERI( & + ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_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 = 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) + implicit none + BEGIN_DOC +! Computes the product of l values of i,j,k,and l + END_DOC + integer, intent(in) :: i,j,k,l + ao_l4 = ao_l(i)*ao_l(j)*ao_l(k)*ao_l(l) +end + + + +subroutine compute_ao_bielec_integrals(j,k,l,sze,buffer_value) + implicit none + use map_module + + BEGIN_DOC + ! Compute AO 1/r12 integrals for all i and fixed j,k,l + END_DOC + + include 'Utils/constants.include.F' + integer, intent(in) :: j,k,l,sze + real(integral_kind), intent(out) :: buffer_value(sze) + double precision :: ao_bielec_integral + + integer :: i + + if (ao_overlap_abs(j,l) < thresh) then + buffer_value = 0._integral_kind + return + endif + if (ao_bielec_integral_schwartz(j,l) < thresh ) then + buffer_value = 0._integral_kind + return + endif + + do i = 1, ao_num + if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then + buffer_value(i) = 0._integral_kind + cycle + endif + if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh ) then + buffer_value(i) = 0._integral_kind + cycle + endif + !DIR$ FORCEINLINE + buffer_value(i) = ao_bielec_integral(i,k,j,l) + enddo + +end + +double precision function general_primitive_integral(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) + + general_primitive_integral = 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 + ! ---------------- + double precision :: p_plus_q + p_plus_q = (p+q) + pq = p_inv*0.5d0*q_inv + + pq_inv = 0.5d0/p_plus_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) + !DIR$ VECTOR ALIGNED + 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)) < thresh) cycle + a = P_new(ix,1) + do jx = 0, iorder_q(1) + d = a*Q_new(jx,1) + if (abs(d) < thresh) cycle + !DEC$ 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) + !DEC$ 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) + !DIR$ VECTOR ALIGNED + 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)) > thresh) then + b = P_new(iy,2) + do jy = 0, iorder_q(2) + e = b*Q_new(jy,2) + if (abs(e) < thresh) cycle + !DEC$ 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) + !DEC$ 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)) > thresh) then + c = P_new(iz,3) + do jz = 0, iorder_q(3) + f = c*Q_new(jz,3) + if (abs(f) < thresh) cycle + !DEC$ 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) + !DEC$ 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 ! le rho qui va bien + 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 + + !DEC$ 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 + + !DEC$ FORCEINLINE + call multiply_poly(d_poly ,n_pt_tmp ,Iz_pol,n_Iz,d1,n_pt_out) + double precision :: rint_sum + accu = accu + rint_sum(n_pt_out,const,d1) + + ! change p+q in dsqrt + general_primitive_integral = fact_p * fact_q * accu *pi_5_2*p_inv*q_inv/dsqrt(p_plus_q) +end + + +double precision function ERI(alpha,beta,delta,gama,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 + ! ATOMIC PRIMTIVE bielectronic integral between the 4 primitives :: + ! primitive_1 = x1**(a_x) y1**(a_y) z1**(a_z) exp(-alpha * r1**2) + ! primitive_2 = x1**(b_x) y1**(b_y) z1**(b_z) exp(- beta * r1**2) + ! primitive_3 = x2**(c_x) y2**(c_y) z2**(c_z) exp(-delta * r2**2) + ! primitive_4 = x2**(d_x) y2**(d_y) z2**(d_z) exp(- gama * r2**2) + END_DOC + double precision, intent(in) :: delta,gama,alpha,beta + integer, intent(in) :: 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_2,b_x_2,c_x_2,d_x_2,a_y_2,b_y_2,c_y_2,d_y_2,a_z_2,b_z_2,c_z_2,d_z_2 + integer :: i,j,k,l,n_pt + integer :: n_pt_sup + double precision :: p,q,denom,coeff + double precision :: I_f + integer :: nx,ny,nz + include 'Utils/constants.include.F' + nx = a_x+b_x+c_x+d_x + if(iand(nx,1) == 1) then + ERI = 0.d0 + return + endif + + ny = a_y+b_y+c_y+d_y + if(iand(ny,1) == 1) then + ERI = 0.d0 + return + endif + + nz = a_z+b_z+c_z+d_z + if(iand(nz,1) == 1) then + ERI = 0.d0 + return + endif + + ASSERT (alpha >= 0.d0) + ASSERT (beta >= 0.d0) + ASSERT (delta >= 0.d0) + ASSERT (gama >= 0.d0) + p = alpha + beta + q = delta + gama + double precision :: p_plus_q + p_plus_q = (p+q) + ASSERT (p+q >= 0.d0) + n_pt = ishft( nx+ny+nz,1 ) + + coeff = pi_5_2 / (p * q * dsqrt(p_plus_q)) + if (n_pt == 0) then + ERI = coeff + return + endif + + call integrale_new(I_f,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z,p,q,n_pt) + + ERI = I_f * coeff +end + + +subroutine integrale_new(I_f,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z,p,q,n_pt) + BEGIN_DOC + ! calculate the integral of the polynom :: + ! I_x1(a_x+b_x, c_x+d_x,p,q) * I_x1(a_y+b_y, c_y+d_y,p,q) * I_x1(a_z+b_z, c_z+d_z,p,q) + ! between ( 0 ; 1) + END_DOC + + + implicit none + include 'Utils/constants.include.F' + double precision :: p,q + 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 :: i, n_iter, n_pt, j + double precision :: I_f, pq_inv, p10_1, p10_2, p01_1, p01_2,rho,pq_inv_2 + integer :: ix,iy,iz, jx,jy,jz, sx,sy,sz + double precision :: p_plus_q + + j = ishft(n_pt,-1) + ASSERT (n_pt > 1) + p_plus_q = (p+q) + pq_inv = 0.5d0/(p_plus_q) + pq_inv_2 = pq_inv + pq_inv + p10_1 = 0.5d0/p + p01_1 = 0.5d0/q + p10_2 = 0.5d0 * q /(p * p_plus_q) + p01_2 = 0.5d0 * p /(q * p_plus_q) + double precision :: B00(n_pt_max_integrals) + double precision :: B10(n_pt_max_integrals), B01(n_pt_max_integrals) + double precision :: t1(n_pt_max_integrals), t2(n_pt_max_integrals) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: t1, t2, B10, B01, B00 + ix = a_x+b_x + jx = c_x+d_x + iy = a_y+b_y + jy = c_y+d_y + iz = a_z+b_z + jz = c_z+d_z + sx = ix+jx + sy = iy+jy + sz = iz+jz + + !DIR$ VECTOR ALIGNED + do i = 1,n_pt + B10(i) = p10_1 - gauleg_t2(i,j)* p10_2 + B01(i) = p01_1 - gauleg_t2(i,j)* p01_2 + B00(i) = gauleg_t2(i,j)*pq_inv + enddo + if (sx > 0) then + call I_x1_new(ix,jx,B10,B01,B00,t1,n_pt) + else + !DIR$ VECTOR ALIGNED + do i = 1,n_pt + t1(i) = 1.d0 + enddo + endif + if (sy > 0) then + call I_x1_new(iy,jy,B10,B01,B00,t2,n_pt) + !DIR$ VECTOR ALIGNED + do i = 1,n_pt + t1(i) = t1(i)*t2(i) + enddo + endif + if (sz > 0) then + call I_x1_new(iz,jz,B10,B01,B00,t2,n_pt) + !DIR$ VECTOR ALIGNED + do i = 1,n_pt + t1(i) = t1(i)*t2(i) + enddo + endif + I_f= 0.d0 + !DIR$ VECTOR ALIGNED + do i = 1,n_pt + I_f += gauleg_w(i,j)*t1(i) + enddo + + + +end + +recursive subroutine I_x1_new(a,c,B_10,B_01,B_00,res,n_pt) + BEGIN_DOC + ! recursive function involved in the bielectronic integral + END_DOC + implicit none + include 'Utils/constants.include.F' + integer, intent(in) :: a,c,n_pt + double precision, intent(in) :: B_10(n_pt_max_integrals),B_01(n_pt_max_integrals),B_00(n_pt_max_integrals) + double precision, intent(out) :: res(n_pt_max_integrals) + double precision :: res2(n_pt_max_integrals) + integer :: i + + if(c<0)then + !DIR$ VECTOR ALIGNED + do i=1,n_pt + res(i) = 0.d0 + enddo + else if (a==0) then + call I_x2_new(c,B_10,B_01,B_00,res,n_pt) + else if (a==1) then + call I_x2_new(c-1,B_10,B_01,B_00,res,n_pt) + !DIR$ VECTOR ALIGNED + do i=1,n_pt + res(i) = c * B_00(i) * res(i) + enddo + else + call I_x1_new(a-2,c,B_10,B_01,B_00,res,n_pt) + call I_x1_new(a-1,c-1,B_10,B_01,B_00,res2,n_pt) + !DIR$ VECTOR ALIGNED + do i=1,n_pt + res(i) = (a-1) * B_10(i) * res(i) & + + c * B_00(i) * res2(i) + enddo + endif +end + +recursive subroutine I_x2_new(c,B_10,B_01,B_00,res,n_pt) + implicit none + BEGIN_DOC + ! recursive function involved in the bielectronic integral + END_DOC + include 'Utils/constants.include.F' + integer, intent(in) :: c, n_pt + double precision, intent(in) :: B_10(n_pt_max_integrals),B_01(n_pt_max_integrals),B_00(n_pt_max_integrals) + double precision, intent(out) :: res(n_pt_max_integrals) + integer :: i + + if(c==1)then + !DIR$ VECTOR ALIGNED + do i=1,n_pt + res(i) = 0.d0 + enddo + elseif(c==0) then + !DIR$ VECTOR ALIGNED + do i=1,n_pt + res(i) = 1.d0 + enddo + else + call I_x1_new(0,c-2,B_10,B_01,B_00,res,n_pt) + !DIR$ VECTOR ALIGNED + do i=1,n_pt + res(i) = (c-1) * B_01(i) * res(i) + enddo + endif +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 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 = 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 + + + + +subroutine give_polynom_mult_center_x(P_center,Q_center,a_x,d_x,p,q,n_pt_in,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,d,n_pt_out) + implicit none + BEGIN_DOC + ! subroutine that returns the explicit polynom in term of the "t" + ! variable of the following polynomw : + ! I_x1(a_x, d_x,p,q) * I_x1(a_y, d_y,p,q) * I_x1(a_z, d_z,p,q) + END_DOC + integer, intent(in) :: n_pt_in + integer,intent(out) :: n_pt_out + integer, intent(in) :: a_x,d_x + double precision, intent(in) :: P_center, Q_center + double precision, intent(in) :: p,q,pq_inv,p10_1,p01_1,p10_2,p01_2,pq_inv_2 + include 'Utils/constants.include.F' + double precision,intent(out) :: d(0:max_dim) + double precision :: accu + accu = 0.d0 + ASSERT (n_pt_in >= 0) + ! pq_inv = 0.5d0/(p+q) + ! pq_inv_2 = 1.d0/(p+q) + ! p10_1 = 0.5d0/p + ! p01_1 = 0.5d0/q + ! p10_2 = 0.5d0 * q /(p * q + p * p) + ! p01_2 = 0.5d0 * p /(q * q + q * p) + double precision :: B10(0:2), B01(0:2), B00(0:2),C00(0:2),D00(0:2) + B10(0) = p10_1 + B10(1) = 0.d0 + B10(2) = - p10_2 + ! B10 = p01_1 - t**2 * p10_2 + B01(0) = p01_1 + B01(1) = 0.d0 + B01(2) = - p01_2 + ! B01 = p01_1- t**2 * pq_inv + B00(0) = 0.d0 + B00(1) = 0.d0 + B00(2) = pq_inv + ! B00 = t**2 * pq_inv + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + integer :: n_pt1,dim,i + n_pt1 = n_pt_in + ! C00 = -q/(p+q)*(Px-Qx) * t^2 + C00(0) = 0.d0 + C00(1) = 0.d0 + C00(2) = -q*(P_center-Q_center) * pq_inv_2 + ! D00 = -p/(p+q)*(Px-Qx) * t^2 + D00(0) = 0.d0 + D00(1) = 0.d0 + D00(2) = -p*(Q_center-P_center) * pq_inv_2 + !D00(2) = -p*(Q_center(1)-P_center(1)) /(p+q) + !DIR$ FORCEINLINE + call I_x1_pol_mult(a_x,d_x,B10,B01,B00,C00,D00,d,n_pt1,n_pt_in) + n_pt_out = n_pt1 + if(n_pt1<0)then + n_pt_out = -1 + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + return + endif + +end + +subroutine I_x1_pol_mult(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) + implicit none + BEGIN_DOC + ! recursive function involved in the bielectronic integral + END_DOC + integer , intent(in) :: n_pt_in + include 'Utils/constants.include.F' + double precision,intent(inout) :: d(0:max_dim) + integer,intent(inout) :: nd + integer, intent(in) :: a,c + double precision, intent(in) :: B_10(0:2),B_01(0:2),B_00(0:2),C_00(0:2),D_00(0:2) + if( (c>=0).and.(nd>=0) )then + + if (a==1) then + call I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) + else if (a==2) then + call I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) + else if (a>2) then + call I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) + else ! a == 0 + + if( c==0 )then + nd = 0 + d(0) = 1.d0 + return + endif + + call I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) + endif + else + nd = -1 + endif +end + +recursive subroutine I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) + implicit none + BEGIN_DOC + ! recursive function involved in the bielectronic integral + END_DOC + integer , intent(in) :: n_pt_in + include 'Utils/constants.include.F' + double precision,intent(inout) :: d(0:max_dim) + integer,intent(inout) :: nd + integer, intent(in) :: a,c + double precision, intent(in) :: B_10(0:2),B_01(0:2),B_00(0:2),C_00(0:2),D_00(0:2) + double precision :: X(0:max_dim) + double precision :: Y(0:max_dim) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X,Y + integer :: nx, ix,iy,ny + + ASSERT (a>2) + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,n_pt_in + X(ix) = 0.d0 + enddo + nx = 0 + if (a==3) then + call I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) + else if (a==4) then + call I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) + else + ASSERT (a>=5) + call I_x1_pol_mult_recurs(a-2,c,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) + endif + + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,nx + X(ix) *= dble(a-1) + enddo + + !DEC$ FORCEINLINE + call multiply_poly(X,nx,B_10,2,d,nd) + + nx = nd + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,n_pt_in + X(ix) = 0.d0 + enddo + + if (c>0) then + if (a==3) then + call I_x1_pol_mult_a2(c-1,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) + else + ASSERT(a >= 4) + call I_x1_pol_mult_recurs(a-1,c-1,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) + endif + if (c>1) then + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,nx + X(ix) *= c + enddo + endif + !DEC$ FORCEINLINE + call multiply_poly(X,nx,B_00,2,d,nd) + endif + + ny=0 + + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,n_pt_in + Y(ix) = 0.d0 + enddo + ASSERT(a > 2) + if (a==3) then + call I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,Y,ny,n_pt_in) + else + ASSERT(a >= 4) + call I_x1_pol_mult_recurs(a-1,c,B_10,B_01,B_00,C_00,D_00,Y,ny,n_pt_in) + endif + + !DEC$ FORCEINLINE + call multiply_poly(Y,ny,C_00,2,d,nd) + +end + +recursive subroutine I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) + implicit none + BEGIN_DOC + ! recursive function involved in the bielectronic integral + END_DOC + integer , intent(in) :: n_pt_in + include 'Utils/constants.include.F' + double precision,intent(inout) :: d(0:max_dim) + integer,intent(inout) :: nd + integer, intent(in) :: c + double precision, intent(in) :: B_10(0:2),B_01(0:2),B_00(0:2),C_00(0:2),D_00(0:2) + double precision :: X(0:max_dim) + double precision :: Y(0:max_dim) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X,Y + integer :: nx, ix,iy,ny + + if( (c<0).or.(nd<0) )then + nd = -1 + return + endif + + nx = nd + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,n_pt_in + X(ix) = 0.d0 + enddo + call I_x2_pol_mult(c-1,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) + + if (c>1) then + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,nx + X(ix) *= dble(c) + enddo + endif + + !DEC$ FORCEINLINE + call multiply_poly(X,nx,B_00,2,d,nd) + + ny=0 + + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,n_pt_in + Y(ix) = 0.d0 + enddo + call I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,Y,ny,n_pt_in) + + !DEC$ FORCEINLINE + call multiply_poly(Y,ny,C_00,2,d,nd) + +end + +recursive subroutine I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) + implicit none + BEGIN_DOC + ! recursive function involved in the bielectronic integral + END_DOC + integer , intent(in) :: n_pt_in + include 'Utils/constants.include.F' + double precision,intent(inout) :: d(0:max_dim) + integer,intent(inout) :: nd + integer, intent(in) :: c + double precision, intent(in) :: B_10(0:2),B_01(0:2),B_00(0:2),C_00(0:2),D_00(0:2) + double precision :: X(0:max_dim) + double precision :: Y(0:max_dim) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X,Y + integer :: nx, ix,iy,ny + + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,n_pt_in + X(ix) = 0.d0 + enddo + nx = 0 + call I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) + + !DEC$ FORCEINLINE + call multiply_poly(X,nx,B_10,2,d,nd) + + nx = nd + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,n_pt_in + X(ix) = 0.d0 + enddo + + !DIR$ FORCEINLINE + call I_x1_pol_mult_a1(c-1,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) + + if (c>1) then + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,nx + X(ix) *= dble(c) + enddo + endif + + !DEC$ FORCEINLINE + call multiply_poly(X,nx,B_00,2,d,nd) + + ny=0 + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,n_pt_in + Y(ix) = 0.d0 + enddo + !DEC$ FORCEINLINE + call I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,Y,ny,n_pt_in) + + !DEC$ FORCEINLINE + call multiply_poly(Y,ny,C_00,2,d,nd) + +end + +recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim) + implicit none + BEGIN_DOC + ! recursive function involved in the bielectronic integral + END_DOC + integer , intent(in) :: dim + include 'Utils/constants.include.F' + double precision :: d(0:max_dim) + integer,intent(inout) :: nd + integer, intent(in) :: c + double precision, intent(in) :: B_10(0:2),B_01(0:2),B_00(0:2),C_00(0:2),D_00(0:2) + integer :: nx, ix,ny + double precision :: X(0:max_dim),Y(0:max_dim) + !DEC$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y + integer :: i + + select case (c) + case (0) + nd = 0 + d(0) = 1.d0 + return + + case (:-1) + nd = -1 + return + + case (1) + nd = 2 + d(0) = D_00(0) + d(1) = D_00(1) + d(2) = D_00(2) + return + + case (2) + nd = 2 + d(0) = B_01(0) + d(1) = B_01(1) + d(2) = B_01(2) + + ny = 2 + Y(0) = D_00(0) + Y(1) = D_00(1) + Y(2) = D_00(2) + + !DEC$ FORCEINLINE + call multiply_poly(Y,ny,D_00,2,d,nd) + return + + case default + + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(6) + do ix=0,c+c + X(ix) = 0.d0 + enddo + nx = 0 + call I_x2_pol_mult(c-2,B_10,B_01,B_00,C_00,D_00,X,nx,dim) + + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(6) + do ix=0,nx + X(ix) *= dble(c-1) + enddo + + !DEC$ FORCEINLINE + call multiply_poly(X,nx,B_01,2,d,nd) + + ny = 0 + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(6) + do ix=0,c+c + Y(ix) = 0.d0 + enddo + call I_x2_pol_mult(c-1,B_10,B_01,B_00,C_00,D_00,Y,ny,dim) + + !DEC$ FORCEINLINE + call multiply_poly(Y,ny,D_00,2,d,nd) + + end select +end + + + + +subroutine compute_ao_integrals_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 :: ao_bielec_integral,cpu_1,cpu_2, wall_1, wall_2 + double precision :: integral, wall_0 + double precision :: thr + integer :: kk, m, j1, i1 + + 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_overlap_abs(i,k)*ao_overlap_abs(j,l) < thr) then + cycle + endif + if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thr ) then + cycle + endif + !DIR$ FORCEINLINE + integral = ao_bielec_integral(i,k,j,l) ! i,k : r1 j,l : r2 + if (abs(integral) < thr) then + cycle + endif + n_integrals += 1 + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,buffer_i(n_integrals)) + buffer_value(n_integrals) = integral + enddo + enddo + +end diff --git a/plugins/Integrals_erf/ao_bi_integrals_erf.irp.f b/plugins/Integrals_erf/ao_bi_integrals_erf.irp.f new file mode 100644 index 00000000..6616883f --- /dev/null +++ b/plugins/Integrals_erf/ao_bi_integrals_erf.irp.f @@ -0,0 +1,570 @@ +double precision function ao_bielec_integral_erf(i,j,k,l) + implicit none + BEGIN_DOC + ! integral of the AO basis 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 '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 :: ao_bielec_integral_schwartz_accel_erf + + if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then + ao_bielec_integral_erf = ao_bielec_integral_schwartz_accel_erf(i,j,k,l) + return + endif + + 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_erf = 0.d0 + + 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 + double precision :: p_inv,q_inv + double precision :: general_primitive_integral_erf + + 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 = general_primitive_integral_erf(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_erf = ao_bielec_integral_erf + 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_erf + + 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) + 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) + integral = ERI_erf( & + ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_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_erf = ao_bielec_integral_erf + coef4 * integral + enddo ! s + enddo ! r + enddo ! q + enddo ! p + + endif + +end + +double precision function ao_bielec_integral_schwartz_accel_erf(i,j,k,l) + implicit none + BEGIN_DOC + ! integral of the AO basis 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 '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 + + 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_erf = 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))) + + + 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 + + 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 + schwartz_kl(s,r) = general_primitive_integral_erf(dim1, & + Q_new,Q_center,fact_q,qq,q_inv,iorder_q, & + 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) + double precision :: coef1 + coef1 = ao_coef_normalized_ordered_transp(p,i) + do q = 1, ao_prim_num(j) + double precision :: coef2 + coef2 = coef1*ao_coef_normalized_ordered_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_ordered_transp(p,i),ao_expo_ordered_transp(q,j), & + I_power,J_power,I_center,J_center,dim1) + p_inv = 1.d0/pp + schwartz_ij = general_primitive_integral_erf(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 < thr) then + cycle + endif + do r = 1, ao_prim_num(k) + if (schwartz_kl(0,r)*schwartz_ij < thr) then + cycle + endif + double precision :: coef3 + coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) + do s = 1, ao_prim_num(l) + double precision :: coef4 + if (schwartz_kl(s,r)*schwartz_ij < thr) then + cycle + endif + coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) + double precision :: general_primitive_integral_erf + 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_erf(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_erf = ao_bielec_integral_schwartz_accel_erf + 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_erf + + 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) + schwartz_kl(s,r) = ERI_erf( & + ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),ao_expo_ordered_transp(r,k),ao_expo_ordered_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)) * & + 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) + schwartz_ij = ERI_erf( & + ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(p,i),ao_expo_ordered_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))*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) + integral = ERI_erf( & + ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_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_erf = ao_bielec_integral_schwartz_accel_erf + coef4 * integral + enddo ! s + enddo ! r + enddo ! q + enddo ! p + + endif + deallocate (schwartz_kl) + +end + + +subroutine compute_ao_bielec_integrals_erf(j,k,l,sze,buffer_value) + implicit none + use map_module + + BEGIN_DOC + ! Compute AO 1/r12 integrals for all i and fixed j,k,l + END_DOC + + include 'Utils/constants.include.F' + integer, intent(in) :: j,k,l,sze + real(integral_kind), intent(out) :: buffer_value(sze) + double precision :: ao_bielec_integral_erf + + integer :: i + + if (ao_overlap_abs(j,l) < thresh) then + buffer_value = 0._integral_kind + return + endif + if (ao_bielec_integral_schwartz(j,l) < thresh ) then + buffer_value = 0._integral_kind + return + endif + + do i = 1, ao_num + if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then + buffer_value(i) = 0._integral_kind + cycle + endif + if (ao_bielec_integral_erf_schwartz(i,k)*ao_bielec_integral_erf_schwartz(j,l) < thresh ) then + buffer_value(i) = 0._integral_kind + cycle + endif + !DIR$ FORCEINLINE + buffer_value(i) = ao_bielec_integral_erf(i,k,j,l) + enddo + +end + +double precision function general_primitive_integral_erf(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) + + general_primitive_integral_erf = 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 + ! ---------------- + double precision :: p_plus_q + p_plus_q = (p+q) * ((p*q)/(p+q) + mu_erf*mu_erf)/(mu_erf*mu_erf) + pq = p_inv*0.5d0*q_inv + + pq_inv = 0.5d0/p_plus_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) + !DIR$ VECTOR ALIGNED + 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)) < thresh) cycle + a = P_new(ix,1) + do jx = 0, iorder_q(1) + d = a*Q_new(jx,1) + if (abs(d) < thresh) cycle + !DEC$ 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) + !DEC$ 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) + !DIR$ VECTOR ALIGNED + 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)) > thresh) then + b = P_new(iy,2) + do jy = 0, iorder_q(2) + e = b*Q_new(jy,2) + if (abs(e) < thresh) cycle + !DEC$ 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) + !DEC$ 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)) > thresh) then + c = P_new(iz,3) + do jz = 0, iorder_q(3) + f = c*Q_new(jz,3) + if (abs(f) < thresh) cycle + !DEC$ 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) + !DEC$ 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 ! le rho qui va bien + 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 + + !DEC$ 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 + + !DEC$ FORCEINLINE + call multiply_poly(d_poly ,n_pt_tmp ,Iz_pol,n_Iz,d1,n_pt_out) + double precision :: rint_sum + accu = accu + rint_sum(n_pt_out,const,d1) + + ! change p+q in dsqrt + general_primitive_integral_erf = fact_p * fact_q * accu *pi_5_2*p_inv*q_inv/dsqrt(p_plus_q) +end + + +double precision function ERI_erf(alpha,beta,delta,gama,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 + ! ATOMIC PRIMTIVE bielectronic integral between the 4 primitives :: + ! primitive_1 = x1**(a_x) y1**(a_y) z1**(a_z) exp(-alpha * r1**2) + ! primitive_2 = x1**(b_x) y1**(b_y) z1**(b_z) exp(- beta * r1**2) + ! primitive_3 = x2**(c_x) y2**(c_y) z2**(c_z) exp(-delta * r2**2) + ! primitive_4 = x2**(d_x) y2**(d_y) z2**(d_z) exp(- gama * r2**2) + END_DOC + double precision, intent(in) :: delta,gama,alpha,beta + integer, intent(in) :: 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_2,b_x_2,c_x_2,d_x_2,a_y_2,b_y_2,c_y_2,d_y_2,a_z_2,b_z_2,c_z_2,d_z_2 + integer :: i,j,k,l,n_pt + integer :: n_pt_sup + double precision :: p,q,denom,coeff + double precision :: I_f + integer :: nx,ny,nz + include 'Utils/constants.include.F' + nx = a_x+b_x+c_x+d_x + if(iand(nx,1) == 1) then + ERI_erf = 0.d0 + return + endif + + ny = a_y+b_y+c_y+d_y + if(iand(ny,1) == 1) then + ERI_erf = 0.d0 + return + endif + + nz = a_z+b_z+c_z+d_z + if(iand(nz,1) == 1) then + ERI_erf = 0.d0 + return + endif + + ASSERT (alpha >= 0.d0) + ASSERT (beta >= 0.d0) + ASSERT (delta >= 0.d0) + ASSERT (gama >= 0.d0) + p = alpha + beta + q = delta + gama + double precision :: p_plus_q + p_plus_q = (p+q) * ((p*q)/(p+q) + mu_erf*mu_erf)/(mu_erf*mu_erf) + ASSERT (p+q >= 0.d0) + n_pt = ishft( nx+ny+nz,1 ) + + coeff = pi_5_2 / (p * q * dsqrt(p_plus_q)) + if (n_pt == 0) then + ERI_erf = coeff + return + endif + + call integrale_new(I_f,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z,p,q,n_pt) + + ERI_erf = I_f * coeff +end + + + + +subroutine compute_ao_integrals_erf_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 :: ao_bielec_integral_erf,cpu_1,cpu_2, wall_1, wall_2 + double precision :: integral, wall_0 + double precision :: thr + integer :: kk, m, j1, i1 + + 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_overlap_abs(i,k)*ao_overlap_abs(j,l) < thr) then + cycle + endif + if (ao_bielec_integral_erf_schwartz(i,k)*ao_bielec_integral_erf_schwartz(j,l) < thr ) then + cycle + endif + !DIR$ FORCEINLINE + integral = ao_bielec_integral_erf(i,k,j,l) ! i,k : r1 j,l : r2 + if (abs(integral) < thr) then + cycle + endif + n_integrals += 1 + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,buffer_i(n_integrals)) + buffer_value(n_integrals) = integral + enddo + enddo + +end diff --git a/plugins/Integrals_erf/ao_bielec_integrals_erf_in_map_slave.irp.f b/plugins/Integrals_erf/ao_bielec_integrals_erf_in_map_slave.irp.f new file mode 100644 index 00000000..36f0e492 --- /dev/null +++ b/plugins/Integrals_erf/ao_bielec_integrals_erf_in_map_slave.irp.f @@ -0,0 +1,175 @@ +subroutine ao_bielec_integrals_erf_in_map_slave_tcp(i) + implicit none + integer, intent(in) :: i + BEGIN_DOC +! Computes a buffer of integrals. i is the ID of the current thread. + END_DOC + call ao_bielec_integrals_erf_in_map_slave(0,i) +end + + +subroutine ao_bielec_integrals_erf_in_map_slave_inproc(i) + implicit none + integer, intent(in) :: i + BEGIN_DOC +! Computes a buffer of integrals. i is the ID of the current thread. + END_DOC + call ao_bielec_integrals_erf_in_map_slave(1,i) +end + + + +subroutine ao_bielec_integrals_erf_in_map_slave(thread,iproc) + use map_module + use f77_zmq + implicit none + BEGIN_DOC +! Computes a buffer of integrals + END_DOC + + integer, intent(in) :: thread, iproc + + integer :: j,l,n_integrals + integer :: rc + real(integral_kind), allocatable :: buffer_value(:) + integer(key_kind), allocatable :: buffer_i(:) + + integer :: worker_id, task_id + character*(512) :: task + + integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + + integer(ZMQ_PTR), external :: new_zmq_push_socket + integer(ZMQ_PTR) :: zmq_socket_push + + character*(64) :: state + + zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() + zmq_socket_push = new_zmq_push_socket(thread) + + allocate ( buffer_i(ao_num*ao_num), buffer_value(ao_num*ao_num) ) + + call connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) + + do + call get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task) + if (task_id == 0) exit + read(task,*) j, l + call compute_ao_integrals_erf_jl(j,l,n_integrals,buffer_i,buffer_value) + call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id) + call push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value, task_id) + enddo + + + call disconnect_from_taskserver(zmq_to_qp_run_socket,zmq_socket_push,worker_id) + deallocate( buffer_i, buffer_value ) + call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) + call end_zmq_push_socket(zmq_socket_push,thread) + +end + + +subroutine ao_bielec_integrals_erf_in_map_collector + use map_module + use f77_zmq + implicit none + BEGIN_DOC +! Collects results from the AO integral calculation + END_DOC + + integer :: j,l,n_integrals + integer :: rc + + real(integral_kind), allocatable :: buffer_value(:) + integer(key_kind), allocatable :: buffer_i(:) + + integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + + integer(ZMQ_PTR), external :: new_zmq_pull_socket + integer(ZMQ_PTR) :: zmq_socket_pull + + integer*8 :: control, accu + integer :: task_id, more, sze + + zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() + zmq_socket_pull = new_zmq_pull_socket() + + sze = ao_num*ao_num + allocate ( buffer_i(sze), buffer_value(sze) ) + + accu = 0_8 + more = 1 + do while (more == 1) + + rc = f77_zmq_recv( zmq_socket_pull, n_integrals, 4, 0) + if (rc == -1) then + n_integrals = 0 + return + endif + if (rc /= 4) then + print *, irp_here, ': f77_zmq_recv( zmq_socket_pull, n_integrals, 4, 0)' + stop 'error' + endif + + if (n_integrals >= 0) then + + if (n_integrals > sze) then + deallocate (buffer_value, buffer_i) + sze = n_integrals + allocate (buffer_value(sze), buffer_i(sze)) + endif + + rc = f77_zmq_recv( zmq_socket_pull, buffer_i, key_kind*n_integrals, 0) + if (rc /= key_kind*n_integrals) then + print *, rc, key_kind, n_integrals + print *, irp_here, ': f77_zmq_recv( zmq_socket_pull, buffer_i, key_kind*n_integrals, 0)' + stop 'error' + endif + + rc = f77_zmq_recv( zmq_socket_pull, buffer_value, integral_kind*n_integrals, 0) + if (rc /= integral_kind*n_integrals) then + print *, irp_here, ': f77_zmq_recv( zmq_socket_pull, buffer_value, integral_kind*n_integrals, 0)' + stop 'error' + endif + + rc = f77_zmq_recv( zmq_socket_pull, task_id, 4, 0) + +! Activate if zmq_socket_pull is a REP + rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0) + if (rc /= 4) then + print *, irp_here, ' : f77_zmq_send (zmq_socket_pull,...' + stop 'error' + endif + + + call insert_into_ao_integrals_erf_map(n_integrals,buffer_i,buffer_value) + accu += n_integrals + if (task_id /= 0) then + call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,more) + endif + endif + + enddo + + deallocate( buffer_i, buffer_value ) + + integer (map_size_kind) :: get_ao_erf_map_size + control = get_ao_erf_map_size(ao_integrals_erf_map) + + if (control /= accu) then + print *, '' + print *, irp_here + print *, 'Control : ', control + print *, 'Accu : ', accu + print *, 'Some integrals were lost during the parallel computation.' + print *, 'Try to reduce the number of threads.' + stop + endif + + call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) + call end_zmq_pull_socket(zmq_socket_pull) + +end + diff --git a/plugins/Integrals_erf/ao_bielec_integrals_in_map_slave.irp.f b/plugins/Integrals_erf/ao_bielec_integrals_in_map_slave.irp.f new file mode 100644 index 00000000..38c78388 --- /dev/null +++ b/plugins/Integrals_erf/ao_bielec_integrals_in_map_slave.irp.f @@ -0,0 +1,225 @@ +subroutine ao_bielec_integrals_in_map_slave_tcp(i) + implicit none + integer, intent(in) :: i + BEGIN_DOC +! Computes a buffer of integrals. i is the ID of the current thread. + END_DOC + call ao_bielec_integrals_in_map_slave(0,i) +end + + +subroutine ao_bielec_integrals_in_map_slave_inproc(i) + implicit none + integer, intent(in) :: i + BEGIN_DOC +! Computes a buffer of integrals. i is the ID of the current thread. + END_DOC + call ao_bielec_integrals_in_map_slave(1,i) +end + + +subroutine push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value, task_id) + use f77_zmq + use map_module + implicit none + BEGIN_DOC +! Push integrals in the push socket + END_DOC + integer(ZMQ_PTR), intent(in) :: zmq_socket_push + integer, intent(in) :: n_integrals + integer(key_kind), intent(in) :: buffer_i(*) + real(integral_kind), intent(in) :: buffer_value(*) + integer, intent(in) :: task_id + integer :: rc + + rc = f77_zmq_send( zmq_socket_push, n_integrals, 4, ZMQ_SNDMORE) + if (rc /= 4) then + print *, irp_here, ': f77_zmq_send( zmq_socket_push, n_integrals, 4, ZMQ_SNDMORE)' + stop 'error' + endif + + rc = f77_zmq_send( zmq_socket_push, buffer_i, key_kind*n_integrals, ZMQ_SNDMORE) + if (rc /= key_kind*n_integrals) then + print *, irp_here, ': f77_zmq_send( zmq_socket_push, buffer_i, key_kind*n_integrals, ZMQ_SNDMORE)' + stop 'error' + endif + + rc = f77_zmq_send( zmq_socket_push, buffer_value, integral_kind*n_integrals, ZMQ_SNDMORE) + if (rc /= integral_kind*n_integrals) then + print *, irp_here, ': f77_zmq_send( zmq_socket_push, buffer_value, integral_kind*n_integrals, 0)' + stop 'error' + endif + + rc = f77_zmq_send( zmq_socket_push, task_id, 4, 0) + if (rc /= 4) then + print *, irp_here, ': f77_zmq_send( zmq_socket_push, task_id, 4, 0)' + stop 'error' + endif + +! Activate is zmq_socket_push is a REQ + integer :: idummy + rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0) + if (rc /= 4) then + print *, irp_here, ': f77_zmq_send( zmq_socket_push, idummy, 4, 0)' + stop 'error' + endif +end + + + + + +subroutine ao_bielec_integrals_in_map_slave(thread,iproc) + use map_module + use f77_zmq + implicit none + BEGIN_DOC +! Computes a buffer of integrals + END_DOC + + integer, intent(in) :: thread, iproc + + integer :: j,l,n_integrals + integer :: rc + real(integral_kind), allocatable :: buffer_value(:) + integer(key_kind), allocatable :: buffer_i(:) + + integer :: worker_id, task_id + character*(512) :: task + + integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + + integer(ZMQ_PTR), external :: new_zmq_push_socket + integer(ZMQ_PTR) :: zmq_socket_push + + character*(64) :: state + + zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() + zmq_socket_push = new_zmq_push_socket(thread) + + allocate ( buffer_i(ao_num*ao_num), buffer_value(ao_num*ao_num) ) + + call connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) + + do + call get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task) + if (task_id == 0) exit + read(task,*) j, l + call compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value) + call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id) + call push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value, task_id) + enddo + + + call disconnect_from_taskserver(zmq_to_qp_run_socket,zmq_socket_push,worker_id) + deallocate( buffer_i, buffer_value ) + call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) + call end_zmq_push_socket(zmq_socket_push,thread) + +end + + +subroutine ao_bielec_integrals_in_map_collector + use map_module + use f77_zmq + implicit none + BEGIN_DOC +! Collects results from the AO integral calculation + END_DOC + + integer :: j,l,n_integrals + integer :: rc + + real(integral_kind), allocatable :: buffer_value(:) + integer(key_kind), allocatable :: buffer_i(:) + + integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + + integer(ZMQ_PTR), external :: new_zmq_pull_socket + integer(ZMQ_PTR) :: zmq_socket_pull + + integer*8 :: control, accu + integer :: task_id, more, sze + + zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() + zmq_socket_pull = new_zmq_pull_socket() + + sze = ao_num*ao_num + allocate ( buffer_i(sze), buffer_value(sze) ) + + accu = 0_8 + more = 1 + do while (more == 1) + + rc = f77_zmq_recv( zmq_socket_pull, n_integrals, 4, 0) + if (rc == -1) then + n_integrals = 0 + return + endif + if (rc /= 4) then + print *, irp_here, ': f77_zmq_recv( zmq_socket_pull, n_integrals, 4, 0)' + stop 'error' + endif + + if (n_integrals >= 0) then + + if (n_integrals > sze) then + deallocate (buffer_value, buffer_i) + sze = n_integrals + allocate (buffer_value(sze), buffer_i(sze)) + endif + + rc = f77_zmq_recv( zmq_socket_pull, buffer_i, key_kind*n_integrals, 0) + if (rc /= key_kind*n_integrals) then + print *, rc, key_kind, n_integrals + print *, irp_here, ': f77_zmq_recv( zmq_socket_pull, buffer_i, key_kind*n_integrals, 0)' + stop 'error' + endif + + rc = f77_zmq_recv( zmq_socket_pull, buffer_value, integral_kind*n_integrals, 0) + if (rc /= integral_kind*n_integrals) then + print *, irp_here, ': f77_zmq_recv( zmq_socket_pull, buffer_value, integral_kind*n_integrals, 0)' + stop 'error' + endif + + rc = f77_zmq_recv( zmq_socket_pull, task_id, 4, 0) + +! Activate if zmq_socket_pull is a REP + rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0) + if (rc /= 4) then + print *, irp_here, ' : f77_zmq_send (zmq_socket_pull,...' + stop 'error' + endif + + + call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value) + accu += n_integrals + if (task_id /= 0) then + call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,more) + endif + endif + + enddo + + deallocate( buffer_i, buffer_value ) + + integer (map_size_kind) :: get_ao_map_size + control = get_ao_map_size(ao_integrals_map) + + if (control /= accu) then + print *, '' + print *, irp_here + print *, 'Control : ', control + print *, 'Accu : ', accu + print *, 'Some integrals were lost during the parallel computation.' + print *, 'Try to reduce the number of threads.' + stop + endif + + call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) + call end_zmq_pull_socket(zmq_socket_pull) + +end + diff --git a/plugins/Integrals_erf/gauss_legendre.irp.f b/plugins/Integrals_erf/gauss_legendre.irp.f new file mode 100644 index 00000000..e62febeb --- /dev/null +++ b/plugins/Integrals_erf/gauss_legendre.irp.f @@ -0,0 +1,66 @@ +BEGIN_PROVIDER [ integer, n_pt_max_integrals_16 ] + implicit none + BEGIN_DOC + ! Aligned n_pt_max_integrals + END_DOC + integer, external :: align_double + n_pt_max_integrals_16 = align_double(n_pt_max_integrals) +END_PROVIDER + + BEGIN_PROVIDER [ double precision, gauleg_t2, (n_pt_max_integrals_16,n_pt_max_integrals/2) ] +&BEGIN_PROVIDER [ double precision, gauleg_w, (n_pt_max_integrals_16,n_pt_max_integrals/2) ] + implicit none + BEGIN_DOC + ! t_w(i,1,k) = w(i) + ! t_w(i,2,k) = t(i) + END_DOC + integer :: i,j,l + l=0 + do i = 2,n_pt_max_integrals,2 + l = l+1 + call gauleg(0.d0,1.d0,gauleg_t2(1,l),gauleg_w(1,l),i) + do j=1,i + gauleg_t2(j,l) *= gauleg_t2(j,l) + enddo + enddo + +END_PROVIDER + +subroutine gauleg(x1,x2,x,w,n) + implicit none + BEGIN_DOC + ! Gauss-Legendre + END_DOC + integer, intent(in) :: n + double precision, intent(in) :: x1, x2 + double precision, intent (out) :: x(n),w(n) + double precision, parameter :: eps=3.d-14 + + integer :: m,i,j + double precision :: xm, xl, z, z1, p1, p2, p3, pp, dn + m=(n+1)/2 + xm=0.5d0*(x2+x1) + xl=0.5d0*(x2-x1) + dn = dble(n) + do i=1,m + z=dcos(3.141592654d0*(dble(i)-.25d0)/(dble(n)+.5d0)) + z1 = z+1.d0 + do while (dabs(z-z1) > eps) + p1=1.d0 + p2=0.d0 + do j=1,n + p3=p2 + p2=p1 + p1=(dble(j+j-1)*z*p2-dble(j-1)*p3)/j + enddo + pp=dn*(z*p1-p2)/(z*z-1.d0) + z1=z + z=z1-p1/pp + end do + x(i)=xm-xl*z + x(n+1-i)=xm+xl*z + w(i)=(xl+xl)/((1.d0-z*z)*pp*pp) + w(n+1-i)=w(i) + enddo +end + diff --git a/plugins/Integrals_erf/integrals_3_index.irp.f b/plugins/Integrals_erf/integrals_3_index.irp.f new file mode 100644 index 00000000..41037b34 --- /dev/null +++ b/plugins/Integrals_erf/integrals_3_index.irp.f @@ -0,0 +1,22 @@ + BEGIN_PROVIDER [double precision, big_array_coulomb_integrals, (mo_tot_num_align,mo_tot_num, mo_tot_num)] +&BEGIN_PROVIDER [double precision, big_array_exchange_integrals,(mo_tot_num_align,mo_tot_num, mo_tot_num)] + implicit none + integer :: i,j,k,l + double precision :: get_mo_bielec_integral + double precision :: integral + + do k = 1, mo_tot_num + do i = 1, mo_tot_num + do j = 1, mo_tot_num + l = j + integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + big_array_coulomb_integrals(j,i,k) = integral + l = j + integral = get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + big_array_exchange_integrals(j,i,k) = integral + enddo + enddo + enddo + + +END_PROVIDER diff --git a/plugins/Integrals_erf/map_integrals.irp.f b/plugins/Integrals_erf/map_integrals.irp.f new file mode 100644 index 00000000..82b89f22 --- /dev/null +++ b/plugins/Integrals_erf/map_integrals.irp.f @@ -0,0 +1,717 @@ +use map_module + +!! AO Map +!! ====== + +BEGIN_PROVIDER [ type(map_type), ao_integrals_map ] + implicit none + BEGIN_DOC + ! AO integrals + END_DOC + integer(key_kind) :: key_max + integer(map_size_kind) :: sze + call bielec_integrals_index(ao_num,ao_num,ao_num,ao_num,key_max) + sze = key_max + call map_init(ao_integrals_map,sze) + print*, 'AO map initialized : ', sze +END_PROVIDER + +subroutine bielec_integrals_index(i,j,k,l,i1) + use map_module + implicit none + integer, intent(in) :: i,j,k,l + integer(key_kind), intent(out) :: i1 + integer(key_kind) :: p,q,r,s,i2 + p = min(i,k) + r = max(i,k) + p = p+ishft(r*r-r,-1) + q = min(j,l) + s = max(j,l) + q = q+ishft(s*s-s,-1) + i1 = min(p,q) + i2 = max(p,q) + i1 = i1+ishft(i2*i2-i2,-1) +end + +subroutine bielec_integrals_index_reverse(i,j,k,l,i1) + use map_module + implicit none + integer, intent(out) :: i(8),j(8),k(8),l(8) + integer(key_kind), intent(in) :: i1 + integer(key_kind) :: i2,i3 + i = 0 + i2 = ceiling(0.5d0*(dsqrt(8.d0*dble(i1)+1.d0)-1.d0)) + l(1) = ceiling(0.5d0*(dsqrt(8.d0*dble(i2)+1.d0)-1.d0)) + i3 = i1 - ishft(i2*i2-i2,-1) + k(1) = ceiling(0.5d0*(dsqrt(8.d0*dble(i3)+1.d0)-1.d0)) + j(1) = int(i2 - ishft(l(1)*l(1)-l(1),-1),4) + i(1) = int(i3 - ishft(k(1)*k(1)-k(1),-1),4) + + !ijkl + i(2) = i(1) !ilkj + j(2) = l(1) + k(2) = k(1) + l(2) = j(1) + + i(3) = k(1) !kjil + j(3) = j(1) + k(3) = i(1) + l(3) = l(1) + + i(4) = k(1) !klij + j(4) = l(1) + k(4) = i(1) + l(4) = j(1) + + i(5) = j(1) !jilk + j(5) = i(1) + k(5) = l(1) + l(5) = k(1) + + i(6) = j(1) !jkli + j(6) = k(1) + k(6) = l(1) + l(6) = i(1) + + i(7) = l(1) !lijk + j(7) = i(1) + k(7) = j(1) + l(7) = k(1) + + i(8) = l(1) !lkji + j(8) = k(1) + k(8) = j(1) + l(8) = i(1) + + integer :: ii, jj + do ii=2,8 + do jj=1,ii-1 + if ( (i(ii) == i(jj)).and. & + (j(ii) == j(jj)).and. & + (k(ii) == k(jj)).and. & + (l(ii) == l(jj)) ) then + i(ii) = 0 + exit + endif + enddo + enddo + do ii=1,8 + if (i(ii) /= 0) then + call bielec_integrals_index(i(ii),j(ii),k(ii),l(ii),i2) + if (i1 /= i2) then + print *, i1, i2 + print *, i(ii), j(jj), k(jj), l(jj) + stop 'bielec_integrals_index_reverse failed' + endif + endif + enddo + + +end + + BEGIN_PROVIDER [ integer, ao_integrals_cache_min ] +&BEGIN_PROVIDER [ integer, ao_integrals_cache_max ] + implicit none + BEGIN_DOC + ! Min and max values of the AOs for which the integrals are in the cache + END_DOC + ao_integrals_cache_min = max(1,ao_num - 63) + ao_integrals_cache_max = ao_num + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_integrals_cache, (0:64*64*64*64) ] + implicit none + BEGIN_DOC + ! Cache of AO integrals for fast access + END_DOC + PROVIDE ao_bielec_integrals_in_map + integer :: i,j,k,l,ii + integer(key_kind) :: idx + real(integral_kind) :: integral + !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral) + do l=ao_integrals_cache_min,ao_integrals_cache_max + do k=ao_integrals_cache_min,ao_integrals_cache_max + do j=ao_integrals_cache_min,ao_integrals_cache_max + do i=ao_integrals_cache_min,ao_integrals_cache_max + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,idx) + !DIR$ FORCEINLINE + call map_get(ao_integrals_map,idx,integral) + ii = l-ao_integrals_cache_min + ii = ior( ishft(ii,6), k-ao_integrals_cache_min) + ii = ior( ishft(ii,6), j-ao_integrals_cache_min) + ii = ior( ishft(ii,6), i-ao_integrals_cache_min) + ao_integrals_cache(ii) = integral + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + + +double precision function get_ao_bielec_integral(i,j,k,l,map) result(result) + use map_module + implicit none + BEGIN_DOC + ! Gets one AO bi-electronic integral from the AO map + END_DOC + integer, intent(in) :: i,j,k,l + integer(key_kind) :: idx + type(map_type), intent(inout) :: map + integer :: ii + real(integral_kind) :: tmp + PROVIDE ao_bielec_integrals_in_map ao_integrals_cache ao_integrals_cache_min + !DIR$ FORCEINLINE + if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then + tmp = 0.d0 + else if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < ao_integrals_threshold) then + tmp = 0.d0 + else + ii = l-ao_integrals_cache_min + ii = ior(ii, k-ao_integrals_cache_min) + ii = ior(ii, j-ao_integrals_cache_min) + ii = ior(ii, i-ao_integrals_cache_min) + if (iand(ii, -64) /= 0) then + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,idx) + !DIR$ FORCEINLINE + call map_get(map,idx,tmp) + tmp = tmp + else + ii = l-ao_integrals_cache_min + ii = ior( ishft(ii,6), k-ao_integrals_cache_min) + ii = ior( ishft(ii,6), j-ao_integrals_cache_min) + ii = ior( ishft(ii,6), i-ao_integrals_cache_min) + tmp = ao_integrals_cache(ii) + endif + endif + result = tmp +end + + +subroutine get_ao_bielec_integrals(j,k,l,sze,out_val) + use map_module + BEGIN_DOC + ! Gets multiple AO bi-electronic integral from the AO map . + ! All i are retrieved for j,k,l fixed. + END_DOC + implicit none + integer, intent(in) :: j,k,l, sze + real(integral_kind), intent(out) :: out_val(sze) + + integer :: i + integer(key_kind) :: hash + double precision :: thresh + PROVIDE ao_bielec_integrals_in_map ao_integrals_map + thresh = ao_integrals_threshold + + if (ao_overlap_abs(j,l) < thresh) then + out_val = 0.d0 + return + endif + + double precision :: get_ao_bielec_integral + do i=1,sze + out_val(i) = get_ao_bielec_integral(i,j,k,l,ao_integrals_map) + enddo + +end + +subroutine get_ao_bielec_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_zero_int) + use map_module + implicit none + BEGIN_DOC + ! Gets multiple AO bi-electronic integral from the AO map . + ! All non-zero i are retrieved for j,k,l fixed. + END_DOC + integer, intent(in) :: j,k,l, sze + real(integral_kind), intent(out) :: out_val(sze) + integer, intent(out) :: out_val_index(sze),non_zero_int + + integer :: i + integer(key_kind) :: hash + double precision :: thresh,tmp + PROVIDE ao_bielec_integrals_in_map + thresh = ao_integrals_threshold + + non_zero_int = 0 + if (ao_overlap_abs(j,l) < thresh) then + out_val = 0.d0 + return + endif + + non_zero_int = 0 + do i=1,sze + integer, external :: ao_l4 + double precision, external :: ao_bielec_integral + !DIR$ FORCEINLINE + if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh) then + cycle + endif + call bielec_integrals_index(i,j,k,l,hash) + call map_get(ao_integrals_map, hash,tmp) + if (dabs(tmp) < thresh ) cycle + non_zero_int = non_zero_int+1 + out_val_index(non_zero_int) = i + out_val(non_zero_int) = tmp + enddo + +end + + +function get_ao_map_size() + implicit none + integer (map_size_kind) :: get_ao_map_size + BEGIN_DOC + ! Returns the number of elements in the AO map + END_DOC + get_ao_map_size = ao_integrals_map % n_elements +end + +subroutine clear_ao_map + implicit none + BEGIN_DOC + ! Frees the memory of the AO map + END_DOC + call map_deinit(ao_integrals_map) + FREE ao_integrals_map +end + + +!! MO Map +!! ====== + +BEGIN_PROVIDER [ type(map_type), mo_integrals_map ] + implicit none + BEGIN_DOC + ! MO integrals + END_DOC + integer(key_kind) :: key_max + integer(map_size_kind) :: sze + call bielec_integrals_index(mo_tot_num,mo_tot_num,mo_tot_num,mo_tot_num,key_max) + sze = key_max + call map_init(mo_integrals_map,sze) + print*, 'MO map initialized' +END_PROVIDER + +subroutine insert_into_ao_integrals_map(n_integrals,buffer_i, buffer_values) + use map_module + implicit none + BEGIN_DOC + ! Create new entry into AO map + END_DOC + + integer, intent(in) :: n_integrals + integer(key_kind), intent(inout) :: buffer_i(n_integrals) + real(integral_kind), intent(inout) :: buffer_values(n_integrals) + + call map_append(ao_integrals_map, buffer_i, buffer_values, n_integrals) +end + +subroutine insert_into_mo_integrals_map(n_integrals, & + buffer_i, buffer_values, thr) + use map_module + implicit none + + BEGIN_DOC + ! Create new entry into MO map, or accumulate in an existing entry + END_DOC + + integer, intent(in) :: n_integrals + integer(key_kind), intent(inout) :: buffer_i(n_integrals) + real(integral_kind), intent(inout) :: buffer_values(n_integrals) + real(integral_kind), intent(in) :: thr + call map_update(mo_integrals_map, buffer_i, buffer_values, n_integrals, thr) +end + + BEGIN_PROVIDER [ integer, mo_integrals_cache_min ] +&BEGIN_PROVIDER [ integer, mo_integrals_cache_max ] + implicit none + BEGIN_DOC + ! Min and max values of the MOs for which the integrals are in the cache + END_DOC + mo_integrals_cache_min = max(1,elec_alpha_num - 31) + mo_integrals_cache_max = min(mo_tot_num,mo_integrals_cache_min+63) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0:64*64*64*64) ] + implicit none + BEGIN_DOC + ! Cache of MO integrals for fast access + END_DOC + PROVIDE mo_bielec_integrals_in_map + integer :: i,j,k,l + integer :: ii + integer(key_kind) :: idx + real(integral_kind) :: integral + FREE ao_integrals_cache + !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral) + do l=mo_integrals_cache_min,mo_integrals_cache_max + do k=mo_integrals_cache_min,mo_integrals_cache_max + do j=mo_integrals_cache_min,mo_integrals_cache_max + do i=mo_integrals_cache_min,mo_integrals_cache_max + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,idx) + !DIR$ FORCEINLINE + call map_get(mo_integrals_map,idx,integral) + ii = l-mo_integrals_cache_min + ii = ior( ishft(ii,6), k-mo_integrals_cache_min) + ii = ior( ishft(ii,6), j-mo_integrals_cache_min) + ii = ior( ishft(ii,6), i-mo_integrals_cache_min) + mo_integrals_cache(ii) = integral + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + + +double precision function get_mo_bielec_integral(i,j,k,l,map) + use map_module + implicit none + BEGIN_DOC + ! Returns one integral in the MO basis + END_DOC + integer, intent(in) :: i,j,k,l + integer(key_kind) :: idx + integer :: ii + type(map_type), intent(inout) :: map + real(integral_kind) :: tmp + PROVIDE mo_bielec_integrals_in_map mo_integrals_cache + ii = l-mo_integrals_cache_min + ii = ior(ii, k-mo_integrals_cache_min) + ii = ior(ii, j-mo_integrals_cache_min) + ii = ior(ii, i-mo_integrals_cache_min) + if (iand(ii, -64) /= 0) then + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,idx) + !DIR$ FORCEINLINE + call map_get(map,idx,tmp) + get_mo_bielec_integral = dble(tmp) + else + ii = l-mo_integrals_cache_min + ii = ior( ishft(ii,6), k-mo_integrals_cache_min) + ii = ior( ishft(ii,6), j-mo_integrals_cache_min) + ii = ior( ishft(ii,6), i-mo_integrals_cache_min) + get_mo_bielec_integral = mo_integrals_cache(ii) + endif +end + + +double precision function mo_bielec_integral(i,j,k,l) + implicit none + BEGIN_DOC + ! Returns one integral in the MO basis + END_DOC + integer, intent(in) :: i,j,k,l + double precision :: get_mo_bielec_integral + PROVIDE mo_bielec_integrals_in_map mo_integrals_cache + !DIR$ FORCEINLINE + PROVIDE mo_bielec_integrals_in_map + mo_bielec_integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + return +end + +subroutine get_mo_bielec_integrals(j,k,l,sze,out_val,map) + use map_module + implicit none + BEGIN_DOC + ! Returns multiple integrals in the MO basis, all + ! i for j,k,l fixed. + END_DOC + integer, intent(in) :: j,k,l, sze + double precision, intent(out) :: out_val(sze) + type(map_type), intent(inout) :: map + integer :: i + integer(key_kind) :: hash(sze) + real(integral_kind) :: tmp_val(sze) + PROVIDE mo_bielec_integrals_in_map + + do i=1,sze + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,hash(i)) + enddo + + if (key_kind == 8) then + call map_get_many(map, hash, out_val, sze) + else + call map_get_many(map, hash, tmp_val, sze) + ! Conversion to double precision + do i=1,sze + out_val(i) = dble(tmp_val(i)) + enddo + endif +end + +subroutine get_mo_bielec_integrals_ij(k,l,sze,out_array,map) + use map_module + implicit none + BEGIN_DOC + ! Returns multiple integrals in the MO basis, all + ! i(1)j(2) 1/r12 k(1)l(2) + ! i, j for k,l fixed. + END_DOC + integer, intent(in) :: k,l, sze + double precision, intent(out) :: out_array(sze,sze) + type(map_type), intent(inout) :: map + integer :: i,j,kk,ll,m + integer(key_kind),allocatable :: hash(:) + integer ,allocatable :: pairs(:,:), iorder(:) + real(integral_kind), allocatable :: tmp_val(:) + + PROVIDE mo_bielec_integrals_in_map + allocate (hash(sze*sze), pairs(2,sze*sze),iorder(sze*sze), & + tmp_val(sze*sze)) + + kk=0 + out_array = 0.d0 + do j=1,sze + do i=1,sze + kk += 1 + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,hash(kk)) + pairs(1,kk) = i + pairs(2,kk) = j + iorder(kk) = kk + enddo + enddo + + logical :: integral_is_in_map + if (key_kind == 8) then + call i8radix_sort(hash,iorder,kk,-1) + else if (key_kind == 4) then + call iradix_sort(hash,iorder,kk,-1) + else if (key_kind == 2) then + call i2radix_sort(hash,iorder,kk,-1) + endif + + call map_get_many(mo_integrals_map, hash, tmp_val, kk) + + do ll=1,kk + m = iorder(ll) + i=pairs(1,m) + j=pairs(2,m) + out_array(i,j) = tmp_val(ll) + enddo + + deallocate(pairs,hash,iorder,tmp_val) +end + +subroutine get_mo_bielec_integrals_coulomb_ii(k,l,sze,out_val,map) + use map_module + implicit none + BEGIN_DOC + ! Returns multiple integrals + ! k(1)i(2) 1/r12 l(1)i(2) :: out_val(i1) + ! for k,l fixed. + END_DOC + integer, intent(in) :: k,l, sze + double precision, intent(out) :: out_val(sze) + type(map_type), intent(inout) :: map + integer :: i + integer(key_kind) :: hash(sze) + real(integral_kind) :: tmp_val(sze) + PROVIDE mo_bielec_integrals_in_map + + integer :: kk + do i=1,sze + !DIR$ FORCEINLINE + call bielec_integrals_index(k,i,l,i,hash(i)) + enddo + + if (key_kind == 8) then + call map_get_many(map, hash, out_val, sze) + else + call map_get_many(map, hash, tmp_val, sze) + ! Conversion to double precision + do i=1,sze + out_val(i) = dble(tmp_val(i)) + enddo + endif +end + +subroutine get_mo_bielec_integrals_exch_ii(k,l,sze,out_val,map) + use map_module + implicit none + BEGIN_DOC + ! Returns multiple integrals + ! k(1)i(2) 1/r12 i(1)l(2) :: out_val(i1) + ! for k,l fixed. + END_DOC + integer, intent(in) :: k,l, sze + double precision, intent(out) :: out_val(sze) + type(map_type), intent(inout) :: map + integer :: i + integer(key_kind) :: hash(sze) + real(integral_kind) :: tmp_val(sze) + PROVIDE mo_bielec_integrals_in_map + + integer :: kk + do i=1,sze + !DIR$ FORCEINLINE + call bielec_integrals_index(k,i,i,l,hash(i)) + enddo + + if (key_kind == 8) then + call map_get_many(map, hash, out_val, sze) + else + call map_get_many(map, hash, tmp_val, sze) + ! Conversion to double precision + do i=1,sze + out_val(i) = dble(tmp_val(i)) + enddo + endif +end + + +integer*8 function get_mo_map_size() + implicit none + BEGIN_DOC + ! Return the number of elements in the MO map + END_DOC + get_mo_map_size = mo_integrals_map % n_elements +end + +BEGIN_TEMPLATE + +subroutine dump_$ao_integrals(filename) + use map_module + implicit none + BEGIN_DOC + ! Save to disk the $ao integrals + END_DOC + character*(*), intent(in) :: filename + integer(cache_key_kind), pointer :: key(:) + real(integral_kind), pointer :: val(:) + integer*8 :: i,j, n + call ezfio_set_work_empty(.False.) + open(unit=66,file=filename,FORM='unformatted') + write(66) integral_kind, key_kind + write(66) $ao_integrals_map%sorted, $ao_integrals_map%map_size, & + $ao_integrals_map%n_elements + do i=0_8,$ao_integrals_map%map_size + write(66) $ao_integrals_map%map(i)%sorted, $ao_integrals_map%map(i)%map_size,& + $ao_integrals_map%map(i)%n_elements + enddo + do i=0_8,$ao_integrals_map%map_size + key => $ao_integrals_map%map(i)%key + val => $ao_integrals_map%map(i)%value + n = $ao_integrals_map%map(i)%n_elements + write(66) (key(j), j=1,n), (val(j), j=1,n) + enddo + close(66) + +end + +IRP_IF COARRAY +subroutine communicate_$ao_integrals() + use map_module + implicit none + BEGIN_DOC + ! Communicate the $ao integrals with co-array + END_DOC + integer(cache_key_kind), pointer :: key(:) + real(integral_kind), pointer :: val(:) + integer*8 :: i,j, k, nmax + integer*8, save :: n[*] + integer :: copy_n + + real(integral_kind), allocatable :: buffer_val(:)[:] + integer(cache_key_kind), allocatable :: buffer_key(:)[:] + real(integral_kind), allocatable :: copy_val(:) + integer(key_kind), allocatable :: copy_key(:) + + n = 0_8 + do i=0_8,$ao_integrals_map%map_size + n = max(n,$ao_integrals_map%map(i)%n_elements) + enddo + sync all + nmax = 0_8 + do j=1,num_images() + nmax = max(nmax,n[j]) + enddo + allocate( buffer_key(nmax)[*], buffer_val(nmax)[*]) + allocate( copy_key(nmax), copy_val(nmax)) + do i=0_8,$ao_integrals_map%map_size + key => $ao_integrals_map%map(i)%key + val => $ao_integrals_map%map(i)%value + n = $ao_integrals_map%map(i)%n_elements + do j=1,n + buffer_key(j) = key(j) + buffer_val(j) = val(j) + enddo + sync all + do j=1,num_images() + if (j /= this_image()) then + copy_n = n[j] + do k=1,copy_n + copy_val(k) = buffer_val(k)[j] + copy_key(k) = buffer_key(k)[j] + copy_key(k) = copy_key(k)+ishft(i,-map_shift) + enddo + call map_append($ao_integrals_map, copy_key, copy_val, copy_n ) + endif + enddo + sync all + enddo + deallocate( buffer_key, buffer_val, copy_val, copy_key) + +end +IRP_ENDIF + + +integer function load_$ao_integrals(filename) + implicit none + BEGIN_DOC + ! Read from disk the $ao integrals + END_DOC + character*(*), intent(in) :: filename + integer*8 :: i + integer(cache_key_kind), pointer :: key(:) + real(integral_kind), pointer :: val(:) + integer :: iknd, kknd + integer*8 :: n, j + load_$ao_integrals = 1 + open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN') + read(66,err=98,end=98) iknd, kknd + if (iknd /= integral_kind) then + print *, 'Wrong integrals kind in file :', iknd + stop 1 + endif + if (kknd /= key_kind) then + print *, 'Wrong key kind in file :', kknd + stop 1 + endif + read(66,err=98,end=98) $ao_integrals_map%sorted, $ao_integrals_map%map_size,& + $ao_integrals_map%n_elements + do i=0_8, $ao_integrals_map%map_size + read(66,err=99,end=99) $ao_integrals_map%map(i)%sorted, & + $ao_integrals_map%map(i)%map_size, $ao_integrals_map%map(i)%n_elements + call cache_map_reallocate($ao_integrals_map%map(i),$ao_integrals_map%map(i)%map_size) + enddo + do i=0_8, $ao_integrals_map%map_size + key => $ao_integrals_map%map(i)%key + val => $ao_integrals_map%map(i)%value + n = $ao_integrals_map%map(i)%n_elements + read(66,err=99,end=99) (key(j), j=1,n), (val(j), j=1,n) + enddo + call map_sort($ao_integrals_map) + load_$ao_integrals = 0 + return + 99 continue + call map_deinit($ao_integrals_map) + 98 continue + stop 'Problem reading $ao_integrals_map file in work/' + +end + +SUBST [ ao_integrals_map, ao_integrals, ao_num ] +ao_integrals_map ; ao_integrals ; ao_num ;; +mo_integrals_map ; mo_integrals ; mo_tot_num ;; +END_TEMPLATE diff --git a/plugins/Integrals_erf/map_integrals_long_range.irp.f b/plugins/Integrals_erf/map_integrals_long_range.irp.f new file mode 100644 index 00000000..2ab2f14f --- /dev/null +++ b/plugins/Integrals_erf/map_integrals_long_range.irp.f @@ -0,0 +1,611 @@ +use map_module + +!! AO Map +!! ====== + +BEGIN_PROVIDER [ type(map_type), ao_integrals_erf_map ] + implicit none + BEGIN_DOC + ! AO integrals + END_DOC + integer(key_kind) :: key_max + integer(map_size_kind) :: sze + call bielec_integrals_index(ao_num,ao_num,ao_num,ao_num,key_max) + sze = key_max + call map_init(ao_integrals_erf_map,sze) + print*, 'AO map initialized : ', sze +END_PROVIDER + + BEGIN_PROVIDER [ integer, ao_integrals_erf_cache_min ] +&BEGIN_PROVIDER [ integer, ao_integrals_erf_cache_max ] + implicit none + BEGIN_DOC + ! Min and max values of the AOs for which the integrals are in the cache + END_DOC + ao_integrals_erf_cache_min = max(1,ao_num - 63) + ao_integrals_erf_cache_max = ao_num + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_integrals_erf_cache, (0:64*64*64*64) ] + implicit none + BEGIN_DOC + ! Cache of AO integrals for fast access + END_DOC + PROVIDE ao_bielec_integrals_erf_in_map + integer :: i,j,k,l,ii + integer(key_kind) :: idx + real(integral_kind) :: integral + !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral) + do l=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max + do k=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max + do j=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max + do i=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max + !DIR$ FORCEINLINE + call bielec_integrals_erf_index(i,j,k,l,idx) + !DIR$ FORCEINLINE + call map_get(ao_integrals_erf_map,idx,integral) + ii = l-ao_integrals_erf_cache_min + ii = ior( ishft(ii,6), k-ao_integrals_erf_cache_min) + ii = ior( ishft(ii,6), j-ao_integrals_erf_cache_min) + ii = ior( ishft(ii,6), i-ao_integrals_erf_cache_min) + ao_integrals_erf_cache(ii) = integral + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + + +double precision function get_ao_bielec_integral_erf(i,j,k,l,map) result(result) + use map_module + implicit none + BEGIN_DOC + ! Gets one AO bi-electronic integral from the AO map + END_DOC + integer, intent(in) :: i,j,k,l + integer(key_kind) :: idx + type(map_type), intent(inout) :: map + integer :: ii + real(integral_kind) :: tmp + PROVIDE ao_bielec_integrals_erf_in_map ao_integrals_erf_cache ao_integrals_erf_cache_min + !DIR$ FORCEINLINE + if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then + tmp = 0.d0 + else if (ao_bielec_integral_erf_schwartz(i,k)*ao_bielec_integral_erf_schwartz(j,l) < ao_integrals_threshold) then + tmp = 0.d0 + else + ii = l-ao_integrals_erf_cache_min + ii = ior(ii, k-ao_integrals_erf_cache_min) + ii = ior(ii, j-ao_integrals_erf_cache_min) + ii = ior(ii, i-ao_integrals_erf_cache_min) + if (iand(ii, -64) /= 0) then + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,idx) + !DIR$ FORCEINLINE + call map_get(map,idx,tmp) + tmp = tmp + else + ii = l-ao_integrals_erf_cache_min + ii = ior( ishft(ii,6), k-ao_integrals_erf_cache_min) + ii = ior( ishft(ii,6), j-ao_integrals_erf_cache_min) + ii = ior( ishft(ii,6), i-ao_integrals_erf_cache_min) + tmp = ao_integrals_erf_cache(ii) + endif + endif + result = tmp +end + + +subroutine get_ao_bielec_integrals_erf(j,k,l,sze,out_val) + use map_module + BEGIN_DOC + ! Gets multiple AO bi-electronic integral from the AO map . + ! All i are retrieved for j,k,l fixed. + END_DOC + implicit none + integer, intent(in) :: j,k,l, sze + real(integral_kind), intent(out) :: out_val(sze) + + integer :: i + integer(key_kind) :: hash + double precision :: thresh + PROVIDE ao_bielec_integrals_erf_in_map ao_integrals_erf_map + thresh = ao_integrals_threshold + + if (ao_overlap_abs(j,l) < thresh) then + out_val = 0.d0 + return + endif + + double precision :: get_ao_bielec_integral_erf + do i=1,sze + out_val(i) = get_ao_bielec_integral_erf(i,j,k,l,ao_integrals_erf_map) + enddo + +end + +subroutine get_ao_bielec_integrals_erf_non_zero(j,k,l,sze,out_val,out_val_index,non_zero_int) + use map_module + implicit none + BEGIN_DOC + ! Gets multiple AO bi-electronic integral from the AO map . + ! All non-zero i are retrieved for j,k,l fixed. + END_DOC + integer, intent(in) :: j,k,l, sze + real(integral_kind), intent(out) :: out_val(sze) + integer, intent(out) :: out_val_index(sze),non_zero_int + + integer :: i + integer(key_kind) :: hash + double precision :: thresh,tmp + PROVIDE ao_bielec_integrals_erf_in_map + thresh = ao_integrals_threshold + + non_zero_int = 0 + if (ao_overlap_abs(j,l) < thresh) then + out_val = 0.d0 + return + endif + + non_zero_int = 0 + do i=1,sze + integer, external :: ao_l4 + double precision, external :: ao_bielec_integral_erf + !DIR$ FORCEINLINE + if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh) then + cycle + endif + call bielec_integrals_index(i,j,k,l,hash) + call map_get(ao_integrals_erf_map, hash,tmp) + if (dabs(tmp) < thresh ) cycle + non_zero_int = non_zero_int+1 + out_val_index(non_zero_int) = i + out_val(non_zero_int) = tmp + enddo + +end + + +function get_ao_erf_map_size() + implicit none + integer (map_size_kind) :: get_ao_erf_map_size + BEGIN_DOC + ! Returns the number of elements in the AO map + END_DOC + get_ao_erf_map_size = ao_integrals_erf_map % n_elements +end + +subroutine clear_ao_erf_map + implicit none + BEGIN_DOC + ! Frees the memory of the AO map + END_DOC + call map_deinit(ao_integrals_erf_map) + FREE ao_integrals_erf_map +end + + + +BEGIN_TEMPLATE + +subroutine dump_$ao_integrals(filename) + use map_module + implicit none + BEGIN_DOC + ! Save to disk the $ao integrals + END_DOC + character*(*), intent(in) :: filename + integer(cache_key_kind), pointer :: key(:) + real(integral_kind), pointer :: val(:) + integer*8 :: i,j, n + call ezfio_set_work_empty(.False.) + open(unit=66,file=filename,FORM='unformatted') + write(66) integral_kind, key_kind + write(66) $ao_integrals_map%sorted, $ao_integrals_map%map_size, & + $ao_integrals_map%n_elements + do i=0_8,$ao_integrals_map%map_size + write(66) $ao_integrals_map%map(i)%sorted, $ao_integrals_map%map(i)%map_size,& + $ao_integrals_map%map(i)%n_elements + enddo + do i=0_8,$ao_integrals_map%map_size + key => $ao_integrals_map%map(i)%key + val => $ao_integrals_map%map(i)%value + n = $ao_integrals_map%map(i)%n_elements + write(66) (key(j), j=1,n), (val(j), j=1,n) + enddo + close(66) + +end + +IRP_IF COARRAY +subroutine communicate_$ao_integrals() + use map_module + implicit none + BEGIN_DOC + ! Communicate the $ao integrals with co-array + END_DOC + integer(cache_key_kind), pointer :: key(:) + real(integral_kind), pointer :: val(:) + integer*8 :: i,j, k, nmax + integer*8, save :: n[*] + integer :: copy_n + + real(integral_kind), allocatable :: buffer_val(:)[:] + integer(cache_key_kind), allocatable :: buffer_key(:)[:] + real(integral_kind), allocatable :: copy_val(:) + integer(key_kind), allocatable :: copy_key(:) + + n = 0_8 + do i=0_8,$ao_integrals_map%map_size + n = max(n,$ao_integrals_map%map(i)%n_elements) + enddo + sync all + nmax = 0_8 + do j=1,num_images() + nmax = max(nmax,n[j]) + enddo + allocate( buffer_key(nmax)[*], buffer_val(nmax)[*]) + allocate( copy_key(nmax), copy_val(nmax)) + do i=0_8,$ao_integrals_map%map_size + key => $ao_integrals_map%map(i)%key + val => $ao_integrals_map%map(i)%value + n = $ao_integrals_map%map(i)%n_elements + do j=1,n + buffer_key(j) = key(j) + buffer_val(j) = val(j) + enddo + sync all + do j=1,num_images() + if (j /= this_image()) then + copy_n = n[j] + do k=1,copy_n + copy_val(k) = buffer_val(k)[j] + copy_key(k) = buffer_key(k)[j] + copy_key(k) = copy_key(k)+ishft(i,-map_shift) + enddo + call map_append($ao_integrals_map, copy_key, copy_val, copy_n ) + endif + enddo + sync all + enddo + deallocate( buffer_key, buffer_val, copy_val, copy_key) + +end +IRP_ENDIF + + +integer function load_$ao_integrals(filename) + implicit none + BEGIN_DOC + ! Read from disk the $ao integrals + END_DOC + character*(*), intent(in) :: filename + integer*8 :: i + integer(cache_key_kind), pointer :: key(:) + real(integral_kind), pointer :: val(:) + integer :: iknd, kknd + integer*8 :: n, j + load_$ao_integrals = 1 + open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN') + read(66,err=98,end=98) iknd, kknd + if (iknd /= integral_kind) then + print *, 'Wrong integrals kind in file :', iknd + stop 1 + endif + if (kknd /= key_kind) then + print *, 'Wrong key kind in file :', kknd + stop 1 + endif + read(66,err=98,end=98) $ao_integrals_map%sorted, $ao_integrals_map%map_size,& + $ao_integrals_map%n_elements + do i=0_8, $ao_integrals_map%map_size + read(66,err=99,end=99) $ao_integrals_map%map(i)%sorted, & + $ao_integrals_map%map(i)%map_size, $ao_integrals_map%map(i)%n_elements + call cache_map_reallocate($ao_integrals_map%map(i),$ao_integrals_map%map(i)%map_size) + enddo + do i=0_8, $ao_integrals_map%map_size + key => $ao_integrals_map%map(i)%key + val => $ao_integrals_map%map(i)%value + n = $ao_integrals_map%map(i)%n_elements + read(66,err=99,end=99) (key(j), j=1,n), (val(j), j=1,n) + enddo + call map_sort($ao_integrals_map) + load_$ao_integrals = 0 + return + 99 continue + call map_deinit($ao_integrals_map) + 98 continue + stop 'Problem reading $ao_integrals_map file in work/' + +end + +SUBST [ ao_integrals_map, ao_integrals, ao_num ] +ao_integrals_erf_map ; ao_integrals_erf ; ao_num ;; +mo_integrals_erf_map ; mo_integrals_erf ; mo_tot_num;; +END_TEMPLATE + + + + +BEGIN_PROVIDER [ type(map_type), mo_integrals_erf_map ] + implicit none + BEGIN_DOC + ! MO integrals + END_DOC + integer(key_kind) :: key_max + integer(map_size_kind) :: sze + call bielec_integrals_index(mo_tot_num,mo_tot_num,mo_tot_num,mo_tot_num,key_max) + sze = key_max + call map_init(mo_integrals_erf_map,sze) + print*, 'MO map initialized' +END_PROVIDER + +subroutine insert_into_mo_integrals_map(n_integrals, & + buffer_i, buffer_values, thr) + use map_module + implicit none + + BEGIN_DOC + ! Create new entry into MO map, or accumulate in an existing entry + END_DOC + + integer, intent(in) :: n_integrals + integer(key_kind), intent(inout) :: buffer_i(n_integrals) + real(integral_kind), intent(inout) :: buffer_values(n_integrals) + real(integral_kind), intent(in) :: thr + call map_update(mo_integrals_erf_map, buffer_i, buffer_values, n_integrals, thr) +end + + BEGIN_PROVIDER [ integer, mo_integrals_erf_cache_min ] +&BEGIN_PROVIDER [ integer, mo_integrals_erf_cache_max ] + implicit none + BEGIN_DOC + ! Min and max values of the MOs for which the integrals are in the cache + END_DOC + mo_integrals_erf_cache_min = max(1,elec_alpha_num - 31) + mo_integrals_erf_cache_max = min(mo_tot_num,mo_integrals_erf_cache_min+63) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, mo_integrals_erf_cache, (0:64*64*64*64) ] + implicit none + BEGIN_DOC + ! Cache of MO integrals for fast access + END_DOC + PROVIDE mo_bielec_integrals_erf_in_map + integer :: i,j,k,l + integer :: ii + integer(key_kind) :: idx + real(integral_kind) :: integral + FREE ao_integrals_cache + !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral) + do l=mo_integrals_erf_cache_min,mo_integrals_erf_cache_max + do k=mo_integrals_erf_cache_min,mo_integrals_erf_cache_max + do j=mo_integrals_erf_cache_min,mo_integrals_erf_cache_max + do i=mo_integrals_erf_cache_min,mo_integrals_erf_cache_max + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,idx) + !DIR$ FORCEINLINE + call map_get(mo_integrals_erf_map,idx,integral) + ii = l-mo_integrals_erf_cache_min + ii = ior( ishft(ii,6), k-mo_integrals_erf_cache_min) + ii = ior( ishft(ii,6), j-mo_integrals_erf_cache_min) + ii = ior( ishft(ii,6), i-mo_integrals_erf_cache_min) + mo_integrals_erf_cache(ii) = integral + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + + +double precision function get_mo_bielec_integral_erf(i,j,k,l,map) + use map_module + implicit none + BEGIN_DOC + ! Returns one integral in the MO basis + END_DOC + integer, intent(in) :: i,j,k,l + integer(key_kind) :: idx + integer :: ii + type(map_type), intent(inout) :: map + real(integral_kind) :: tmp + PROVIDE mo_bielec_integrals_erf_in_map mo_integrals_erf_cache + ii = l-mo_integrals_erf_cache_min + ii = ior(ii, k-mo_integrals_erf_cache_min) + ii = ior(ii, j-mo_integrals_erf_cache_min) + ii = ior(ii, i-mo_integrals_erf_cache_min) + if (iand(ii, -64) /= 0) then + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,idx) + !DIR$ FORCEINLINE + call map_get(map,idx,tmp) + get_mo_bielec_integral_erf = dble(tmp) + else + ii = l-mo_integrals_erf_cache_min + ii = ior( ishft(ii,6), k-mo_integrals_erf_cache_min) + ii = ior( ishft(ii,6), j-mo_integrals_erf_cache_min) + ii = ior( ishft(ii,6), i-mo_integrals_erf_cache_min) + get_mo_bielec_integral_erf = mo_integrals_erf_cache(ii) + endif +end + + +double precision function mo_bielec_integral_erf(i,j,k,l) + implicit none + BEGIN_DOC + ! Returns one integral in the MO basis + END_DOC + integer, intent(in) :: i,j,k,l + double precision :: get_mo_bielec_integral_erf + PROVIDE mo_bielec_integrals_erf_in_map mo_integrals_erf_cache + !DIR$ FORCEINLINE + PROVIDE mo_bielec_integrals_erf_in_map + mo_bielec_integral_erf = get_mo_bielec_integral_erf(i,j,k,l,mo_integrals_erf_map) + return +end + +subroutine get_mo_bielec_integrals_erf(j,k,l,sze,out_val,map) + use map_module + implicit none + BEGIN_DOC + ! Returns multiple integrals in the MO basis, all + ! i for j,k,l fixed. + END_DOC + integer, intent(in) :: j,k,l, sze + double precision, intent(out) :: out_val(sze) + type(map_type), intent(inout) :: map + integer :: i + integer(key_kind) :: hash(sze) + real(integral_kind) :: tmp_val(sze) + PROVIDE mo_bielec_integrals_in_map + + do i=1,sze + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,hash(i)) + enddo + + if (key_kind == 8) then + call map_get_many(map, hash, out_val, sze) + else + call map_get_many(map, hash, tmp_val, sze) + ! Conversion to double precision + do i=1,sze + out_val(i) = dble(tmp_val(i)) + enddo + endif +end + +subroutine get_mo_bielec_integrals_erf_ij(k,l,sze,out_array,map) + use map_module + implicit none + BEGIN_DOC + ! Returns multiple integrals in the MO basis, all + ! i(1)j(2) 1/r12 k(1)l(2) + ! i, j for k,l fixed. + END_DOC + integer, intent(in) :: k,l, sze + double precision, intent(out) :: out_array(sze,sze) + type(map_type), intent(inout) :: map + integer :: i,j,kk,ll,m + integer(key_kind),allocatable :: hash(:) + integer ,allocatable :: pairs(:,:), iorder(:) + real(integral_kind), allocatable :: tmp_val(:) + + PROVIDE mo_bielec_integrals_erf_in_map + allocate (hash(sze*sze), pairs(2,sze*sze),iorder(sze*sze), & + tmp_val(sze*sze)) + + kk=0 + out_array = 0.d0 + do j=1,sze + do i=1,sze + kk += 1 + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,hash(kk)) + pairs(1,kk) = i + pairs(2,kk) = j + iorder(kk) = kk + enddo + enddo + + logical :: integral_is_in_map + if (key_kind == 8) then + call i8radix_sort(hash,iorder,kk,-1) + else if (key_kind == 4) then + call iradix_sort(hash,iorder,kk,-1) + else if (key_kind == 2) then + call i2radix_sort(hash,iorder,kk,-1) + endif + + call map_get_many(mo_integrals_erf_map, hash, tmp_val, kk) + + do ll=1,kk + m = iorder(ll) + i=pairs(1,m) + j=pairs(2,m) + out_array(i,j) = tmp_val(ll) + enddo + + deallocate(pairs,hash,iorder,tmp_val) +end + +subroutine get_mo_bielec_integrals_erf_coulomb_ii(k,l,sze,out_val,map) + use map_module + implicit none + BEGIN_DOC + ! Returns multiple integrals + ! k(1)i(2) 1/r12 l(1)i(2) :: out_val(i1) + ! for k,l fixed. + END_DOC + integer, intent(in) :: k,l, sze + double precision, intent(out) :: out_val(sze) + type(map_type), intent(inout) :: map + integer :: i + integer(key_kind) :: hash(sze) + real(integral_kind) :: tmp_val(sze) + PROVIDE mo_bielec_integrals_in_map + + integer :: kk + do i=1,sze + !DIR$ FORCEINLINE + call bielec_integrals_index(k,i,l,i,hash(i)) + enddo + + if (key_kind == 8) then + call map_get_many(map, hash, out_val, sze) + else + call map_get_many(map, hash, tmp_val, sze) + ! Conversion to double precision + do i=1,sze + out_val(i) = dble(tmp_val(i)) + enddo + endif +end + +subroutine get_mo_bielec_integrals_erf_exch_ii(k,l,sze,out_val,map) + use map_module + implicit none + BEGIN_DOC + ! Returns multiple integrals + ! k(1)i(2) 1/r12 i(1)l(2) :: out_val(i1) + ! for k,l fixed. + END_DOC + integer, intent(in) :: k,l, sze + double precision, intent(out) :: out_val(sze) + type(map_type), intent(inout) :: map + integer :: i + integer(key_kind) :: hash(sze) + real(integral_kind) :: tmp_val(sze) + PROVIDE mo_bielec_integrals_erf_in_map + + integer :: kk + do i=1,sze + !DIR$ FORCEINLINE + call bielec_integrals_index(k,i,i,l,hash(i)) + enddo + + if (key_kind == 8) then + call map_get_many(map, hash, out_val, sze) + else + call map_get_many(map, hash, tmp_val, sze) + ! Conversion to double precision + do i=1,sze + out_val(i) = dble(tmp_val(i)) + enddo + endif +end + + +integer*8 function get_mo_erf_map_size() + implicit none + BEGIN_DOC + ! Return the number of elements in the MO map + END_DOC + get_mo_erf_map_size = mo_integrals_erf_map % n_elements +end diff --git a/plugins/Integrals_erf/mo_bi_integrals.irp.f b/plugins/Integrals_erf/mo_bi_integrals.irp.f new file mode 100644 index 00000000..f4a20cfc --- /dev/null +++ b/plugins/Integrals_erf/mo_bi_integrals.irp.f @@ -0,0 +1,1400 @@ +subroutine mo_bielec_integrals_index(i,j,k,l,i1) + use map_module + implicit none + BEGIN_DOC + ! Computes an unique index for i,j,k,l integrals + END_DOC + integer, intent(in) :: i,j,k,l + integer(key_kind), intent(out) :: i1 + integer(key_kind) :: p,q,r,s,i2 + p = min(i,k) + r = max(i,k) + p = p+ishft(r*r-r,-1) + q = min(j,l) + s = max(j,l) + q = q+ishft(s*s-s,-1) + i1 = min(p,q) + i2 = max(p,q) + i1 = i1+ishft(i2*i2-i2,-1) +end + + +BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ] + use map_module + implicit none + integer(bit_kind) :: mask_ijkl(N_int,4) + integer(bit_kind) :: mask_ijk(N_int,3) + + BEGIN_DOC + ! If True, the map of MO bielectronic integrals is provided + END_DOC + + mo_bielec_integrals_in_map = .True. + if (read_mo_integrals) then + print*,'Reading the MO integrals' + call map_load_from_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map) + print*, 'MO integrals provided' + return + else + PROVIDE ao_bielec_integrals_in_map + endif + + if(no_vvvv_integrals)then + integer :: i,j,k,l + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I I I !!!!!!!!!!!!!!!!!!!! + ! (core+inact+act) ^ 4 + ! + print*, '' + print*, '' + do i = 1,N_int + mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,2) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,3) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,4) = core_inact_act_bitmask_4(i,1) + enddo + call add_integrals_to_map(mask_ijkl) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I V V !!!!!!!!!!!!!!!!!!!! + ! (core+inact+act) ^ 2 (virt) ^2 + ! = J_iv + print*, '' + print*, '' + do i = 1,N_int + mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,2) = virt_bitmask(i,1) + mask_ijkl(i,3) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,4) = virt_bitmask(i,1) + enddo + call add_integrals_to_map(mask_ijkl) + + ! (core+inact+act) ^ 2 (virt) ^2 + ! = (iv|iv) + print*, '' + print*, '' + do i = 1,N_int + mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,2) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,3) = virt_bitmask(i,1) + mask_ijkl(i,4) = virt_bitmask(i,1) + enddo + call add_integrals_to_map(mask_ijkl) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! V V V !!!!!!!!!!!!!!!!!!!!!!! + if(.not.no_vvv_integrals)then + print*, '' + print*, ' and ' + do i = 1,N_int + mask_ijk(i,1) = virt_bitmask(i,1) + mask_ijk(i,2) = virt_bitmask(i,1) + mask_ijk(i,3) = virt_bitmask(i,1) + enddo + call add_integrals_to_map_three_indices(mask_ijk) + endif + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I I V !!!!!!!!!!!!!!!!!!!! + ! (core+inact+act) ^ 3 (virt) ^1 + ! + print*, '' + print*, '' + do i = 1,N_int + mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,2) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,3) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,4) = virt_bitmask(i,1) + enddo + call add_integrals_to_map(mask_ijkl) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I V V V !!!!!!!!!!!!!!!!!!!! + ! (core+inact+act) ^ 1 (virt) ^3 + ! + if(.not.no_ivvv_integrals)then + print*, '' + print*, '' + do i = 1,N_int + mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,2) = virt_bitmask(i,1) + mask_ijkl(i,3) = virt_bitmask(i,1) + mask_ijkl(i,4) = virt_bitmask(i,1) + enddo + call add_integrals_to_map_no_exit_34(mask_ijkl) + endif + + else + call add_integrals_to_map(full_ijkl_bitmask_4) + endif + if (write_mo_integrals) then + call ezfio_set_work_empty(.False.) + call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map) + call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read") + endif + +END_PROVIDER + +subroutine set_integrals_jj_into_map + use bitmasks + implicit none + integer :: i,j,n_integrals,i0,j0 + double precision :: buffer_value(mo_tot_num * mo_tot_num) + integer(key_kind) :: buffer_i(mo_tot_num*mo_tot_num) + n_integrals = 0 + do j0 = 1, n_virt_orb + j = list_virt(j0) + do i0 = j0, n_virt_orb + i = list_virt(i0) + n_integrals += 1 + ! mo_bielec_integral_jj_exchange(i,j) = mo_bielec_integral_vv_exchange_from_ao(i,j) + call mo_bielec_integrals_index(i,j,i,j,buffer_i(n_integrals)) + buffer_value(n_integrals) = mo_bielec_integral_vv_from_ao(i,j) + enddo + enddo + call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + call map_unique(mo_integrals_map) +end + +subroutine set_integrals_exchange_jj_into_map + use bitmasks + implicit none + integer :: i,j,n_integrals,i0,j0 + double precision :: buffer_value(mo_tot_num * mo_tot_num) + integer(key_kind) :: buffer_i(mo_tot_num*mo_tot_num) + n_integrals = 0 + do j0 = 1, n_virt_orb + j = list_virt(j0) + do i0 = j0+1, n_virt_orb + i = list_virt(i0) + n_integrals += 1 + call mo_bielec_integrals_index(i,j,j,i,buffer_i(n_integrals)) + buffer_value(n_integrals) = mo_bielec_integral_vv_exchange_from_ao(i,j) + enddo + enddo + call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + call map_unique(mo_integrals_map) + +end + +subroutine add_integrals_to_map(mask_ijkl) + use bitmasks + implicit none + + BEGIN_DOC + ! Adds integrals to tha MO map according to some bitmask + END_DOC + + integer(bit_kind), intent(in) :: mask_ijkl(N_int,4) + + integer :: i,j,k,l + integer :: i0,j0,k0,l0 + double precision :: c, cpu_1, cpu_2, wall_1, wall_2, wall_0 + + integer, allocatable :: list_ijkl(:,:) + integer :: n_i, n_j, n_k, n_l + integer, allocatable :: bielec_tmp_0_idx(:) + real(integral_kind), allocatable :: bielec_tmp_0(:,:) + double precision, allocatable :: bielec_tmp_1(:) + double precision, allocatable :: bielec_tmp_2(:,:) + double precision, allocatable :: bielec_tmp_3(:,:,:) + !DEC$ ATTRIBUTES ALIGN : 64 :: bielec_tmp_1, bielec_tmp_2, bielec_tmp_3 + + integer :: n_integrals + integer :: size_buffer + integer(key_kind),allocatable :: buffer_i(:) + real(integral_kind),allocatable :: buffer_value(:) + real :: map_mb + + integer :: i1,j1,k1,l1, ii1, kmax, thread_num + integer :: i2,i3,i4 + double precision,parameter :: thr_coef = 1.d-10 + + PROVIDE ao_bielec_integrals_in_map mo_coef + + !Get list of MOs for i,j,k and l + !------------------------------- + + allocate(list_ijkl(mo_tot_num,4)) + call bitstring_to_list( mask_ijkl(1,1), list_ijkl(1,1), n_i, N_int ) + call bitstring_to_list( mask_ijkl(1,2), list_ijkl(1,2), n_j, N_int ) + call bitstring_to_list( mask_ijkl(1,3), list_ijkl(1,3), n_k, N_int ) + call bitstring_to_list( mask_ijkl(1,4), list_ijkl(1,4), n_l, N_int ) + character*(2048) :: output(1) + print*, 'i' + call bitstring_to_str( output(1), mask_ijkl(1,1), N_int ) + print *, trim(output(1)) + j = 0 + do i = 1, N_int + j += popcnt(mask_ijkl(i,1)) + enddo + if(j==0)then + return + endif + + print*, 'j' + call bitstring_to_str( output(1), mask_ijkl(1,2), N_int ) + print *, trim(output(1)) + j = 0 + do i = 1, N_int + j += popcnt(mask_ijkl(i,2)) + enddo + if(j==0)then + return + endif + + print*, 'k' + call bitstring_to_str( output(1), mask_ijkl(1,3), N_int ) + print *, trim(output(1)) + j = 0 + do i = 1, N_int + j += popcnt(mask_ijkl(i,3)) + enddo + if(j==0)then + return + endif + + print*, 'l' + call bitstring_to_str( output(1), mask_ijkl(1,4), N_int ) + print *, trim(output(1)) + j = 0 + do i = 1, N_int + j += popcnt(mask_ijkl(i,4)) + enddo + if(j==0)then + return + endif + + size_buffer = min(ao_num*ao_num*ao_num,16000000) + print*, 'Providing the molecular integrals ' + print*, 'Buffers : ', 8.*(mo_tot_num_align*(n_j)*(n_k+1) + mo_tot_num_align +& + ao_num+ao_num*ao_num+ size_buffer*3)/(1024*1024), 'MB / core' + + call wall_time(wall_1) + call cpu_time(cpu_1) + double precision :: accu_bis + accu_bis = 0.d0 + + !$OMP PARALLEL PRIVATE(l1,k1,j1,i1,i2,i3,i4,i,j,k,l,c, ii1,kmax, & + !$OMP bielec_tmp_0_idx, bielec_tmp_0, bielec_tmp_1,bielec_tmp_2,bielec_tmp_3,& + !$OMP buffer_i,buffer_value,n_integrals,wall_2,i0,j0,k0,l0, & + !$OMP wall_0,thread_num,accu_bis) & + !$OMP DEFAULT(NONE) & + !$OMP SHARED(size_buffer,ao_num,mo_tot_num,n_i,n_j,n_k,n_l,mo_tot_num_align,& + !$OMP mo_coef_transp, & + !$OMP mo_coef_transp_is_built, list_ijkl, & + !$OMP mo_coef_is_built, wall_1, & + !$OMP mo_coef,mo_integrals_threshold,mo_integrals_map) + n_integrals = 0 + wall_0 = wall_1 + allocate(bielec_tmp_3(mo_tot_num_align, n_j, n_k), & + bielec_tmp_1(mo_tot_num_align), & + bielec_tmp_0(ao_num,ao_num), & + bielec_tmp_0_idx(ao_num), & + bielec_tmp_2(mo_tot_num_align, n_j), & + buffer_i(size_buffer), & + buffer_value(size_buffer) ) + + thread_num = 0 + !$ thread_num = omp_get_thread_num() + !$OMP DO SCHEDULE(guided) + do l1 = 1,ao_num + !DEC$ VECTOR ALIGNED + bielec_tmp_3 = 0.d0 + do k1 = 1,ao_num + !DEC$ VECTOR ALIGNED + bielec_tmp_2 = 0.d0 + do j1 = 1,ao_num + call get_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1)) + ! call compute_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1)) + enddo + do j1 = 1,ao_num + kmax = 0 + do i1 = 1,ao_num + c = bielec_tmp_0(i1,j1) + if (c == 0.d0) then + cycle + endif + kmax += 1 + bielec_tmp_0(kmax,j1) = c + bielec_tmp_0_idx(kmax) = i1 + enddo + + if (kmax==0) then + cycle + endif + + !DEC$ VECTOR ALIGNED + bielec_tmp_1 = 0.d0 + ii1=1 + do ii1 = 1,kmax-4,4 + i1 = bielec_tmp_0_idx(ii1) + i2 = bielec_tmp_0_idx(ii1+1) + i3 = bielec_tmp_0_idx(ii1+2) + i4 = bielec_tmp_0_idx(ii1+3) + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_1(i) = bielec_tmp_1(i) + & + mo_coef_transp(i,i1) * bielec_tmp_0(ii1,j1) + & + mo_coef_transp(i,i2) * bielec_tmp_0(ii1+1,j1) + & + mo_coef_transp(i,i3) * bielec_tmp_0(ii1+2,j1) + & + mo_coef_transp(i,i4) * bielec_tmp_0(ii1+3,j1) + enddo ! i + enddo ! ii1 + + i2 = ii1 + do ii1 = i2,kmax + i1 = bielec_tmp_0_idx(ii1) + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_1(i) = bielec_tmp_1(i) + mo_coef_transp(i,i1) * bielec_tmp_0(ii1,j1) + enddo ! i + enddo ! ii1 + c = 0.d0 + + do i = list_ijkl(1,1), list_ijkl(n_i,1) + c = max(c,abs(bielec_tmp_1(i))) + if (c>mo_integrals_threshold) exit + enddo + if ( c < mo_integrals_threshold ) then + cycle + endif + + do j0 = 1, n_j + j = list_ijkl(j0,2) + c = mo_coef_transp(j,j1) + if (abs(c) < thr_coef) then + cycle + endif + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_2(i,j0) = bielec_tmp_2(i,j0) + c * bielec_tmp_1(i) + enddo ! i + enddo ! j + enddo !j1 + if ( maxval(abs(bielec_tmp_2)) < mo_integrals_threshold ) then + cycle + endif + + + do k0 = 1, n_k + k = list_ijkl(k0,3) + c = mo_coef_transp(k,k1) + if (abs(c) < thr_coef) then + cycle + endif + + do j0 = 1, n_j + j = list_ijkl(j0,2) + do i = list_ijkl(1,1), k + bielec_tmp_3(i,j0,k0) = bielec_tmp_3(i,j0,k0) + c* bielec_tmp_2(i,j0) + enddo!i + enddo !j + + enddo !k + enddo !k1 + + + + do l0 = 1,n_l + l = list_ijkl(l0,4) + c = mo_coef_transp(l,l1) + if (abs(c) < thr_coef) then + cycle + endif + j1 = ishft((l*l-l),-1) + do j0 = 1, n_j + j = list_ijkl(j0,2) + if (j > l) then + exit + endif + j1 += 1 + do k0 = 1, n_k + k = list_ijkl(k0,3) + i1 = ishft((k*k-k),-1) + if (i1<=j1) then + continue + else + exit + endif + bielec_tmp_1 = 0.d0 + do i0 = 1, n_i + i = list_ijkl(i0,1) + if (i>k) then + exit + endif + bielec_tmp_1(i) = c*bielec_tmp_3(i,j0,k0) + ! i1+=1 + enddo + + do i0 = 1, n_i + i = list_ijkl(i0,1) + if(i> min(k,j1-i1+list_ijkl(1,1)-1))then + exit + endif + if (abs(bielec_tmp_1(i)) < mo_integrals_threshold) then + cycle + endif + n_integrals += 1 + buffer_value(n_integrals) = bielec_tmp_1(i) + !DEC$ FORCEINLINE + call mo_bielec_integrals_index(i,j,k,l,buffer_i(n_integrals)) + if (n_integrals == size_buffer) then + call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + n_integrals = 0 + endif + enddo + enddo + enddo + enddo + + call wall_time(wall_2) + if (thread_num == 0) then + if (wall_2 - wall_0 > 1.d0) then + wall_0 = wall_2 + print*, 100.*float(l1)/float(ao_num), '% in ', & + wall_2-wall_1, 's', map_mb(mo_integrals_map) ,'MB' + endif + endif + enddo + !$OMP END DO NOWAIT + deallocate (bielec_tmp_1,bielec_tmp_2,bielec_tmp_3) + + integer :: index_needed + + call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + deallocate(buffer_i, buffer_value) + !$OMP END PARALLEL + call map_unique(mo_integrals_map) + + call wall_time(wall_2) + call cpu_time(cpu_2) + integer*8 :: get_mo_map_size, mo_map_size + mo_map_size = get_mo_map_size() + + deallocate(list_ijkl) + + + print*,'Molecular integrals provided:' + print*,' Size of MO map ', map_mb(mo_integrals_map) ,'MB' + print*,' Number of MO integrals: ', mo_map_size + print*,' cpu time :',cpu_2 - cpu_1, 's' + print*,' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')' + +end + + +subroutine add_integrals_to_map_three_indices(mask_ijk) + use bitmasks + implicit none + + BEGIN_DOC + ! Adds integrals to tha MO map according to some bitmask + END_DOC + + integer(bit_kind), intent(in) :: mask_ijk(N_int,3) + + integer :: i,j,k,l + integer :: i0,j0,k0,l0 + double precision :: c, cpu_1, cpu_2, wall_1, wall_2, wall_0 + + integer, allocatable :: list_ijkl(:,:) + integer :: n_i, n_j, n_k + integer :: m + integer, allocatable :: bielec_tmp_0_idx(:) + real(integral_kind), allocatable :: bielec_tmp_0(:,:) + double precision, allocatable :: bielec_tmp_1(:) + double precision, allocatable :: bielec_tmp_2(:,:) + double precision, allocatable :: bielec_tmp_3(:,:,:) + !DEC$ ATTRIBUTES ALIGN : 64 :: bielec_tmp_1, bielec_tmp_2, bielec_tmp_3 + + integer :: n_integrals + integer :: size_buffer + integer(key_kind),allocatable :: buffer_i(:) + real(integral_kind),allocatable :: buffer_value(:) + real :: map_mb + + integer :: i1,j1,k1,l1, ii1, kmax, thread_num + integer :: i2,i3,i4 + double precision,parameter :: thr_coef = 1.d-10 + + PROVIDE ao_bielec_integrals_in_map mo_coef + + !Get list of MOs for i,j,k and l + !------------------------------- + + allocate(list_ijkl(mo_tot_num,4)) + call bitstring_to_list( mask_ijk(1,1), list_ijkl(1,1), n_i, N_int ) + call bitstring_to_list( mask_ijk(1,2), list_ijkl(1,2), n_j, N_int ) + call bitstring_to_list( mask_ijk(1,3), list_ijkl(1,3), n_k, N_int ) + character*(2048) :: output(1) + print*, 'i' + call bitstring_to_str( output(1), mask_ijk(1,1), N_int ) + print *, trim(output(1)) + j = 0 + do i = 1, N_int + j += popcnt(mask_ijk(i,1)) + enddo + if(j==0)then + return + endif + + print*, 'j' + call bitstring_to_str( output(1), mask_ijk(1,2), N_int ) + print *, trim(output(1)) + j = 0 + do i = 1, N_int + j += popcnt(mask_ijk(i,2)) + enddo + if(j==0)then + return + endif + + print*, 'k' + call bitstring_to_str( output(1), mask_ijk(1,3), N_int ) + print *, trim(output(1)) + j = 0 + do i = 1, N_int + j += popcnt(mask_ijk(i,3)) + enddo + if(j==0)then + return + endif + + size_buffer = min(ao_num*ao_num*ao_num,16000000) + print*, 'Providing the molecular integrals ' + print*, 'Buffers : ', 8.*(mo_tot_num_align*(n_j)*(n_k+1) + mo_tot_num_align +& + ao_num+ao_num*ao_num+ size_buffer*3)/(1024*1024), 'MB / core' + + call wall_time(wall_1) + call cpu_time(cpu_1) + double precision :: accu_bis + accu_bis = 0.d0 + !$OMP PARALLEL PRIVATE(m,l1,k1,j1,i1,i2,i3,i4,i,j,k,l,c, ii1,kmax, & + !$OMP bielec_tmp_0_idx, bielec_tmp_0, bielec_tmp_1,bielec_tmp_2,bielec_tmp_3,& + !$OMP buffer_i,buffer_value,n_integrals,wall_2,i0,j0,k0,l0, & + !$OMP wall_0,thread_num,accu_bis) & + !$OMP DEFAULT(NONE) & + !$OMP SHARED(size_buffer,ao_num,mo_tot_num,n_i,n_j,n_k,mo_tot_num_align,& + !$OMP mo_coef_transp, & + !$OMP mo_coef_transp_is_built, list_ijkl, & + !$OMP mo_coef_is_built, wall_1, & + !$OMP mo_coef,mo_integrals_threshold,mo_integrals_map) + n_integrals = 0 + wall_0 = wall_1 + allocate(bielec_tmp_3(mo_tot_num_align, n_j, n_k), & + bielec_tmp_1(mo_tot_num_align), & + bielec_tmp_0(ao_num,ao_num), & + bielec_tmp_0_idx(ao_num), & + bielec_tmp_2(mo_tot_num_align, n_j), & + buffer_i(size_buffer), & + buffer_value(size_buffer) ) + + thread_num = 0 + !$ thread_num = omp_get_thread_num() + !$OMP DO SCHEDULE(guided) + do l1 = 1,ao_num + !DEC$ VECTOR ALIGNED + bielec_tmp_3 = 0.d0 + do k1 = 1,ao_num + !DEC$ VECTOR ALIGNED + bielec_tmp_2 = 0.d0 + do j1 = 1,ao_num + call get_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1)) + enddo + do j1 = 1,ao_num + kmax = 0 + do i1 = 1,ao_num + c = bielec_tmp_0(i1,j1) + if (c == 0.d0) then + cycle + endif + kmax += 1 + bielec_tmp_0(kmax,j1) = c + bielec_tmp_0_idx(kmax) = i1 + enddo + + if (kmax==0) then + cycle + endif + + !DEC$ VECTOR ALIGNED + bielec_tmp_1 = 0.d0 + ii1=1 + do ii1 = 1,kmax-4,4 + i1 = bielec_tmp_0_idx(ii1) + i2 = bielec_tmp_0_idx(ii1+1) + i3 = bielec_tmp_0_idx(ii1+2) + i4 = bielec_tmp_0_idx(ii1+3) + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_1(i) = bielec_tmp_1(i) + & + mo_coef_transp(i,i1) * bielec_tmp_0(ii1,j1) + & + mo_coef_transp(i,i2) * bielec_tmp_0(ii1+1,j1) + & + mo_coef_transp(i,i3) * bielec_tmp_0(ii1+2,j1) + & + mo_coef_transp(i,i4) * bielec_tmp_0(ii1+3,j1) + enddo ! i + enddo ! ii1 + + i2 = ii1 + do ii1 = i2,kmax + i1 = bielec_tmp_0_idx(ii1) + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_1(i) = bielec_tmp_1(i) + mo_coef_transp(i,i1) * bielec_tmp_0(ii1,j1) + enddo ! i + enddo ! ii1 + c = 0.d0 + + do i = list_ijkl(1,1), list_ijkl(n_i,1) + c = max(c,abs(bielec_tmp_1(i))) + if (c>mo_integrals_threshold) exit + enddo + if ( c < mo_integrals_threshold ) then + cycle + endif + + do j0 = 1, n_j + j = list_ijkl(j0,2) + c = mo_coef_transp(j,j1) + if (abs(c) < thr_coef) then + cycle + endif + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_2(i,j0) = bielec_tmp_2(i,j0) + c * bielec_tmp_1(i) + enddo ! i + enddo ! j + enddo !j1 + if ( maxval(abs(bielec_tmp_2)) < mo_integrals_threshold ) then + cycle + endif + + + do k0 = 1, n_k + k = list_ijkl(k0,3) + c = mo_coef_transp(k,k1) + if (abs(c) < thr_coef) then + cycle + endif + + do j0 = 1, n_j + j = list_ijkl(j0,2) + do i = list_ijkl(1,1), k + bielec_tmp_3(i,j0,k0) = bielec_tmp_3(i,j0,k0) + c* bielec_tmp_2(i,j0) + enddo!i + enddo !j + + enddo !k + enddo !k1 + + + + do l0 = 1,n_j + l = list_ijkl(l0,2) + c = mo_coef_transp(l,l1) + if (abs(c) < thr_coef) then + cycle + endif + do k0 = 1, n_k + k = list_ijkl(k0,3) + i1 = ishft((k*k-k),-1) + bielec_tmp_1 = 0.d0 + j0 = l0 + j = list_ijkl(j0,2) + do i0 = 1, n_i + i = list_ijkl(i0,1) + if (i>k) then + exit + endif + bielec_tmp_1(i) = c*bielec_tmp_3(i,j0,k0) + enddo + + do i0 = 1, n_i + i = list_ijkl(i0,1) + if (i>k) then !min(k,j1-i1) + exit + endif + if (abs(bielec_tmp_1(i)) < mo_integrals_threshold) then + cycle + endif + n_integrals += 1 + buffer_value(n_integrals) = bielec_tmp_1(i) + if(i==k .and. j==l .and. i.ne.j)then + buffer_value(n_integrals) = buffer_value(n_integrals) *0.5d0 + endif + !DEC$ FORCEINLINE + call mo_bielec_integrals_index(i,j,k,l,buffer_i(n_integrals)) + if (n_integrals == size_buffer) then + call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + n_integrals = 0 + endif + enddo + enddo + enddo + + do l0 = 1,n_j + l = list_ijkl(l0,2) + c = mo_coef_transp(l,l1) + if (abs(c) < thr_coef) then + cycle + endif + do k0 = 1, n_k + k = list_ijkl(k0,3) + i1 = ishft((k*k-k),-1) + bielec_tmp_1 = 0.d0 + j0 = k0 + j = list_ijkl(k0,2) + i0 = l0 + i = list_ijkl(i0,2) + if (k==l) then + cycle + endif + bielec_tmp_1(i) = c*bielec_tmp_3(i,j0,k0) + + n_integrals += 1 + buffer_value(n_integrals) = bielec_tmp_1(i) + !DEC$ FORCEINLINE + call mo_bielec_integrals_index(i,j,k,l,buffer_i(n_integrals)) + if (n_integrals == size_buffer) then + call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + n_integrals = 0 + endif + enddo + enddo + + call wall_time(wall_2) + if (thread_num == 0) then + if (wall_2 - wall_0 > 1.d0) then + wall_0 = wall_2 + print*, 100.*float(l1)/float(ao_num), '% in ', & + wall_2-wall_1, 's', map_mb(mo_integrals_map) ,'MB' + endif + endif + enddo + !$OMP END DO NOWAIT + deallocate (bielec_tmp_1,bielec_tmp_2,bielec_tmp_3) + + integer :: index_needed + + call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + deallocate(buffer_i, buffer_value) + !$OMP END PARALLEL + call map_unique(mo_integrals_map) + + call wall_time(wall_2) + call cpu_time(cpu_2) + integer*8 :: get_mo_map_size, mo_map_size + mo_map_size = get_mo_map_size() + + deallocate(list_ijkl) + + + print*,'Molecular integrals provided:' + print*,' Size of MO map ', map_mb(mo_integrals_map) ,'MB' + print*,' Number of MO integrals: ', mo_map_size + print*,' cpu time :',cpu_2 - cpu_1, 's' + print*,' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')' + +end + + +subroutine add_integrals_to_map_no_exit_34(mask_ijkl) + use bitmasks + implicit none + + BEGIN_DOC + ! Adds integrals to tha MO map according to some bitmask + END_DOC + + integer(bit_kind), intent(in) :: mask_ijkl(N_int,4) + + integer :: i,j,k,l + integer :: i0,j0,k0,l0 + double precision :: c, cpu_1, cpu_2, wall_1, wall_2, wall_0 + + integer, allocatable :: list_ijkl(:,:) + integer :: n_i, n_j, n_k, n_l + integer, allocatable :: bielec_tmp_0_idx(:) + real(integral_kind), allocatable :: bielec_tmp_0(:,:) + double precision, allocatable :: bielec_tmp_1(:) + double precision, allocatable :: bielec_tmp_2(:,:) + double precision, allocatable :: bielec_tmp_3(:,:,:) + !DEC$ ATTRIBUTES ALIGN : 64 :: bielec_tmp_1, bielec_tmp_2, bielec_tmp_3 + + integer :: n_integrals + integer :: size_buffer + integer(key_kind),allocatable :: buffer_i(:) + real(integral_kind),allocatable :: buffer_value(:) + real :: map_mb + + integer :: i1,j1,k1,l1, ii1, kmax, thread_num + integer :: i2,i3,i4 + double precision,parameter :: thr_coef = 1.d-10 + + PROVIDE ao_bielec_integrals_in_map mo_coef + + !Get list of MOs for i,j,k and l + !------------------------------- + + allocate(list_ijkl(mo_tot_num,4)) + call bitstring_to_list( mask_ijkl(1,1), list_ijkl(1,1), n_i, N_int ) + call bitstring_to_list( mask_ijkl(1,2), list_ijkl(1,2), n_j, N_int ) + call bitstring_to_list( mask_ijkl(1,3), list_ijkl(1,3), n_k, N_int ) + call bitstring_to_list( mask_ijkl(1,4), list_ijkl(1,4), n_l, N_int ) + + size_buffer = min(ao_num*ao_num*ao_num,16000000) + print*, 'Providing the molecular integrals ' + print*, 'Buffers : ', 8.*(mo_tot_num_align*(n_j)*(n_k+1) + mo_tot_num_align +& + ao_num+ao_num*ao_num+ size_buffer*3)/(1024*1024), 'MB / core' + + call wall_time(wall_1) + call cpu_time(cpu_1) + + !$OMP PARALLEL PRIVATE(l1,k1,j1,i1,i2,i3,i4,i,j,k,l,c, ii1,kmax, & + !$OMP bielec_tmp_0_idx, bielec_tmp_0, bielec_tmp_1,bielec_tmp_2,bielec_tmp_3,& + !$OMP buffer_i,buffer_value,n_integrals,wall_2,i0,j0,k0,l0, & + !$OMP wall_0,thread_num) & + !$OMP DEFAULT(NONE) & + !$OMP SHARED(size_buffer,ao_num,mo_tot_num,n_i,n_j,n_k,n_l,mo_tot_num_align,& + !$OMP mo_coef_transp, & + !$OMP mo_coef_transp_is_built, list_ijkl, & + !$OMP mo_coef_is_built, wall_1, & + !$OMP mo_coef,mo_integrals_threshold,mo_integrals_map) + n_integrals = 0 + wall_0 = wall_1 + allocate(bielec_tmp_3(mo_tot_num_align, n_j, n_k), & + bielec_tmp_1(mo_tot_num_align), & + bielec_tmp_0(ao_num,ao_num), & + bielec_tmp_0_idx(ao_num), & + bielec_tmp_2(mo_tot_num_align, n_j), & + buffer_i(size_buffer), & + buffer_value(size_buffer) ) + + thread_num = 0 + !$ thread_num = omp_get_thread_num() + !$OMP DO SCHEDULE(guided) + do l1 = 1,ao_num + !IRP_IF COARRAY + ! if (mod(l1-this_image(),num_images()) /= 0 ) then + ! cycle + ! endif + !IRP_ENDIF + !DEC$ VECTOR ALIGNED + bielec_tmp_3 = 0.d0 + do k1 = 1,ao_num + !DEC$ VECTOR ALIGNED + bielec_tmp_2 = 0.d0 + do j1 = 1,ao_num + call get_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1)) + ! call compute_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1)) + enddo + do j1 = 1,ao_num + kmax = 0 + do i1 = 1,ao_num + c = bielec_tmp_0(i1,j1) + if (c == 0.d0) then + cycle + endif + kmax += 1 + bielec_tmp_0(kmax,j1) = c + bielec_tmp_0_idx(kmax) = i1 + enddo + + if (kmax==0) then + cycle + endif + + !DEC$ VECTOR ALIGNED + bielec_tmp_1 = 0.d0 + ii1=1 + do ii1 = 1,kmax-4,4 + i1 = bielec_tmp_0_idx(ii1) + i2 = bielec_tmp_0_idx(ii1+1) + i3 = bielec_tmp_0_idx(ii1+2) + i4 = bielec_tmp_0_idx(ii1+3) + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_1(i) = bielec_tmp_1(i) + & + mo_coef_transp(i,i1) * bielec_tmp_0(ii1,j1) + & + mo_coef_transp(i,i2) * bielec_tmp_0(ii1+1,j1) + & + mo_coef_transp(i,i3) * bielec_tmp_0(ii1+2,j1) + & + mo_coef_transp(i,i4) * bielec_tmp_0(ii1+3,j1) + enddo ! i + enddo ! ii1 + + i2 = ii1 + do ii1 = i2,kmax + i1 = bielec_tmp_0_idx(ii1) + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_1(i) = bielec_tmp_1(i) + mo_coef_transp(i,i1) * bielec_tmp_0(ii1,j1) + enddo ! i + enddo ! ii1 + c = 0.d0 + + do i = list_ijkl(1,1), list_ijkl(n_i,1) + c = max(c,abs(bielec_tmp_1(i))) + if (c>mo_integrals_threshold) exit + enddo + if ( c < mo_integrals_threshold ) then + cycle + endif + + do j0 = 1, n_j + j = list_ijkl(j0,2) + c = mo_coef_transp(j,j1) + if (abs(c) < thr_coef) then + cycle + endif + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_2(i,j0) = bielec_tmp_2(i,j0) + c * bielec_tmp_1(i) + enddo ! i + enddo ! j + enddo !j1 + if ( maxval(abs(bielec_tmp_2)) < mo_integrals_threshold ) then + cycle + endif + + + do k0 = 1, n_k + k = list_ijkl(k0,3) + c = mo_coef_transp(k,k1) + if (abs(c) < thr_coef) then + cycle + endif + + do j0 = 1, n_j + j = list_ijkl(j0,2) + do i = list_ijkl(1,1), k + bielec_tmp_3(i,j0,k0) = bielec_tmp_3(i,j0,k0) + c* bielec_tmp_2(i,j0) + enddo!i + enddo !j + + enddo !k + enddo !k1 + + + + do l0 = 1,n_l + l = list_ijkl(l0,4) + c = mo_coef_transp(l,l1) + if (abs(c) < thr_coef) then + cycle + endif + j1 = ishft((l*l-l),-1) + do j0 = 1, n_j + j = list_ijkl(j0,2) + if (j > l) then + exit + endif + j1 += 1 + do k0 = 1, n_k + k = list_ijkl(k0,3) + i1 = ishft((k*k-k),-1) + bielec_tmp_1 = 0.d0 + do i0 = 1, n_i + i = list_ijkl(i0,1) + if (i>k) then + exit + endif + bielec_tmp_1(i) = c*bielec_tmp_3(i,j0,k0) + enddo + + do i0 = 1, n_i + i = list_ijkl(i0,1) + if(i> k)then + exit + endif + + if (abs(bielec_tmp_1(i)) < mo_integrals_threshold) then + cycle + endif + n_integrals += 1 + buffer_value(n_integrals) = bielec_tmp_1(i) + !DEC$ FORCEINLINE + call mo_bielec_integrals_index(i,j,k,l,buffer_i(n_integrals)) + if (n_integrals == size_buffer) then + call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + n_integrals = 0 + endif + enddo + enddo + enddo + enddo + + call wall_time(wall_2) + if (thread_num == 0) then + if (wall_2 - wall_0 > 1.d0) then + wall_0 = wall_2 + print*, 100.*float(l1)/float(ao_num), '% in ', & + wall_2-wall_1, 's', map_mb(mo_integrals_map) ,'MB' + endif + endif + enddo + !$OMP END DO NOWAIT + deallocate (bielec_tmp_1,bielec_tmp_2,bielec_tmp_3) + + call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + deallocate(buffer_i, buffer_value) + !$OMP END PARALLEL + !IRP_IF COARRAY + ! print*, 'Communicating the map' + ! call communicate_mo_integrals() + !IRP_ENDIF + call map_unique(mo_integrals_map) + + call wall_time(wall_2) + call cpu_time(cpu_2) + integer*8 :: get_mo_map_size, mo_map_size + mo_map_size = get_mo_map_size() + + deallocate(list_ijkl) + + + print*,'Molecular integrals provided:' + print*,' Size of MO map ', map_mb(mo_integrals_map) ,'MB' + print*,' Number of MO integrals: ', mo_map_size + print*,' cpu time :',cpu_2 - cpu_1, 's' + print*,' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')' + + +end + + + + BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_from_ao, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_exchange_from_ao, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_anti_from_ao, (mo_tot_num_align,mo_tot_num) ] + BEGIN_DOC + ! mo_bielec_integral_jj_from_ao(i,j) = J_ij + ! mo_bielec_integral_jj_exchange_from_ao(i,j) = J_ij + ! mo_bielec_integral_jj_anti_from_ao(i,j) = J_ij - K_ij + END_DOC + implicit none + long_range = .false. + touch long_range +! call routine_mo_bielec_integral_jj_from_ao + +! END_PROVIDER +! +! subroutine routine_mo_bielec_integral_jj_from_ao + integer :: i,j,p,q,r,s + double precision :: c + real(integral_kind) :: integral + integer :: n, pp + real(integral_kind), allocatable :: int_value(:) + integer, allocatable :: int_idx(:) + + double precision, allocatable :: iqrs(:,:), iqsr(:,:), iqis(:), iqri(:) + + if (.not.do_direct_integrals) then + PROVIDE ao_bielec_integrals_in_map mo_coef + endif + + mo_bielec_integral_jj_from_ao = 0.d0 + mo_bielec_integral_jj_exchange_from_ao = 0.d0 + + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: iqrs, iqsr + + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE (i,j,p,q,r,s,integral,c,n,pp,int_value,int_idx, & + !$OMP iqrs, iqsr,iqri,iqis) & + !$OMP SHARED(mo_tot_num,mo_coef_transp,mo_tot_num_align,ao_num,& + !$OMP ao_integrals_threshold,do_direct_integrals) & + !$OMP REDUCTION(+:mo_bielec_integral_jj_from_ao,mo_bielec_integral_jj_exchange_from_ao) + + allocate( int_value(ao_num), int_idx(ao_num), & + iqrs(mo_tot_num_align,ao_num), iqis(mo_tot_num), iqri(mo_tot_num),& + iqsr(mo_tot_num_align,ao_num) ) + + !$OMP DO SCHEDULE (guided) + do s=1,ao_num + do q=1,ao_num + + do j=1,ao_num + !DIR$ VECTOR ALIGNED + do i=1,mo_tot_num + iqrs(i,j) = 0.d0 + iqsr(i,j) = 0.d0 + enddo + enddo + + if (do_direct_integrals) then + double precision :: ao_bielec_integral + do r=1,ao_num + call compute_ao_bielec_integrals(q,r,s,ao_num,int_value) + do p=1,ao_num + integral = int_value(p) + if (abs(integral) > ao_integrals_threshold) then + !DIR$ VECTOR ALIGNED + do i=1,mo_tot_num + iqrs(i,r) += mo_coef_transp(i,p) * integral + enddo + endif + enddo + call compute_ao_bielec_integrals(q,s,r,ao_num,int_value) + do p=1,ao_num + integral = int_value(p) + if (abs(integral) > ao_integrals_threshold) then + !DIR$ VECTOR ALIGNED + do i=1,mo_tot_num + iqsr(i,r) += mo_coef_transp(i,p) * integral + enddo + endif + enddo + enddo + + else + + do r=1,ao_num + call get_ao_bielec_integrals_non_zero(q,r,s,ao_num,int_value,int_idx,n) + do pp=1,n + p = int_idx(pp) + integral = int_value(pp) + if (abs(integral) > ao_integrals_threshold) then + !DIR$ VECTOR ALIGNED + do i=1,mo_tot_num + iqrs(i,r) += mo_coef_transp(i,p) * integral + enddo + endif + enddo + call get_ao_bielec_integrals_non_zero(q,s,r,ao_num,int_value,int_idx,n) + do pp=1,n + p = int_idx(pp) + integral = int_value(pp) + if (abs(integral) > ao_integrals_threshold) then + !DIR$ VECTOR ALIGNED + do i=1,mo_tot_num + iqsr(i,r) += mo_coef_transp(i,p) * integral + enddo + endif + enddo + enddo + endif + iqis = 0.d0 + iqri = 0.d0 + do r=1,ao_num + !DIR$ VECTOR ALIGNED + do i=1,mo_tot_num + iqis(i) += mo_coef_transp(i,r) * iqrs(i,r) + iqri(i) += mo_coef_transp(i,r) * iqsr(i,r) + enddo + enddo + do i=1,mo_tot_num + !DIR$ VECTOR ALIGNED + do j=1,mo_tot_num + c = mo_coef_transp(j,q)*mo_coef_transp(j,s) + mo_bielec_integral_jj_from_ao(j,i) += c * iqis(i) + mo_bielec_integral_jj_exchange_from_ao(j,i) += c * iqri(i) + enddo + enddo + + enddo + enddo + !$OMP END DO NOWAIT + deallocate(iqrs,iqsr,int_value,int_idx) + !$OMP END PARALLEL + + mo_bielec_integral_jj_anti_from_ao = mo_bielec_integral_jj_from_ao - mo_bielec_integral_jj_exchange_from_ao + + +! end +END_PROVIDER + BEGIN_PROVIDER [ double precision, mo_bielec_integral_vv_from_ao, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, mo_bielec_integral_vv_exchange_from_ao, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, mo_bielec_integral_vv_anti_from_ao, (mo_tot_num_align,mo_tot_num) ] + implicit none + long_range = .false. + touch long_range + BEGIN_DOC + ! mo_bielec_integral_vv_from_ao(i,j) = J_ij + ! mo_bielec_integral_vv_exchange_from_ao(i,j) = J_ij + ! mo_bielec_integral_vv_anti_from_ao(i,j) = J_ij - K_ij + ! but only for the virtual orbitals + END_DOC + +! END_PROVIDER + +! subroutine routine_mo_bielec_integral_vv_from_ao + + integer :: i,j,p,q,r,s + integer :: i0,j0 + double precision :: c + real(integral_kind) :: integral + integer :: n, pp + real(integral_kind), allocatable :: int_value(:) + integer, allocatable :: int_idx(:) + + double precision, allocatable :: iqrs(:,:), iqsr(:,:), iqis(:), iqri(:) + + if (.not.do_direct_integrals) then + PROVIDE ao_bielec_integrals_in_map mo_coef + endif + + mo_bielec_integral_vv_from_ao = 0.d0 + mo_bielec_integral_vv_exchange_from_ao = 0.d0 + + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: iqrs, iqsr + + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE (i0,j0,i,j,p,q,r,s,integral,c,n,pp,int_value,int_idx,& + !$OMP iqrs, iqsr,iqri,iqis) & + !$OMP SHARED(n_virt_orb,mo_tot_num,list_virt,mo_coef_transp,mo_tot_num_align,ao_num,& + !$OMP ao_integrals_threshold,do_direct_integrals) & + !$OMP REDUCTION(+:mo_bielec_integral_vv_from_ao,mo_bielec_integral_vv_exchange_from_ao) + + allocate( int_value(ao_num), int_idx(ao_num), & + iqrs(mo_tot_num_align,ao_num), iqis(mo_tot_num), iqri(mo_tot_num),& + iqsr(mo_tot_num_align,ao_num) ) + + !$OMP DO SCHEDULE (guided) + do s=1,ao_num + do q=1,ao_num + + do j=1,ao_num + !DIR$ VECTOR ALIGNED + do i0=1,n_virt_orb + i = list_virt(i0) + iqrs(i,j) = 0.d0 + iqsr(i,j) = 0.d0 + enddo + enddo + + if (do_direct_integrals) then + double precision :: ao_bielec_integral + do r=1,ao_num + call compute_ao_bielec_integrals(q,r,s,ao_num,int_value) + do p=1,ao_num + integral = int_value(p) + if (abs(integral) > ao_integrals_threshold) then + !DIR$ VECTOR ALIGNED + do i0=1,n_virt_orb + i = list_virt(i0) + iqrs(i,r) += mo_coef_transp(i,p) * integral + enddo + endif + enddo + call compute_ao_bielec_integrals(q,s,r,ao_num,int_value) + do p=1,ao_num + integral = int_value(p) + if (abs(integral) > ao_integrals_threshold) then + !DIR$ VECTOR ALIGNED + do i0=1,n_virt_orb + i =list_virt(i0) + iqsr(i,r) += mo_coef_transp(i,p) * integral + enddo + endif + enddo + enddo + + else + + do r=1,ao_num + call get_ao_bielec_integrals_non_zero(q,r,s,ao_num,int_value,int_idx,n) + do pp=1,n + p = int_idx(pp) + integral = int_value(pp) + if (abs(integral) > ao_integrals_threshold) then + !DIR$ VECTOR ALIGNED + do i0=1,n_virt_orb + i =list_virt(i0) + iqrs(i,r) += mo_coef_transp(i,p) * integral + enddo + endif + enddo + call get_ao_bielec_integrals_non_zero(q,s,r,ao_num,int_value,int_idx,n) + do pp=1,n + p = int_idx(pp) + integral = int_value(pp) + if (abs(integral) > ao_integrals_threshold) then + !DIR$ VECTOR ALIGNED + do i0=1,n_virt_orb + i = list_virt(i0) + iqsr(i,r) += mo_coef_transp(i,p) * integral + enddo + endif + enddo + enddo + endif + iqis = 0.d0 + iqri = 0.d0 + do r=1,ao_num + !DIR$ VECTOR ALIGNED + do i0=1,n_virt_orb + i = list_virt(i0) + iqis(i) += mo_coef_transp(i,r) * iqrs(i,r) + iqri(i) += mo_coef_transp(i,r) * iqsr(i,r) + enddo + enddo + do i0=1,n_virt_orb + i= list_virt(i0) + !DIR$ VECTOR ALIGNED + do j0=1,n_virt_orb + j = list_virt(j0) + c = mo_coef_transp(j,q)*mo_coef_transp(j,s) + mo_bielec_integral_vv_from_ao(j,i) += c * iqis(i) + mo_bielec_integral_vv_exchange_from_ao(j,i) += c * iqri(i) + enddo + enddo + + enddo + enddo + !$OMP END DO NOWAIT + deallocate(iqrs,iqsr,int_value,int_idx) + !$OMP END PARALLEL + + mo_bielec_integral_vv_anti_from_ao = mo_bielec_integral_vv_from_ao - mo_bielec_integral_vv_exchange_from_ao + ! print*, '**********' + ! do i0 =1, n_virt_orb + ! i = list_virt(i0) + ! print*, mo_bielec_integral_vv_from_ao(i,i) + ! enddo + ! print*, '**********' + + +!end +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_exchange, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_anti, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! mo_bielec_integral_jj(i,j) = J_ij + ! mo_bielec_integral_jj_exchange(i,j) = K_ij + ! mo_bielec_integral_jj_anti(i,j) = J_ij - K_ij + END_DOC + + integer :: i,j + double precision :: get_mo_bielec_integral + + PROVIDE mo_bielec_integrals_in_map + mo_bielec_integral_jj = 0.d0 + mo_bielec_integral_jj_exchange = 0.d0 + + do j=1,mo_tot_num + do i=1,mo_tot_num + mo_bielec_integral_jj(i,j) = get_mo_bielec_integral(i,j,i,j,mo_integrals_map) + mo_bielec_integral_jj_exchange(i,j) = get_mo_bielec_integral(i,j,j,i,mo_integrals_map) + mo_bielec_integral_jj_anti(i,j) = mo_bielec_integral_jj(i,j) - mo_bielec_integral_jj_exchange(i,j) + enddo + enddo + +END_PROVIDER + + +subroutine clear_mo_map + implicit none + BEGIN_DOC + ! Frees the memory of the MO map + END_DOC + call map_deinit(mo_integrals_map) + FREE mo_integrals_map mo_bielec_integral_jj mo_bielec_integral_jj_anti + FREE mo_bielec_integral_jj_exchange mo_bielec_integrals_in_map + + +end + +subroutine provide_all_mo_integrals + implicit none + provide mo_integrals_map mo_bielec_integral_jj mo_bielec_integral_jj_anti + provide mo_bielec_integral_jj_exchange mo_bielec_integrals_in_map + +end diff --git a/plugins/Integrals_erf/mo_bi_integrals_erf.irp.f b/plugins/Integrals_erf/mo_bi_integrals_erf.irp.f new file mode 100644 index 00000000..c7193819 --- /dev/null +++ b/plugins/Integrals_erf/mo_bi_integrals_erf.irp.f @@ -0,0 +1,616 @@ +subroutine mo_bielec_integrals_erf_index(i,j,k,l,i1) + use map_module + implicit none + BEGIN_DOC + ! Computes an unique index for i,j,k,l integrals + END_DOC + integer, intent(in) :: i,j,k,l + integer(key_kind), intent(out) :: i1 + integer(key_kind) :: p,q,r,s,i2 + p = min(i,k) + r = max(i,k) + p = p+ishft(r*r-r,-1) + q = min(j,l) + s = max(j,l) + q = q+ishft(s*s-s,-1) + i1 = min(p,q) + i2 = max(p,q) + i1 = i1+ishft(i2*i2-i2,-1) +end + + +BEGIN_PROVIDER [ logical, mo_bielec_integrals_erf_in_map ] + use map_module + implicit none + integer(bit_kind) :: mask_ijkl(N_int,4) + integer(bit_kind) :: mask_ijk(N_int,3) + + BEGIN_DOC + ! If True, the map of MO bielectronic integrals is provided + END_DOC + + mo_bielec_integrals_erf_in_map = .True. +! if (read_mo_integrals) then +! print*,'Reading the MO integrals' +! call map_load_from_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map) +! print*, 'MO integrals provided' +! return +! else + PROVIDE ao_bielec_integrals_in_map +! endif + + !if(no_vvvv_integrals)then + ! integer :: i,j,k,l + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I I I !!!!!!!!!!!!!!!!!!!! + ! ! (core+inact+act) ^ 4 + ! ! + ! print*, '' + ! print*, '' + ! do i = 1,N_int + ! mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + ! mask_ijkl(i,2) = core_inact_act_bitmask_4(i,1) + ! mask_ijkl(i,3) = core_inact_act_bitmask_4(i,1) + ! mask_ijkl(i,4) = core_inact_act_bitmask_4(i,1) + ! enddo + ! call add_integrals_to_map(mask_ijkl) + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I V V !!!!!!!!!!!!!!!!!!!! + ! ! (core+inact+act) ^ 2 (virt) ^2 + ! ! = J_iv + ! print*, '' + ! print*, '' + ! do i = 1,N_int + ! mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + ! mask_ijkl(i,2) = virt_bitmask(i,1) + ! mask_ijkl(i,3) = core_inact_act_bitmask_4(i,1) + ! mask_ijkl(i,4) = virt_bitmask(i,1) + ! enddo + ! call add_integrals_to_map(mask_ijkl) + ! + ! ! (core+inact+act) ^ 2 (virt) ^2 + ! ! = (iv|iv) + ! print*, '' + ! print*, '' + ! do i = 1,N_int + ! mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + ! mask_ijkl(i,2) = core_inact_act_bitmask_4(i,1) + ! mask_ijkl(i,3) = virt_bitmask(i,1) + ! mask_ijkl(i,4) = virt_bitmask(i,1) + ! enddo + ! call add_integrals_to_map(mask_ijkl) + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! V V V !!!!!!!!!!!!!!!!!!!!!!! + ! if(.not.no_vvv_integrals)then + ! print*, '' + ! print*, ' and ' + ! do i = 1,N_int + ! mask_ijk(i,1) = virt_bitmask(i,1) + ! mask_ijk(i,2) = virt_bitmask(i,1) + ! mask_ijk(i,3) = virt_bitmask(i,1) + ! enddo + ! call add_integrals_to_map_three_indices(mask_ijk) + ! endif + ! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I I V !!!!!!!!!!!!!!!!!!!! + ! ! (core+inact+act) ^ 3 (virt) ^1 + ! ! + ! print*, '' + ! print*, '' + ! do i = 1,N_int + ! mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + ! mask_ijkl(i,2) = core_inact_act_bitmask_4(i,1) + ! mask_ijkl(i,3) = core_inact_act_bitmask_4(i,1) + ! mask_ijkl(i,4) = virt_bitmask(i,1) + ! enddo + ! call add_integrals_to_map(mask_ijkl) + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I V V V !!!!!!!!!!!!!!!!!!!! + ! ! (core+inact+act) ^ 1 (virt) ^3 + ! ! + ! if(.not.no_ivvv_integrals)then + ! print*, '' + ! print*, '' + ! do i = 1,N_int + ! mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + ! mask_ijkl(i,2) = virt_bitmask(i,1) + ! mask_ijkl(i,3) = virt_bitmask(i,1) + ! mask_ijkl(i,4) = virt_bitmask(i,1) + ! enddo + ! call add_integrals_to_map_no_exit_34(mask_ijkl) + ! endif + ! + !else + call add_integrals_to_map(full_ijkl_bitmask_4) + !endif + !if (write_mo_integrals) then + ! call ezfio_set_work_empty(.False.) + ! call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map) + ! call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read") + !endif + +END_PROVIDER + +subroutine add_integrals_erf_to_map(mask_ijkl) + use bitmasks + implicit none + + BEGIN_DOC + ! Adds integrals to tha MO map according to some bitmask + END_DOC + + integer(bit_kind), intent(in) :: mask_ijkl(N_int,4) + + integer :: i,j,k,l + integer :: i0,j0,k0,l0 + double precision :: c, cpu_1, cpu_2, wall_1, wall_2, wall_0 + + integer, allocatable :: list_ijkl(:,:) + integer :: n_i, n_j, n_k, n_l + integer, allocatable :: bielec_tmp_0_idx(:) + real(integral_kind), allocatable :: bielec_tmp_0(:,:) + double precision, allocatable :: bielec_tmp_1(:) + double precision, allocatable :: bielec_tmp_2(:,:) + double precision, allocatable :: bielec_tmp_3(:,:,:) + !DEC$ ATTRIBUTES ALIGN : 64 :: bielec_tmp_1, bielec_tmp_2, bielec_tmp_3 + + integer :: n_integrals + integer :: size_buffer + integer(key_kind),allocatable :: buffer_i(:) + real(integral_kind),allocatable :: buffer_value(:) + real :: map_mb + + integer :: i1,j1,k1,l1, ii1, kmax, thread_num + integer :: i2,i3,i4 + double precision,parameter :: thr_coef = 1.d-10 + + PROVIDE ao_bielec_integrals_in_map mo_coef + + !Get list of MOs for i,j,k and l + !------------------------------- + + allocate(list_ijkl(mo_tot_num,4)) + call bitstring_to_list( mask_ijkl(1,1), list_ijkl(1,1), n_i, N_int ) + call bitstring_to_list( mask_ijkl(1,2), list_ijkl(1,2), n_j, N_int ) + call bitstring_to_list( mask_ijkl(1,3), list_ijkl(1,3), n_k, N_int ) + call bitstring_to_list( mask_ijkl(1,4), list_ijkl(1,4), n_l, N_int ) + character*(2048) :: output(1) + print*, 'i' + call bitstring_to_str( output(1), mask_ijkl(1,1), N_int ) + print *, trim(output(1)) + j = 0 + do i = 1, N_int + j += popcnt(mask_ijkl(i,1)) + enddo + if(j==0)then + return + endif + + print*, 'j' + call bitstring_to_str( output(1), mask_ijkl(1,2), N_int ) + print *, trim(output(1)) + j = 0 + do i = 1, N_int + j += popcnt(mask_ijkl(i,2)) + enddo + if(j==0)then + return + endif + + print*, 'k' + call bitstring_to_str( output(1), mask_ijkl(1,3), N_int ) + print *, trim(output(1)) + j = 0 + do i = 1, N_int + j += popcnt(mask_ijkl(i,3)) + enddo + if(j==0)then + return + endif + + print*, 'l' + call bitstring_to_str( output(1), mask_ijkl(1,4), N_int ) + print *, trim(output(1)) + j = 0 + do i = 1, N_int + j += popcnt(mask_ijkl(i,4)) + enddo + if(j==0)then + return + endif + + size_buffer = min(ao_num*ao_num*ao_num,16000000) + print*, 'Providing the molecular integrals ' + print*, 'Buffers : ', 8.*(mo_tot_num_align*(n_j)*(n_k+1) + mo_tot_num_align +& + ao_num+ao_num*ao_num+ size_buffer*3)/(1024*1024), 'MB / core' + + call wall_time(wall_1) + call cpu_time(cpu_1) + double precision :: accu_bis + accu_bis = 0.d0 + + !$OMP PARALLEL PRIVATE(l1,k1,j1,i1,i2,i3,i4,i,j,k,l,c, ii1,kmax, & + !$OMP bielec_tmp_0_idx, bielec_tmp_0, bielec_tmp_1,bielec_tmp_2,bielec_tmp_3,& + !$OMP buffer_i,buffer_value,n_integrals,wall_2,i0,j0,k0,l0, & + !$OMP wall_0,thread_num,accu_bis) & + !$OMP DEFAULT(NONE) & + !$OMP SHARED(size_buffer,ao_num,mo_tot_num,n_i,n_j,n_k,n_l,mo_tot_num_align,& + !$OMP mo_coef_transp, & + !$OMP mo_coef_transp_is_built, list_ijkl, & + !$OMP mo_coef_is_built, wall_1, & + !$OMP mo_coef,mo_integrals_threshold,mo_integrals_erf_map) + n_integrals = 0 + wall_0 = wall_1 + allocate(bielec_tmp_3(mo_tot_num_align, n_j, n_k), & + bielec_tmp_1(mo_tot_num_align), & + bielec_tmp_0(ao_num,ao_num), & + bielec_tmp_0_idx(ao_num), & + bielec_tmp_2(mo_tot_num_align, n_j), & + buffer_i(size_buffer), & + buffer_value(size_buffer) ) + + thread_num = 0 + !$ thread_num = omp_get_thread_num() + !$OMP DO SCHEDULE(guided) + do l1 = 1,ao_num + !DEC$ VECTOR ALIGNED + bielec_tmp_3 = 0.d0 + do k1 = 1,ao_num + !DEC$ VECTOR ALIGNED + bielec_tmp_2 = 0.d0 + do j1 = 1,ao_num + call get_ao_bielec_integrals_erf(j1,k1,l1,ao_num,bielec_tmp_0(1,j1)) + ! call compute_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1)) + enddo + do j1 = 1,ao_num + kmax = 0 + do i1 = 1,ao_num + c = bielec_tmp_0(i1,j1) + if (c == 0.d0) then + cycle + endif + kmax += 1 + bielec_tmp_0(kmax,j1) = c + bielec_tmp_0_idx(kmax) = i1 + enddo + + if (kmax==0) then + cycle + endif + + !DEC$ VECTOR ALIGNED + bielec_tmp_1 = 0.d0 + ii1=1 + do ii1 = 1,kmax-4,4 + i1 = bielec_tmp_0_idx(ii1) + i2 = bielec_tmp_0_idx(ii1+1) + i3 = bielec_tmp_0_idx(ii1+2) + i4 = bielec_tmp_0_idx(ii1+3) + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_1(i) = bielec_tmp_1(i) + & + mo_coef_transp(i,i1) * bielec_tmp_0(ii1,j1) + & + mo_coef_transp(i,i2) * bielec_tmp_0(ii1+1,j1) + & + mo_coef_transp(i,i3) * bielec_tmp_0(ii1+2,j1) + & + mo_coef_transp(i,i4) * bielec_tmp_0(ii1+3,j1) + enddo ! i + enddo ! ii1 + + i2 = ii1 + do ii1 = i2,kmax + i1 = bielec_tmp_0_idx(ii1) + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_1(i) = bielec_tmp_1(i) + mo_coef_transp(i,i1) * bielec_tmp_0(ii1,j1) + enddo ! i + enddo ! ii1 + c = 0.d0 + + do i = list_ijkl(1,1), list_ijkl(n_i,1) + c = max(c,abs(bielec_tmp_1(i))) + if (c>mo_integrals_threshold) exit + enddo + if ( c < mo_integrals_threshold ) then + cycle + endif + + do j0 = 1, n_j + j = list_ijkl(j0,2) + c = mo_coef_transp(j,j1) + if (abs(c) < thr_coef) then + cycle + endif + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_2(i,j0) = bielec_tmp_2(i,j0) + c * bielec_tmp_1(i) + enddo ! i + enddo ! j + enddo !j1 + if ( maxval(abs(bielec_tmp_2)) < mo_integrals_threshold ) then + cycle + endif + + + do k0 = 1, n_k + k = list_ijkl(k0,3) + c = mo_coef_transp(k,k1) + if (abs(c) < thr_coef) then + cycle + endif + + do j0 = 1, n_j + j = list_ijkl(j0,2) + do i = list_ijkl(1,1), k + bielec_tmp_3(i,j0,k0) = bielec_tmp_3(i,j0,k0) + c* bielec_tmp_2(i,j0) + enddo!i + enddo !j + + enddo !k + enddo !k1 + + + + do l0 = 1,n_l + l = list_ijkl(l0,4) + c = mo_coef_transp(l,l1) + if (abs(c) < thr_coef) then + cycle + endif + j1 = ishft((l*l-l),-1) + do j0 = 1, n_j + j = list_ijkl(j0,2) + if (j > l) then + exit + endif + j1 += 1 + do k0 = 1, n_k + k = list_ijkl(k0,3) + i1 = ishft((k*k-k),-1) + if (i1<=j1) then + continue + else + exit + endif + bielec_tmp_1 = 0.d0 + do i0 = 1, n_i + i = list_ijkl(i0,1) + if (i>k) then + exit + endif + bielec_tmp_1(i) = c*bielec_tmp_3(i,j0,k0) + ! i1+=1 + enddo + + do i0 = 1, n_i + i = list_ijkl(i0,1) + if(i> min(k,j1-i1+list_ijkl(1,1)-1))then + exit + endif + if (abs(bielec_tmp_1(i)) < mo_integrals_threshold) then + cycle + endif + n_integrals += 1 + buffer_value(n_integrals) = bielec_tmp_1(i) + !DEC$ FORCEINLINE + call mo_bielec_integrals_erf_index(i,j,k,l,buffer_i(n_integrals)) + if (n_integrals == size_buffer) then + call insert_into_mo_integrals_erf_map(n_integrals,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + n_integrals = 0 + endif + enddo + enddo + enddo + enddo + + call wall_time(wall_2) + if (thread_num == 0) then + if (wall_2 - wall_0 > 1.d0) then + wall_0 = wall_2 + print*, 100.*float(l1)/float(ao_num), '% in ', & + wall_2-wall_1, 's', map_mb(mo_integrals_erf_map) ,'MB' + endif + endif + enddo + !$OMP END DO NOWAIT + deallocate (bielec_tmp_1,bielec_tmp_2,bielec_tmp_3) + + integer :: index_needed + + call insert_into_mo_integrals_erf_map(n_integrals,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + deallocate(buffer_i, buffer_value) + !$OMP END PARALLEL + call map_unique(mo_integrals_erf_map) + + call wall_time(wall_2) + call cpu_time(cpu_2) + integer*8 :: get_mo_erf_map_size, mo_erf_map_size + mo_erf_map_size = get_mo_erf_map_size() + + deallocate(list_ijkl) + + + print*,'Molecular integrals provided:' + print*,' Size of MO map ', map_mb(mo_integrals_erf_map) ,'MB' + print*,' Number of MO integrals: ', mo_erf_map_size + print*,' cpu time :',cpu_2 - cpu_1, 's' + print*,' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')' + +end + + + + BEGIN_PROVIDER [ double precision, mo_bielec_integral_erf_jj_from_ao, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, mo_bielec_integral_erf_jj_exchange_from_ao, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, mo_bielec_integral_erf_jj_anti_from_ao, (mo_tot_num_align,mo_tot_num) ] + BEGIN_DOC + ! mo_bielec_integral_jj_from_ao(i,j) = J_ij + ! mo_bielec_integral_jj_exchange_from_ao(i,j) = J_ij + ! mo_bielec_integral_jj_anti_from_ao(i,j) = J_ij - K_ij + END_DOC + implicit none + integer :: i,j,p,q,r,s + double precision :: c + real(integral_kind) :: integral + integer :: n, pp + real(integral_kind), allocatable :: int_value(:) + integer, allocatable :: int_idx(:) + + double precision, allocatable :: iqrs(:,:), iqsr(:,:), iqis(:), iqri(:) + + if (.not.do_direct_integrals) then + PROVIDE ao_bielec_integrals_erf_in_map mo_coef + endif + + mo_bielec_integral_erf_jj_from_ao = 0.d0 + mo_bielec_integral_erf_jj_exchange_from_ao = 0.d0 + + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: iqrs, iqsr + + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE (i,j,p,q,r,s,integral,c,n,pp,int_value,int_idx, & + !$OMP iqrs, iqsr,iqri,iqis) & + !$OMP SHARED(mo_tot_num,mo_coef_transp,mo_tot_num_align,ao_num,& + !$OMP ao_integrals_threshold,do_direct_integrals) & + !$OMP REDUCTION(+:mo_bielec_integral_erf_jj_from_ao,mo_bielec_integral_erf_jj_exchange_from_ao) + + allocate( int_value(ao_num), int_idx(ao_num), & + iqrs(mo_tot_num_align,ao_num), iqis(mo_tot_num), iqri(mo_tot_num),& + iqsr(mo_tot_num_align,ao_num) ) + + !$OMP DO SCHEDULE (guided) + do s=1,ao_num + do q=1,ao_num + + do j=1,ao_num + !DIR$ VECTOR ALIGNED + do i=1,mo_tot_num + iqrs(i,j) = 0.d0 + iqsr(i,j) = 0.d0 + enddo + enddo + + if (do_direct_integrals) then + double precision :: ao_bielec_integral_erf + do r=1,ao_num + call compute_ao_bielec_integrals_erf(q,r,s,ao_num,int_value) + do p=1,ao_num + integral = int_value(p) + if (abs(integral) > ao_integrals_threshold) then + !DIR$ VECTOR ALIGNED + do i=1,mo_tot_num + iqrs(i,r) += mo_coef_transp(i,p) * integral + enddo + endif + enddo + call compute_ao_bielec_integrals_erf(q,s,r,ao_num,int_value) + do p=1,ao_num + integral = int_value(p) + if (abs(integral) > ao_integrals_threshold) then + !DIR$ VECTOR ALIGNED + do i=1,mo_tot_num + iqsr(i,r) += mo_coef_transp(i,p) * integral + enddo + endif + enddo + enddo + + else + + do r=1,ao_num + call get_ao_bielec_integrals_erf_non_zero(q,r,s,ao_num,int_value,int_idx,n) + do pp=1,n + p = int_idx(pp) + integral = int_value(pp) + if (abs(integral) > ao_integrals_threshold) then + !DIR$ VECTOR ALIGNED + do i=1,mo_tot_num + iqrs(i,r) += mo_coef_transp(i,p) * integral + enddo + endif + enddo + call get_ao_bielec_integrals_erf_non_zero(q,s,r,ao_num,int_value,int_idx,n) + do pp=1,n + p = int_idx(pp) + integral = int_value(pp) + if (abs(integral) > ao_integrals_threshold) then + !DIR$ VECTOR ALIGNED + do i=1,mo_tot_num + iqsr(i,r) += mo_coef_transp(i,p) * integral + enddo + endif + enddo + enddo + endif + iqis = 0.d0 + iqri = 0.d0 + do r=1,ao_num + !DIR$ VECTOR ALIGNED + do i=1,mo_tot_num + iqis(i) += mo_coef_transp(i,r) * iqrs(i,r) + iqri(i) += mo_coef_transp(i,r) * iqsr(i,r) + enddo + enddo + do i=1,mo_tot_num + !DIR$ VECTOR ALIGNED + do j=1,mo_tot_num + c = mo_coef_transp(j,q)*mo_coef_transp(j,s) + mo_bielec_integral_erf_jj_from_ao(j,i) += c * iqis(i) + mo_bielec_integral_erf_jj_exchange_from_ao(j,i) += c * iqri(i) + enddo + enddo + + enddo + enddo + !$OMP END DO NOWAIT + deallocate(iqrs,iqsr,int_value,int_idx) + !$OMP END PARALLEL + + mo_bielec_integral_erf_jj_anti_from_ao = mo_bielec_integral_erf_jj_from_ao - mo_bielec_integral_erf_jj_exchange_from_ao + + +! end +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, mo_bielec_integral_erf_jj, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, mo_bielec_integral_erf_jj_exchange, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, mo_bielec_integral_erf_jj_anti, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! mo_bielec_integral_jj(i,j) = J_ij + ! mo_bielec_integral_jj_exchange(i,j) = K_ij + ! mo_bielec_integral_jj_anti(i,j) = J_ij - K_ij + END_DOC + + integer :: i,j + double precision :: get_mo_bielec_integral_erf + + PROVIDE mo_bielec_integrals_erf_in_map + mo_bielec_integral_erf_jj = 0.d0 + mo_bielec_integral_erf_jj_exchange = 0.d0 + + do j=1,mo_tot_num + do i=1,mo_tot_num + mo_bielec_integral_erf_jj(i,j) = get_mo_bielec_integral_erf(i,j,i,j,mo_integrals_map) + mo_bielec_integral_erf_jj_exchange(i,j) = get_mo_bielec_integral_erf(i,j,j,i,mo_integrals_map) + mo_bielec_integral_erf_jj_anti(i,j) = mo_bielec_integral_erf_jj(i,j) - mo_bielec_integral_erf_jj_exchange(i,j) + enddo + enddo + +END_PROVIDER + + +subroutine clear_mo_erf_map + implicit none + BEGIN_DOC + ! Frees the memory of the MO map + END_DOC + call map_deinit(mo_integrals_erf_map) + FREE mo_integrals_erf_map mo_bielec_integral_erf_jj mo_bielec_integral_erf_jj_anti + FREE mo_bielec_integral_Erf_jj_exchange mo_bielec_integrals_erf_in_map + + +end + +subroutine provide_all_mo_integrals_erf + implicit none + provide mo_integrals_erf_map mo_bielec_integral_erf_jj mo_bielec_integral_erf_jj_anti + provide mo_bielec_integral_erf_jj_exchange mo_bielec_integrals_erf_in_map + +end diff --git a/plugins/Integrals_erf/providers_ao_erf.irp.f b/plugins/Integrals_erf/providers_ao_erf.irp.f new file mode 100644 index 00000000..cff3e49c --- /dev/null +++ b/plugins/Integrals_erf/providers_ao_erf.irp.f @@ -0,0 +1,119 @@ + +BEGIN_PROVIDER [ logical, ao_bielec_integrals_erf_in_map ] + implicit none + use f77_zmq + use map_module + BEGIN_DOC + ! Map of Atomic integrals + ! i(r1) j(r2) 1/r12 k(r1) l(r2) + END_DOC + + integer :: i,j,k,l + double precision :: ao_bielec_integral_erf,cpu_1,cpu_2, wall_1, wall_2 + double precision :: integral, wall_0 + include 'Utils/constants.include.F' + + ! For integrals file + integer(key_kind),allocatable :: buffer_i(:) + integer,parameter :: size_buffer = 1024*64 + real(integral_kind),allocatable :: buffer_value(:) + + integer :: n_integrals, rc + integer :: kk, m, j1, i1, lmax + character*(64) :: fmt + + integral = ao_bielec_integral_erf(1,1,1,1) + + real :: map_mb +! PROVIDE read_ao_integrals disk_access_ao_integrals +! if (read_ao_integrals) then +! print*,'Reading the AO integrals' +! call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) +! print*, 'AO integrals provided' +! ao_bielec_integrals_in_map = .True. +! return +! endif + + print*, 'Providing the AO integrals' + call wall_time(wall_0) + call wall_time(wall_1) + call cpu_time(cpu_1) + + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + call new_parallel_job(zmq_to_qp_run_socket,'ao_integrals') + + character(len=:), allocatable :: task + allocate(character(len=ao_num*12) :: task) + write(fmt,*) '(', ao_num, '(I5,X,I5,''|''))' + do l=1,ao_num + write(task,fmt) (i,l, i=1,l) + call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task)) + enddo + deallocate(task) + + call zmq_set_running(zmq_to_qp_run_socket) + + PROVIDE nproc + !$OMP PARALLEL DEFAULT(private) num_threads(nproc+1) + i = omp_get_thread_num() + if (i==0) then + call ao_bielec_integrals_erf_in_map_collector(i) + else + call ao_bielec_integrals_erf_in_map_slave_inproc(i) + endif + !$OMP END PARALLEL + + call end_parallel_job(zmq_to_qp_run_socket, 'ao_integrals') + + + print*, 'Sorting the map' + call map_sort(ao_integrals_erf_map) + call cpu_time(cpu_2) + call wall_time(wall_2) + integer(map_size_kind) :: get_ao_erf_map_size, ao_erf_map_size + ao_erf_map_size = get_ao_erf_map_size() + + print*, 'AO integrals provided:' + print*, ' Size of AO map : ', map_mb(ao_integrals_erf_map) ,'MB' + print*, ' Number of AO integrals :', ao_erf_map_size + print*, ' cpu time :',cpu_2 - cpu_1, 's' + print*, ' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1+tiny(1.d0)), ' )' + + ao_bielec_integrals_erf_in_map = .True. + +! if (write_ao_integrals) then +! call ezfio_set_work_empty(.False.) +! call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_erf_map) +! call ezfio_set_integrals_bielec_disk_access_ao_integrals("Read") +! endif + +END_PROVIDER + + + + +BEGIN_PROVIDER [ double precision, ao_bielec_integral_erf_schwartz,(ao_num,ao_num) ] + implicit none + BEGIN_DOC + ! Needed to compute Schwartz inequalities + END_DOC + + integer :: i,k + double precision :: ao_bielec_integral_erf,cpu_1,cpu_2, wall_1, wall_2 + + ao_bielec_integral_erf_schwartz(1,1) = ao_bielec_integral_erf(1,1,1,1) + !$OMP PARALLEL DO PRIVATE(i,k) & + !$OMP DEFAULT(NONE) & + !$OMP SHARED (ao_num,ao_bielec_integral_erf_schwartz) & + !$OMP SCHEDULE(dynamic) + do i=1,ao_num + do k=1,i + ao_bielec_integral_erf_schwartz(i,k) = dsqrt(ao_bielec_integral_erf(i,k,i,k)) + ao_bielec_integral_erf_schwartz(k,i) = ao_bielec_integral_erf_schwartz(i,k) + enddo + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + + diff --git a/plugins/Integrals_erf/providers_ao_standard.irp.f b/plugins/Integrals_erf/providers_ao_standard.irp.f new file mode 100644 index 00000000..c0b043da --- /dev/null +++ b/plugins/Integrals_erf/providers_ao_standard.irp.f @@ -0,0 +1,116 @@ + +BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] + implicit none + use f77_zmq + use map_module + BEGIN_DOC + ! Map of Atomic integrals + ! i(r1) j(r2) 1/r12 k(r1) l(r2) + END_DOC + + integer :: i,j,k,l + double precision :: ao_bielec_integral,cpu_1,cpu_2, wall_1, wall_2 + double precision :: integral, wall_0 + include 'Utils/constants.include.F' + + ! For integrals file + integer(key_kind),allocatable :: buffer_i(:) + integer,parameter :: size_buffer = 1024*64 + real(integral_kind),allocatable :: buffer_value(:) + + integer :: n_integrals, rc + integer :: kk, m, j1, i1, lmax + character*(64) :: fmt + + integral = ao_bielec_integral(1,1,1,1) + + real :: map_mb + PROVIDE read_ao_integrals disk_access_ao_integrals + if (read_ao_integrals) then + print*,'Reading the AO integrals' + call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) + print*, 'AO integrals provided' + ao_bielec_integrals_in_map = .True. + return + endif + + print*, 'Providing the AO integrals' + call wall_time(wall_0) + call wall_time(wall_1) + call cpu_time(cpu_1) + + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + call new_parallel_job(zmq_to_qp_run_socket,'ao_integrals') + + character(len=:), allocatable :: task + allocate(character(len=ao_num*12) :: task) + write(fmt,*) '(', ao_num, '(I5,X,I5,''|''))' + do l=1,ao_num + write(task,fmt) (i,l, i=1,l) + call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task)) + enddo + deallocate(task) + + call zmq_set_running(zmq_to_qp_run_socket) + + PROVIDE nproc + !$OMP PARALLEL DEFAULT(private) num_threads(nproc+1) + i = omp_get_thread_num() + if (i==0) then + call ao_bielec_integrals_in_map_collector(i) + else + call ao_bielec_integrals_in_map_slave_inproc(i) + endif + !$OMP END PARALLEL + + call end_parallel_job(zmq_to_qp_run_socket, 'ao_integrals') + + + print*, 'Sorting the map' + call map_sort(ao_integrals_map) + call cpu_time(cpu_2) + call wall_time(wall_2) + integer(map_size_kind) :: get_ao_map_size, ao_map_size + ao_map_size = get_ao_map_size() + + print*, 'AO integrals provided:' + print*, ' Size of AO map : ', map_mb(ao_integrals_map) ,'MB' + print*, ' Number of AO integrals :', ao_map_size + print*, ' cpu time :',cpu_2 - cpu_1, 's' + print*, ' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1+tiny(1.d0)), ' )' + + ao_bielec_integrals_in_map = .True. + + if (write_ao_integrals) then + call ezfio_set_work_empty(.False.) + call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) + call ezfio_set_integrals_bielec_disk_access_ao_integrals("Read") + endif + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_bielec_integral_schwartz,(ao_num,ao_num) ] + implicit none + BEGIN_DOC + ! Needed to compute Schwartz inequalities + END_DOC + + integer :: i,k + double precision :: ao_bielec_integral,cpu_1,cpu_2, wall_1, wall_2 + + ao_bielec_integral_schwartz(1,1) = ao_bielec_integral(1,1,1,1) + !$OMP PARALLEL DO PRIVATE(i,k) & + !$OMP DEFAULT(NONE) & + !$OMP SHARED (ao_num,ao_bielec_integral_schwartz) & + !$OMP SCHEDULE(dynamic) + do i=1,ao_num + do k=1,i + ao_bielec_integral_schwartz(i,k) = dsqrt(ao_bielec_integral(i,k,i,k)) + ao_bielec_integral_schwartz(k,i) = ao_bielec_integral_schwartz(i,k) + enddo + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + + diff --git a/plugins/Integrals_erf/qp_ao_ints.irp.f b/plugins/Integrals_erf/qp_ao_ints.irp.f new file mode 100644 index 00000000..93f62a7d --- /dev/null +++ b/plugins/Integrals_erf/qp_ao_ints.irp.f @@ -0,0 +1,32 @@ +program qp_ao_ints + use omp_lib + implicit none + BEGIN_DOC +! Increments a running calculation to compute AO integrals + END_DOC + integer :: i + + call switch_qp_run_to_master + + zmq_context = f77_zmq_ctx_new () + + ! Set the state of the ZMQ + zmq_state = 'ao_integrals' + + ! Provide everything needed + double precision :: integral, ao_bielec_integral + integral = ao_bielec_integral(1,1,1,1) + + character*(64) :: state + call wait_for_state(zmq_state,state) + do while (state /= 'Stopped') + !$OMP PARALLEL DEFAULT(PRIVATE) PRIVATE(i) + i = omp_get_thread_num() + call ao_bielec_integrals_in_map_slave_tcp(i) + !$OMP END PARALLEL + call wait_for_state(zmq_state,state) + enddo + + print *, 'Done' +end + diff --git a/plugins/Integrals_erf/read_write.irp.f b/plugins/Integrals_erf/read_write.irp.f new file mode 100644 index 00000000..5b2b7f3e --- /dev/null +++ b/plugins/Integrals_erf/read_write.irp.f @@ -0,0 +1,47 @@ +BEGIN_PROVIDER [ logical, read_ao_integrals ] +&BEGIN_PROVIDER [ logical, read_mo_integrals ] +&BEGIN_PROVIDER [ logical, write_ao_integrals ] +&BEGIN_PROVIDER [ logical, write_mo_integrals ] + + BEGIN_DOC +! One level of abstraction for disk_access_ao_integrals and disk_access_mo_integrals + END_DOC +implicit none + + if (disk_access_ao_integrals.EQ.'Read') then + read_ao_integrals = .True. + write_ao_integrals = .False. + + else if (disk_access_ao_integrals.EQ.'Write') then + read_ao_integrals = .False. + write_ao_integrals = .True. + + else if (disk_access_ao_integrals.EQ.'None') then + read_ao_integrals = .False. + write_ao_integrals = .False. + + else + print *, 'bielec_integrals/disk_access_ao_integrals has a wrong type' + stop 1 + + endif + + if (disk_access_mo_integrals.EQ.'Read') then + read_mo_integrals = .True. + write_mo_integrals = .False. + + else if (disk_access_mo_integrals.EQ.'Write') then + read_mo_integrals = .False. + write_mo_integrals = .True. + + else if (disk_access_mo_integrals.EQ.'None') then + read_mo_integrals = .False. + write_mo_integrals = .False. + + else + print *, 'bielec_integrals/disk_access_mo_integrals has a wrong type' + stop 1 + + endif + +END_PROVIDER diff --git a/plugins/Integrals_erf/test_bitmasks_integrals.irp.f b/plugins/Integrals_erf/test_bitmasks_integrals.irp.f new file mode 100644 index 00000000..d0eba3f7 --- /dev/null +++ b/plugins/Integrals_erf/test_bitmasks_integrals.irp.f @@ -0,0 +1,33 @@ +program pouet + implicit none + + call routine + +end + +subroutine routine + implicit none + integer(bit_kind) :: mask_ijkl(N_int,4) + integer, allocatable :: list_ijkl(:,:) + integer :: i,j + integer :: n_i,n_j,n_k,n_l + do i = 1,N_int + mask_ijkl(i,1) = inact_bitmask(i,1) + mask_ijkl(i,2) = inact_bitmask(i,1) + mask_ijkl(i,3) = inact_bitmask(i,1) + mask_ijkl(i,4) = inact_bitmask(i,1) + enddo + allocate(list_ijkl(mo_tot_num,4)) + call bitstring_to_list( mask_ijkl(1,1), list_ijkl(1,1), n_i, N_int ) + call bitstring_to_list( mask_ijkl(1,2), list_ijkl(1,2), n_j, N_int ) + call bitstring_to_list( mask_ijkl(1,3), list_ijkl(1,3), n_k, N_int ) + call bitstring_to_list( mask_ijkl(1,4), list_ijkl(1,4), n_l, N_int ) + print*,'n_i,n_j = ',n_i,n_j + print*,'n_k,n_l = ',n_k,n_l + do i =1, n_i + print*,list_ijkl(i,1), list_ijkl(i,2) + enddo + deallocate(list_ijkl) + + +end diff --git a/plugins/core_integrals/AO_Basis b/plugins/core_integrals/AO_Basis deleted file mode 120000 index bee933ae..00000000 --- a/plugins/core_integrals/AO_Basis +++ /dev/null @@ -1 +0,0 @@ -/home/giner/qp_bis/quantum_package/src/AO_Basis \ No newline at end of file diff --git a/plugins/core_integrals/Bitmask b/plugins/core_integrals/Bitmask deleted file mode 120000 index 2f693e1c..00000000 --- a/plugins/core_integrals/Bitmask +++ /dev/null @@ -1 +0,0 @@ -/home/giner/qp_bis/quantum_package/src/Bitmask \ No newline at end of file diff --git a/plugins/core_integrals/Electrons b/plugins/core_integrals/Electrons deleted file mode 120000 index c29d5ab7..00000000 --- a/plugins/core_integrals/Electrons +++ /dev/null @@ -1 +0,0 @@ -/home/giner/qp_bis/quantum_package/src/Electrons \ No newline at end of file diff --git a/plugins/core_integrals/Ezfio_files b/plugins/core_integrals/Ezfio_files deleted file mode 120000 index 07172795..00000000 --- a/plugins/core_integrals/Ezfio_files +++ /dev/null @@ -1 +0,0 @@ -/home/giner/qp_bis/quantum_package/src/Ezfio_files \ No newline at end of file diff --git a/plugins/core_integrals/Integrals_Bielec b/plugins/core_integrals/Integrals_Bielec deleted file mode 120000 index fe90d314..00000000 --- a/plugins/core_integrals/Integrals_Bielec +++ /dev/null @@ -1 +0,0 @@ -/home/giner/qp_bis/quantum_package/src/Integrals_Bielec \ No newline at end of file diff --git a/plugins/core_integrals/Integrals_Monoelec b/plugins/core_integrals/Integrals_Monoelec deleted file mode 120000 index ef392ff9..00000000 --- a/plugins/core_integrals/Integrals_Monoelec +++ /dev/null @@ -1 +0,0 @@ -/home/giner/qp_bis/quantum_package/src/Integrals_Monoelec \ No newline at end of file diff --git a/plugins/core_integrals/MO_Basis b/plugins/core_integrals/MO_Basis deleted file mode 120000 index 55d01a67..00000000 --- a/plugins/core_integrals/MO_Basis +++ /dev/null @@ -1 +0,0 @@ -/home/giner/qp_bis/quantum_package/src/MO_Basis \ No newline at end of file diff --git a/plugins/core_integrals/Makefile b/plugins/core_integrals/Makefile deleted file mode 100644 index bfcd2493..00000000 --- a/plugins/core_integrals/Makefile +++ /dev/null @@ -1,16 +0,0 @@ -IRPF90 = irpf90 #-a -d -FC = gfortran -FCFLAGS= -O2 -ffree-line-length-none -I . -NINJA = ninja -AR = ar -RANLIB = ranlib - -SRC= -OBJ= -LIB= - -include irpf90.make -export - -irpf90.make: $(filter-out IRPF90_temp/%, $(wildcard */*.irp.f)) $(wildcard *.irp.f) $(wildcard *.inc.f) Makefile - $(IRPF90) diff --git a/plugins/core_integrals/Nuclei b/plugins/core_integrals/Nuclei deleted file mode 120000 index a21d375b..00000000 --- a/plugins/core_integrals/Nuclei +++ /dev/null @@ -1 +0,0 @@ -/home/giner/qp_bis/quantum_package/src/Nuclei \ No newline at end of file diff --git a/plugins/core_integrals/Pseudo b/plugins/core_integrals/Pseudo deleted file mode 120000 index 69623db7..00000000 --- a/plugins/core_integrals/Pseudo +++ /dev/null @@ -1 +0,0 @@ -/home/giner/qp_bis/quantum_package/src/Pseudo \ No newline at end of file diff --git a/plugins/core_integrals/Utils b/plugins/core_integrals/Utils deleted file mode 120000 index 22bcf3f5..00000000 --- a/plugins/core_integrals/Utils +++ /dev/null @@ -1 +0,0 @@ -/home/giner/qp_bis/quantum_package/src/Utils \ No newline at end of file diff --git a/plugins/core_integrals/ZMQ b/plugins/core_integrals/ZMQ deleted file mode 120000 index cbbfca4c..00000000 --- a/plugins/core_integrals/ZMQ +++ /dev/null @@ -1 +0,0 @@ -/home/giner/qp_bis/quantum_package/src/ZMQ \ No newline at end of file diff --git a/plugins/core_integrals/core_integrals.main b/plugins/core_integrals/core_integrals.main deleted file mode 100755 index de1348c4..00000000 Binary files a/plugins/core_integrals/core_integrals.main and /dev/null differ