From 3e12b2f359e7f56e1acf42fca6c63a203d66d4da Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 17 Apr 2017 12:41:00 +0200 Subject: [PATCH] Handling of two different mo and ao integrals map --- plugins/DFT_Utils/grid_density.irp.f | 3 +- plugins/Integrals_erf/EZFIO.cfg | 42 - plugins/Integrals_erf/NEEDED_CHILDREN_MODULES | 2 +- plugins/Integrals_erf/ao_bi_integrals.irp.f | 1121 ------------- .../Integrals_erf/ao_bi_integrals_erf.irp.f | 2 +- .../ao_bielec_integrals_in_map_slave.irp.f | 225 --- plugins/Integrals_erf/gauss_legendre.irp.f | 66 - plugins/Integrals_erf/integrals_3_index.irp.f | 22 - .../Integrals_erf/integrals_3_index_erf.irp.f | 22 + plugins/Integrals_erf/map_integrals.irp.f | 717 --------- ...ng_range.irp.f => map_integrals_erf.irp.f} | 8 +- plugins/Integrals_erf/mo_bi_integrals.irp.f | 1396 ----------------- .../Integrals_erf/mo_bi_integrals_erf.irp.f | 4 +- .../Integrals_erf/providers_ao_standard.irp.f | 116 -- ...{qp_ao_ints.irp.f => qp_ao_erf_ints.irp.f} | 10 +- plugins/Integrals_erf/read_write.irp.f | 48 - .../test_bitmasks_integrals.irp.f | 33 - .../NEEDED_CHILDREN_MODULES | 1 + plugins/Integrals_restart_DFT/README.rst | 12 + .../short_range_coulomb.irp.f | 79 + .../write_integrals_restart_dft.irp.f | 18 + .../Slater_rules_DFT/NEEDED_CHILDREN_MODULES | 1 + plugins/Slater_rules_DFT/README.rst | 12 + .../Slater_rules_DFT.main.irp.f | 38 + plugins/Slater_rules_DFT/energy.irp.f | 7 + .../Slater_rules_DFT/slater_rules_erf.irp.f | 445 ++++++ src/Integrals_Monoelec/EZFIO.cfg | 8 + src/Integrals_Monoelec/mo_mono_ints.irp.f | 5 +- src/Integrals_Monoelec/read_write.irp.f | 9 +- 29 files changed, 668 insertions(+), 3804 deletions(-) delete mode 100644 plugins/Integrals_erf/ao_bi_integrals.irp.f delete mode 100644 plugins/Integrals_erf/ao_bielec_integrals_in_map_slave.irp.f delete mode 100644 plugins/Integrals_erf/gauss_legendre.irp.f delete mode 100644 plugins/Integrals_erf/integrals_3_index.irp.f create mode 100644 plugins/Integrals_erf/integrals_3_index_erf.irp.f delete mode 100644 plugins/Integrals_erf/map_integrals.irp.f rename plugins/Integrals_erf/{map_integrals_long_range.irp.f => map_integrals_erf.irp.f} (98%) delete mode 100644 plugins/Integrals_erf/mo_bi_integrals.irp.f delete mode 100644 plugins/Integrals_erf/providers_ao_standard.irp.f rename plugins/Integrals_erf/{qp_ao_ints.irp.f => qp_ao_erf_ints.irp.f} (65%) delete mode 100644 plugins/Integrals_erf/test_bitmasks_integrals.irp.f create mode 100644 plugins/Integrals_restart_DFT/NEEDED_CHILDREN_MODULES create mode 100644 plugins/Integrals_restart_DFT/README.rst create mode 100644 plugins/Integrals_restart_DFT/short_range_coulomb.irp.f create mode 100644 plugins/Integrals_restart_DFT/write_integrals_restart_dft.irp.f create mode 100644 plugins/Slater_rules_DFT/NEEDED_CHILDREN_MODULES create mode 100644 plugins/Slater_rules_DFT/README.rst create mode 100644 plugins/Slater_rules_DFT/Slater_rules_DFT.main.irp.f create mode 100644 plugins/Slater_rules_DFT/energy.irp.f create mode 100644 plugins/Slater_rules_DFT/slater_rules_erf.irp.f diff --git a/plugins/DFT_Utils/grid_density.irp.f b/plugins/DFT_Utils/grid_density.irp.f index ad5e0d51..7c9d2c05 100644 --- a/plugins/DFT_Utils/grid_density.irp.f +++ b/plugins/DFT_Utils/grid_density.irp.f @@ -5,7 +5,7 @@ BEGIN_PROVIDER [integer, n_points_radial_grid] implicit none - n_points_radial_grid = 10000 + n_points_radial_grid = 100 END_PROVIDER @@ -188,6 +188,7 @@ END_PROVIDER call give_all_mos_at_r(r,mos_array) do m = 1, mo_tot_num do i = 1, mo_tot_num + if(dabs(one_body_dm_mo_alpha(i,m,i_state)).lt.1.d-10)cycle contrib = mos_array(i) * mos_array(m) one_body_dm_mo_alpha_at_grid_points(l,k,j,i_state) += one_body_dm_mo_alpha(i,m,i_state) * contrib one_body_dm_mo_beta_at_grid_points(l,k,j,i_state) += one_body_dm_mo_beta(i,m,i_state) * contrib diff --git a/plugins/Integrals_erf/EZFIO.cfg b/plugins/Integrals_erf/EZFIO.cfg index ca798f68..916bcd34 100644 --- a/plugins/Integrals_erf/EZFIO.cfg +++ b/plugins/Integrals_erf/EZFIO.cfg @@ -1,37 +1,3 @@ -[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_erf] type: Disk_access doc: Read/Write AO integrals with the long range interaction from/to disk [ Write | Read | None ] @@ -45,14 +11,6 @@ doc: Read/Write MO integrals with the long range interaction from/to disk [ Writ 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 diff --git a/plugins/Integrals_erf/NEEDED_CHILDREN_MODULES b/plugins/Integrals_erf/NEEDED_CHILDREN_MODULES index 152711f3..8361b2eb 100644 --- a/plugins/Integrals_erf/NEEDED_CHILDREN_MODULES +++ b/plugins/Integrals_erf/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Pseudo Bitmask ZMQ +Pseudo Bitmask ZMQ Integrals_Bielec diff --git a/plugins/Integrals_erf/ao_bi_integrals.irp.f b/plugins/Integrals_erf/ao_bi_integrals.irp.f deleted file mode 100644 index b4c9adb5..00000000 --- a/plugins/Integrals_erf/ao_bi_integrals.irp.f +++ /dev/null @@ -1,1121 +0,0 @@ -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 index 6616883f..2b4b2fad 100644 --- a/plugins/Integrals_erf/ao_bi_integrals_erf.irp.f +++ b/plugins/Integrals_erf/ao_bi_integrals_erf.irp.f @@ -294,7 +294,7 @@ subroutine compute_ao_bielec_integrals_erf(j,k,l,sze,buffer_value) buffer_value = 0._integral_kind return endif - if (ao_bielec_integral_schwartz(j,l) < thresh ) then + if (ao_bielec_integral_erf_schwartz(j,l) < thresh ) then buffer_value = 0._integral_kind return endif 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 deleted file mode 100644 index 38c78388..00000000 --- a/plugins/Integrals_erf/ao_bielec_integrals_in_map_slave.irp.f +++ /dev/null @@ -1,225 +0,0 @@ -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 deleted file mode 100644 index e62febeb..00000000 --- a/plugins/Integrals_erf/gauss_legendre.irp.f +++ /dev/null @@ -1,66 +0,0 @@ -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 deleted file mode 100644 index 41037b34..00000000 --- a/plugins/Integrals_erf/integrals_3_index.irp.f +++ /dev/null @@ -1,22 +0,0 @@ - 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/integrals_3_index_erf.irp.f b/plugins/Integrals_erf/integrals_3_index_erf.irp.f new file mode 100644 index 00000000..d9b1e9f7 --- /dev/null +++ b/plugins/Integrals_erf/integrals_3_index_erf.irp.f @@ -0,0 +1,22 @@ + BEGIN_PROVIDER [double precision, big_array_coulomb_integrals_erf, (mo_tot_num_align,mo_tot_num, mo_tot_num)] +&BEGIN_PROVIDER [double precision, big_array_exchange_integrals_erf,(mo_tot_num_align,mo_tot_num, mo_tot_num)] + implicit none + integer :: i,j,k,l + double precision :: get_mo_bielec_integral_erf + 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_erf(i,j,k,l,mo_integrals_erf_map) + big_array_coulomb_integrals_erf(j,i,k) = integral + l = j + integral = get_mo_bielec_integral_erf(i,j,l,k,mo_integrals_erf_map) + big_array_exchange_integrals_erf(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 deleted file mode 100644 index 82b89f22..00000000 --- a/plugins/Integrals_erf/map_integrals.irp.f +++ /dev/null @@ -1,717 +0,0 @@ -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_erf.irp.f similarity index 98% rename from plugins/Integrals_erf/map_integrals_long_range.irp.f rename to plugins/Integrals_erf/map_integrals_erf.irp.f index 8f24692a..ecf72282 100644 --- a/plugins/Integrals_erf/map_integrals_long_range.irp.f +++ b/plugins/Integrals_erf/map_integrals_erf.irp.f @@ -156,7 +156,7 @@ subroutine get_ao_bielec_integrals_erf_non_zero(j,k,l,sze,out_val,out_val_index, 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 + if (ao_bielec_integral_erf_schwartz(i,k)*ao_bielec_integral_erf_schwartz(j,l) < thresh) then cycle endif call bielec_integrals_index(i,j,k,l,hash) @@ -395,7 +395,7 @@ BEGIN_PROVIDER [ double precision, mo_integrals_erf_cache, (0:64*64*64*64) ] integer :: ii integer(key_kind) :: idx real(integral_kind) :: integral - FREE ao_integrals_cache + FREE ao_integrals_erf_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 @@ -478,7 +478,7 @@ subroutine get_mo_bielec_integrals_erf(j,k,l,sze,out_val,map) integer :: i integer(key_kind) :: hash(sze) real(integral_kind) :: tmp_val(sze) - PROVIDE mo_bielec_integrals_in_map + PROVIDE mo_bielec_integrals_erf_in_map do i=1,sze !DIR$ FORCEINLINE @@ -564,7 +564,7 @@ subroutine get_mo_bielec_integrals_erf_coulomb_ii(k,l,sze,out_val,map) integer :: i integer(key_kind) :: hash(sze) real(integral_kind) :: tmp_val(sze) - PROVIDE mo_bielec_integrals_in_map + PROVIDE mo_bielec_integrals_erf_in_map integer :: kk do i=1,sze diff --git a/plugins/Integrals_erf/mo_bi_integrals.irp.f b/plugins/Integrals_erf/mo_bi_integrals.irp.f deleted file mode 100644 index 9cca2e12..00000000 --- a/plugins/Integrals_erf/mo_bi_integrals.irp.f +++ /dev/null @@ -1,1396 +0,0 @@ -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_erf_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 -! 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 - 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 index ecb55757..b0c954c1 100644 --- a/plugins/Integrals_erf/mo_bi_integrals_erf.irp.f +++ b/plugins/Integrals_erf/mo_bi_integrals_erf.irp.f @@ -587,8 +587,8 @@ END_PROVIDER 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(i,j) = get_mo_bielec_integral_erf(i,j,i,j,mo_integrals_erf_map) + mo_bielec_integral_erf_jj_exchange(i,j) = get_mo_bielec_integral_erf(i,j,j,i,mo_integrals_erf_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 diff --git a/plugins/Integrals_erf/providers_ao_standard.irp.f b/plugins/Integrals_erf/providers_ao_standard.irp.f deleted file mode 100644 index 973c0e58..00000000 --- a/plugins/Integrals_erf/providers_ao_standard.irp.f +++ /dev/null @@ -1,116 +0,0 @@ - -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_erf_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_erf_ints.irp.f similarity index 65% rename from plugins/Integrals_erf/qp_ao_ints.irp.f rename to plugins/Integrals_erf/qp_ao_erf_ints.irp.f index 93f62a7d..df6d8d16 100644 --- a/plugins/Integrals_erf/qp_ao_ints.irp.f +++ b/plugins/Integrals_erf/qp_ao_erf_ints.irp.f @@ -2,7 +2,7 @@ program qp_ao_ints use omp_lib implicit none BEGIN_DOC -! Increments a running calculation to compute AO integrals +! Increments a running calculation to compute AO integral_erfs END_DOC integer :: i @@ -11,18 +11,18 @@ program qp_ao_ints zmq_context = f77_zmq_ctx_new () ! Set the state of the ZMQ - zmq_state = 'ao_integrals' + zmq_state = 'ao_integral_erfs' ! Provide everything needed - double precision :: integral, ao_bielec_integral - integral = ao_bielec_integral(1,1,1,1) + double precision :: integral_erf, ao_bielec_integral_erf + integral_erf = ao_bielec_integral_erf(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) + call ao_bielec_integrals_erf_in_map_slave_tcp(i) !$OMP END PARALLEL call wait_for_state(zmq_state,state) enddo diff --git a/plugins/Integrals_erf/read_write.irp.f b/plugins/Integrals_erf/read_write.irp.f index e9fc0f91..12bbf0bc 100644 --- a/plugins/Integrals_erf/read_write.irp.f +++ b/plugins/Integrals_erf/read_write.irp.f @@ -1,51 +1,3 @@ -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 - BEGIN_PROVIDER [ logical, read_ao_integrals_erf ] &BEGIN_PROVIDER [ logical, read_mo_integrals_erf ] &BEGIN_PROVIDER [ logical, write_ao_integrals_erf ] diff --git a/plugins/Integrals_erf/test_bitmasks_integrals.irp.f b/plugins/Integrals_erf/test_bitmasks_integrals.irp.f deleted file mode 100644 index d0eba3f7..00000000 --- a/plugins/Integrals_erf/test_bitmasks_integrals.irp.f +++ /dev/null @@ -1,33 +0,0 @@ -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/Integrals_restart_DFT/NEEDED_CHILDREN_MODULES b/plugins/Integrals_restart_DFT/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..08317b5e --- /dev/null +++ b/plugins/Integrals_restart_DFT/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Integrals_Monoelec Integrals_erf Determinants DFT_Utils diff --git a/plugins/Integrals_restart_DFT/README.rst b/plugins/Integrals_restart_DFT/README.rst new file mode 100644 index 00000000..589e0a00 --- /dev/null +++ b/plugins/Integrals_restart_DFT/README.rst @@ -0,0 +1,12 @@ +============== +core_integrals +============== + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/plugins/Integrals_restart_DFT/short_range_coulomb.irp.f b/plugins/Integrals_restart_DFT/short_range_coulomb.irp.f new file mode 100644 index 00000000..aeb2589c --- /dev/null +++ b/plugins/Integrals_restart_DFT/short_range_coulomb.irp.f @@ -0,0 +1,79 @@ +BEGIN_PROVIDER [double precision, density_matrix_read, (mo_tot_num, mo_tot_num)] + implicit none + integer :: i,j,k,l + logical :: exists + call ezfio_has_determinants_density_matrix_mo_disk(exists) + if(exists)then + print*, 'reading the density matrix from input' + call ezfio_get_determinants_density_matrix_mo_disk(exists) + print*, 'reading done' + else + print*, 'no density matrix found in EZFIO file ...' + print*, 'stopping ..' + stop + endif + +END_PROVIDER + + +BEGIN_PROVIDER [double precision, effective_short_range_operator, (mo_tot_num,mo_tot_num)] + implicit none + integer :: i,j,k,l,m,n + double precision :: get_mo_bielec_integral,get_mo_bielec_integral_erf + double precision :: integral, integral_erf + effective_short_range_operator = 0.d0 + do i = 1, mo_tot_num + do j = 1, mo_tot_num + if(dabs(one_body_dm_mo(i,j)).le.1.d-10)cycle + do k = 1, mo_tot_num + do l = 1, mo_tot_num + integral = get_mo_bielec_integral(i,k,j,l,mo_integrals_map) +! integral_erf = get_mo_bielec_integral_erf(i,k,j,l,mo_integrals_erf_map) + effective_short_range_operator(l,k) += one_body_dm_mo(i,j) * integral + enddo + enddo + enddo + enddo +END_PROVIDER + + +BEGIN_PROVIDER [double precision, effective_one_e_potential, (mo_tot_num_align, mo_tot_num,N_states)] + implicit none + integer :: i,j,i_state + effective_one_e_potential = 0.d0 + do i_state = 1, N_states + do i = 1, mo_tot_num + do j = 1, mo_tot_num + effective_one_e_potential(i,j,i_state) = effective_short_range_operator(i,j) + mo_nucl_elec_integral(i,j) + mo_kinetic_integral(i,j) & + + 0.5d0 * (lda_ex_potential_alpha_ao(i,j,i_state) + lda_ex_potential_beta_ao(i,j,i_state)) + enddo + enddo + enddo + +END_PROVIDER + +subroutine save_one_e_effective_potential + implicit none + double precision, allocatable :: tmp(:,:) + allocate(tmp(size(effective_one_e_potential,1),size(effective_one_e_potential,2))) + integer :: i,j + do i = 1, mo_tot_num + do j = 1, mo_tot_num + tmp(i,j) = effective_one_e_potential(i,j,1) + enddo + enddo + call write_one_e_integrals('mo_one_integral', tmp, & + size(tmp,1), size(tmp,2)) + call ezfio_set_integrals_monoelec_disk_access_only_mo_one_integrals("Read") + deallocate(tmp) + +end + +subroutine save_erf_bi_elec_integrals + implicit none + integer :: i,j,k,l + PROVIDE mo_bielec_integrals_erf_in_map + call ezfio_set_work_empty(.False.) + call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_erf_map) + call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read") +end diff --git a/plugins/Integrals_restart_DFT/write_integrals_restart_dft.irp.f b/plugins/Integrals_restart_DFT/write_integrals_restart_dft.irp.f new file mode 100644 index 00000000..d89b965d --- /dev/null +++ b/plugins/Integrals_restart_DFT/write_integrals_restart_dft.irp.f @@ -0,0 +1,18 @@ +program write_integrals + implicit none + read_wf = .true. + touch read_wf + disk_access_only_mo_one_integrals = "None" + touch disk_access_only_mo_one_integrals + disk_access_mo_integrals = "None" + touch disk_access_mo_integrals + call routine + +end + +subroutine routine + implicit none + call save_one_e_effective_potential + call save_erf_bi_elec_integrals + +end diff --git a/plugins/Slater_rules_DFT/NEEDED_CHILDREN_MODULES b/plugins/Slater_rules_DFT/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..994f4bf6 --- /dev/null +++ b/plugins/Slater_rules_DFT/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Determinants Integrals_restart_DFT Davidson diff --git a/plugins/Slater_rules_DFT/README.rst b/plugins/Slater_rules_DFT/README.rst new file mode 100644 index 00000000..f492095e --- /dev/null +++ b/plugins/Slater_rules_DFT/README.rst @@ -0,0 +1,12 @@ +================ +Slater_rules_DFT +================ + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/plugins/Slater_rules_DFT/Slater_rules_DFT.main.irp.f b/plugins/Slater_rules_DFT/Slater_rules_DFT.main.irp.f new file mode 100644 index 00000000..3d99e376 --- /dev/null +++ b/plugins/Slater_rules_DFT/Slater_rules_DFT.main.irp.f @@ -0,0 +1,38 @@ +program Slater_rules_DFT + implicit none + BEGIN_DOC +! TODO + END_DOC + print *, ' _/ ' + print *, ' -:\_?, _Jm####La ' + print *, 'J"(:" > _]#AZ#Z#UUZ##, ' + print *, '_,::./ %(|i%12XmX1*1XL _?, ' + print *, ' \..\ _\(vmWQwodY+ia%lnL _",/ ( ' + print *, ' .:< ]J=mQD?WXn|,)nr" ' + print *, ' 4XZ#Xov1v}=)vnXAX1nnv;1n" ' + print *, ' ]XX#ZXoovvvivnnnlvvo2*i7 ' + print *, ' "23Z#1S2oo2XXSnnnoSo2>v" ' + print *, ' miX#L -~`""!!1}oSoe|i7 ' + print *, ' 4cn#m, v221=|v[ ' + print *, ' ]hI3Zma,;..__wXSe=+vo ' + print *, ' ]Zov*XSUXXZXZXSe||vo2 ' + print *, ' ]Z#>=|< ' + print *, ' -ziiiii||||||+||==+> ' + print *, ' -%|+++||=|=+|=|==/ ' + print *, ' -a>====+|====-:- ' + print *, ' "~,- -- /- ' + print *, ' -. )> ' + print *, ' .~ +- ' + print *, ' . .... : . ' + print *, ' -------~ ' + print *, '' +end diff --git a/plugins/Slater_rules_DFT/energy.irp.f b/plugins/Slater_rules_DFT/energy.irp.f new file mode 100644 index 00000000..7734d73e --- /dev/null +++ b/plugins/Slater_rules_DFT/energy.irp.f @@ -0,0 +1,7 @@ +! BEGIN_PROVIDER [double precision, energy_total] +!&BEGIN_PROVIDER [double precision, energy_one_electron] +!&BEGIN_PROVIDER [double precision, energy_hartree] +!&BEGIN_PROVIDER [double precision, energy] +! implicit none +! +!END_PROVIDER diff --git a/plugins/Slater_rules_DFT/slater_rules_erf.irp.f b/plugins/Slater_rules_DFT/slater_rules_erf.irp.f new file mode 100644 index 00000000..64d5d217 --- /dev/null +++ b/plugins/Slater_rules_DFT/slater_rules_erf.irp.f @@ -0,0 +1,445 @@ + +subroutine i_H_j_erf(key_i,key_j,Nint,hij) + use bitmasks + implicit none + BEGIN_DOC + ! Returns where i and j are determinants + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hij + + integer :: exc(0:2,2,2) + integer :: degree + double precision :: get_mo_bielec_integral_erf + integer :: m,n,p,q + integer :: i,j,k + integer :: occ(Nint*bit_kind_size,2) + double precision :: diag_H_mat_elem_erf, phase,phase_2 + integer :: n_occ_ab(2) + PROVIDE mo_bielec_integrals_erf_in_map mo_integrals_erf_map big_array_exchange_integrals_erf + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num) + ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num) + + hij = 0.d0 + !DIR$ FORCEINLINE + call get_excitation_degree(key_i,key_j,degree,Nint) + integer :: spin + select case (degree) + case (2) + call get_double_excitation(key_i,key_j,exc,phase,Nint) + if (exc(0,1,1) == 1) then + ! Mono alpha, mono beta + if(exc(1,1,1) == exc(1,2,2) )then + hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1)) + else if (exc(1,2,1) ==exc(1,1,2))then + hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2)) + else + hij = phase*get_mo_bielec_integral_erf( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_erf_map) + endif + else if (exc(0,1,1) == 2) then + ! Double alpha + hij = phase*(get_mo_bielec_integral_erf( & + exc(1,1,1), & + exc(2,1,1), & + exc(1,2,1), & + exc(2,2,1) ,mo_integrals_erf_map) - & + get_mo_bielec_integral_erf( & + exc(1,1,1), & + exc(2,1,1), & + exc(2,2,1), & + exc(1,2,1) ,mo_integrals_erf_map) ) + else if (exc(0,1,2) == 2) then + ! Double beta + hij = phase*(get_mo_bielec_integral_erf( & + exc(1,1,2), & + exc(2,1,2), & + exc(1,2,2), & + exc(2,2,2) ,mo_integrals_erf_map) - & + get_mo_bielec_integral_erf( & + exc(1,1,2), & + exc(2,1,2), & + exc(2,2,2), & + exc(1,2,2) ,mo_integrals_erf_map) ) + endif + case (1) + call get_mono_excitation(key_i,key_j,exc,phase,Nint) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) + if (exc(0,1,1) == 1) then + ! Mono alpha + m = exc(1,1,1) + p = exc(1,2,1) + spin = 1 + do i = 1, n_occ_ab(1) + hij += -big_array_exchange_integrals_erf(occ(i,1),m,p) + big_array_coulomb_integrals_erf(occ(i,1),m,p) + enddo + do i = 1, n_occ_ab(2) + hij += big_array_coulomb_integrals_erf(occ(i,2),m,p) + enddo + else + ! Mono beta + m = exc(1,1,2) + p = exc(1,2,2) + spin = 2 + do i = 1, n_occ_ab(2) + hij += -big_array_exchange_integrals_erf(occ(i,2),m,p) + big_array_coulomb_integrals_erf(occ(i,2),m,p) + enddo + do i = 1, n_occ_ab(1) + hij += big_array_coulomb_integrals_erf(occ(i,1),m,p) + enddo + endif + hij = hij + mo_nucl_elec_integral(m,p) + mo_kinetic_integral(m,p) + hij = hij * phase + case (0) + hij = diag_H_mat_elem_erf(key_i,Nint) + end select +end + +double precision function diag_H_mat_elem_erf(key_i,Nint) + implicit none + integer(bit_kind), intent(in) :: key_i(N_int,2) + integer, intent(in) :: Nint + integer :: i,j + integer :: occ(Nint*bit_kind_size,2) + integer :: n_occ_ab(2) + call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) + diag_H_mat_elem_erf = 0.d0 + ! alpha - alpha + do i = 1, n_occ_ab(1) + diag_H_mat_elem_erf += mo_nucl_elec_integral(occ(i,1),mo_nucl_elec_integral(i,1)) + do j = i+1, n_occ_ab(1) + diag_H_mat_elem_erf += mo_bielec_integral_erf_jj_anti(occ(i,1),occ(j,1)) + enddo + enddo + + ! beta - beta + do i = 1, n_occ_ab(2) + diag_H_mat_elem_erf += mo_nucl_elec_integral(occ(i,2),mo_nucl_elec_integral(i,2)) + do j = i+1, n_occ_ab(2) + diag_H_mat_elem_erf += mo_bielec_integral_erf_jj_anti(occ(i,2),occ(j,2)) + enddo + enddo + + ! alpha - beta + do i = 1, n_occ_ab(1) + do j = 1, n_occ_ab(2) + diag_H_mat_elem_erf += mo_bielec_integral_erf_jj(occ(i,1),occ(j,2)) + enddo + enddo + +end + + + +subroutine i_H_j_erf_and_short_coulomb(key_i,key_j,Nint,hij) + use bitmasks + implicit none + BEGIN_DOC + ! Returns where i and j are determinants + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hij + + integer :: exc(0:2,2,2) + integer :: degree + double precision :: get_mo_bielec_integral_erf + integer :: m,n,p,q + integer :: i,j,k + integer :: occ(Nint*bit_kind_size,2) + double precision :: diag_H_mat_elem_erf, phase,phase_2 + integer :: n_occ_ab(2) + PROVIDE mo_bielec_integrals_erf_in_map mo_integrals_erf_map big_array_exchange_integrals_erf + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num) + ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num) + + hij = 0.d0 + !DIR$ FORCEINLINE + call get_excitation_degree(key_i,key_j,degree,Nint) + integer :: spin + select case (degree) + case (2) + call get_double_excitation(key_i,key_j,exc,phase,Nint) + if (exc(0,1,1) == 1) then + ! Mono alpha, mono beta + if(exc(1,1,1) == exc(1,2,2) )then + hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1)) + else if (exc(1,2,1) ==exc(1,1,2))then + hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2)) + else + hij = phase*get_mo_bielec_integral_erf( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_erf_map) + endif + else if (exc(0,1,1) == 2) then + ! Double alpha + hij = phase*(get_mo_bielec_integral_erf( & + exc(1,1,1), & + exc(2,1,1), & + exc(1,2,1), & + exc(2,2,1) ,mo_integrals_erf_map) - & + get_mo_bielec_integral_erf( & + exc(1,1,1), & + exc(2,1,1), & + exc(2,2,1), & + exc(1,2,1) ,mo_integrals_erf_map) ) + else if (exc(0,1,2) == 2) then + ! Double beta + hij = phase*(get_mo_bielec_integral_erf( & + exc(1,1,2), & + exc(2,1,2), & + exc(1,2,2), & + exc(2,2,2) ,mo_integrals_erf_map) - & + get_mo_bielec_integral_erf( & + exc(1,1,2), & + exc(2,1,2), & + exc(2,2,2), & + exc(1,2,2) ,mo_integrals_erf_map) ) + endif + case (1) + call get_mono_excitation(key_i,key_j,exc,phase,Nint) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) + if (exc(0,1,1) == 1) then + ! Mono alpha + m = exc(1,1,1) + p = exc(1,2,1) + spin = 1 + do i = 1, n_occ_ab(1) + hij += -big_array_exchange_integrals_erf(occ(i,1),m,p) + big_array_coulomb_integrals_erf(occ(i,1),m,p) + enddo + do i = 1, n_occ_ab(2) + hij += big_array_coulomb_integrals_erf(occ(i,2),m,p) + enddo + else + ! Mono beta + m = exc(1,1,2) + p = exc(1,2,2) + spin = 2 + do i = 1, n_occ_ab(2) + hij += -big_array_exchange_integrals_erf(occ(i,2),m,p) + big_array_coulomb_integrals_erf(occ(i,2),m,p) + enddo + do i = 1, n_occ_ab(1) + hij += big_array_coulomb_integrals_erf(occ(i,1),m,p) + enddo + endif + hij = hij + mo_nucl_elec_integral(m,p) + mo_kinetic_integral(m,p) + effective_short_range_operator(m,p) + hij = hij * phase + case (0) + hij = diag_H_mat_elem_erf(key_i,Nint) + end select +end + +double precision function diag_H_mat_elem_erf_and_short_coulomb(key_i,Nint) + implicit none + integer(bit_kind), intent(in) :: key_i(N_int,2) + integer, intent(in) :: Nint + integer :: i,j + integer :: occ(Nint*bit_kind_size,2) + integer :: n_occ_ab(2) + + call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) + diag_H_mat_elem_erf_and_short_coulomb = 0.d0 + ! alpha - alpha + do i = 1, n_occ_ab(1) + diag_H_mat_elem_erf_and_short_coulomb += mo_nucl_elec_integral(occ(i,1),mo_nucl_elec_integral(i,1)) + mo_kinetic_integral(occ(i,1),mo_nucl_elec_integral(i,1)) & + + effective_short_range_operator(occ(i,1),occ(i,1)) + do j = i+1, n_occ_ab(1) + diag_H_mat_elem_erf_and_short_coulomb += mo_bielec_integral_erf_jj_anti(occ(i,1),occ(j,1)) + enddo + enddo + + ! beta - beta + do i = 1, n_occ_ab(2) + diag_H_mat_elem_erf_and_short_coulomb += mo_nucl_elec_integral(occ(i,2),mo_nucl_elec_integral(i,2)) + mo_kinetic_integral(occ(i,2),mo_nucl_elec_integral(i,2)) & + + effective_short_range_operator(occ(i,2),occ(i,2)) + do j = i+1, n_occ_ab(2) + diag_H_mat_elem_erf_and_short_coulomb += mo_bielec_integral_erf_jj_anti(occ(i,2),occ(j,2)) + enddo + enddo + + ! alpha - beta + do i = 1, n_occ_ab(1) + do j = 1, n_occ_ab(2) + diag_H_mat_elem_erf_and_short_coulomb += mo_bielec_integral_erf_jj(occ(i,1),occ(j,2)) + enddo + enddo + +end + + +subroutine i_H_j_erf_component(key_i,key_j,Nint,hij_core,hij_hartree,hij_erf,hij_total) + use bitmasks + implicit none + BEGIN_DOC + ! Returns where i and j are determinants + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hij_core + double precision, intent(out) :: hij_hartree + double precision, intent(out) :: hij_erf + double precision, intent(out) :: hij_total + + integer :: exc(0:2,2,2) + integer :: degree + double precision :: get_mo_bielec_integral_erf + integer :: m,n,p,q + integer :: i,j,k + integer :: occ(Nint*bit_kind_size,2) + double precision :: diag_H_mat_elem_erf, phase,phase_2 + integer :: n_occ_ab(2) + PROVIDE mo_bielec_integrals_erf_in_map mo_integrals_erf_map big_array_exchange_integrals_erf + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num) + ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num) + + hij_core = 0.d0 + hij_hartree = 0.d0 + hij_erf = 0.d0 + + !DIR$ FORCEINLINE + call get_excitation_degree(key_i,key_j,degree,Nint) + integer :: spin + select case (degree) + case (2) + call get_double_excitation(key_i,key_j,exc,phase,Nint) + if (exc(0,1,1) == 1) then + ! Mono alpha, mono beta + if(exc(1,1,1) == exc(1,2,2) )then + hij_erf = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1)) + else if (exc(1,2,1) ==exc(1,1,2))then + hij_erf = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2)) + else + hij_erf = phase*get_mo_bielec_integral_erf( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_erf_map) + endif + else if (exc(0,1,1) == 2) then + ! Double alpha + hij_erf = phase*(get_mo_bielec_integral_erf( & + exc(1,1,1), & + exc(2,1,1), & + exc(1,2,1), & + exc(2,2,1) ,mo_integrals_erf_map) - & + get_mo_bielec_integral_erf( & + exc(1,1,1), & + exc(2,1,1), & + exc(2,2,1), & + exc(1,2,1) ,mo_integrals_erf_map) ) + else if (exc(0,1,2) == 2) then + ! Double beta + hij_erf = phase*(get_mo_bielec_integral_erf( & + exc(1,1,2), & + exc(2,1,2), & + exc(1,2,2), & + exc(2,2,2) ,mo_integrals_erf_map) - & + get_mo_bielec_integral_erf( & + exc(1,1,2), & + exc(2,1,2), & + exc(2,2,2), & + exc(1,2,2) ,mo_integrals_erf_map) ) + endif + case (1) + call get_mono_excitation(key_i,key_j,exc,phase,Nint) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) + if (exc(0,1,1) == 1) then + ! Mono alpha + m = exc(1,1,1) + p = exc(1,2,1) + spin = 1 + do i = 1, n_occ_ab(1) + hij_erf += -big_array_exchange_integrals_erf(occ(i,1),m,p) + big_array_coulomb_integrals_erf(occ(i,1),m,p) + enddo + do i = 1, n_occ_ab(2) + hij_erf += big_array_coulomb_integrals_erf(occ(i,2),m,p) + enddo + else + ! Mono beta + m = exc(1,1,2) + p = exc(1,2,2) + spin = 2 + do i = 1, n_occ_ab(2) + hij_erf += -big_array_exchange_integrals_erf(occ(i,2),m,p) + big_array_coulomb_integrals_erf(occ(i,2),m,p) + enddo + do i = 1, n_occ_ab(1) + hij_erf += big_array_coulomb_integrals_erf(occ(i,1),m,p) + enddo + endif + hij_core = mo_nucl_elec_integral(m,p) + mo_kinetic_integral(m,p) + hij_hartree = effective_short_range_operator(m,p) + hij_total = (hij_erf + hij_core + hij_hartree) * phase + case (0) + call diag_H_mat_elem_erf_component(key_i,hij_core,hij_hartree,hij_erf,hij_total,Nint) + end select +end + +subroutine diag_H_mat_elem_erf_component(key_i,hij_core,hij_hartree,hij_erf,hij_total,Nint) + implicit none + integer(bit_kind), intent(in) :: key_i(N_int,2) + integer, intent(in) :: Nint + double precision, intent(out) :: hij_core + double precision, intent(out) :: hij_hartree + double precision, intent(out) :: hij_erf + double precision, intent(out) :: hij_total + integer :: i,j + integer :: occ(Nint*bit_kind_size,2) + integer :: n_occ_ab(2) + + call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) + hij_core = 0.d0 + hij_hartree = 0.d0 + hij_erf = 0.d0 + ! alpha - alpha + do i = 1, n_occ_ab(1) + hij_core += mo_nucl_elec_integral(occ(i,1),mo_nucl_elec_integral(i,1)) + mo_kinetic_integral(occ(i,1),mo_nucl_elec_integral(i,1)) + hij_hartree += effective_short_range_operator(occ(i,1),occ(i,1)) + do j = i+1, n_occ_ab(1) + hij_erf += mo_bielec_integral_erf_jj_anti(occ(i,1),occ(j,1)) + enddo + enddo + + ! beta - beta + do i = 1, n_occ_ab(2) + hij_core += mo_nucl_elec_integral(occ(i,2),mo_nucl_elec_integral(i,2)) + mo_kinetic_integral(occ(i,2),mo_nucl_elec_integral(i,2)) + hij_hartree += effective_short_range_operator(occ(i,2),occ(i,2)) + do j = i+1, n_occ_ab(2) + hij_erf += mo_bielec_integral_erf_jj_anti(occ(i,2),occ(j,2)) + enddo + enddo + + ! alpha - beta + do i = 1, n_occ_ab(1) + do j = 1, n_occ_ab(2) + hij_erf += mo_bielec_integral_erf_jj(occ(i,1),occ(j,2)) + enddo + enddo + hij_total = hij_erf + hij_hartree + hij_core + +end + + diff --git a/src/Integrals_Monoelec/EZFIO.cfg b/src/Integrals_Monoelec/EZFIO.cfg index 04e49ec1..c8a8eaef 100644 --- a/src/Integrals_Monoelec/EZFIO.cfg +++ b/src/Integrals_Monoelec/EZFIO.cfg @@ -4,6 +4,14 @@ doc: Read/Write MO one-electron integrals from/to disk [ Write | Read | None ] interface: ezfio,provider,ocaml default: None + +[disk_access_only_mo_one_integrals] +type: Disk_access +doc: Read/Write MO for only the total one-electron integrals which can be anything [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + + [disk_access_ao_one_integrals] type: Disk_access doc: Read/Write AO one-electron integrals from/to disk [ Write | Read | None ] diff --git a/src/Integrals_Monoelec/mo_mono_ints.irp.f b/src/Integrals_Monoelec/mo_mono_ints.irp.f index 20aae1fd..816dd277 100644 --- a/src/Integrals_Monoelec/mo_mono_ints.irp.f +++ b/src/Integrals_Monoelec/mo_mono_ints.irp.f @@ -6,7 +6,8 @@ BEGIN_PROVIDER [ double precision, mo_mono_elec_integral,(mo_tot_num_align,mo_to ! sum of the kinetic and nuclear electronic potential END_DOC print*,'Providing the mono electronic integrals' - if (read_mo_one_integrals) then + if (read_only_mo_one_integrals) then + print*, 'Reading the mono electronic integrals from disk' call read_one_e_integrals('mo_one_integral', mo_mono_elec_integral, & size(mo_mono_elec_integral,1), size(mo_mono_elec_integral,2)) print *, 'MO N-e integrals read from disk' @@ -14,7 +15,7 @@ BEGIN_PROVIDER [ double precision, mo_mono_elec_integral,(mo_tot_num_align,mo_to do j = 1, mo_tot_num do i = 1, mo_tot_num mo_mono_elec_integral(i,j) = mo_nucl_elec_integral(i,j) + & - mo_kinetic_integral(i,j) + mo_pseudo_integral(i,j) + mo_kinetic_integral(i,j) + mo_pseudo_integral(i,j) enddo enddo endif diff --git a/src/Integrals_Monoelec/read_write.irp.f b/src/Integrals_Monoelec/read_write.irp.f index 697bf356..0e758740 100644 --- a/src/Integrals_Monoelec/read_write.irp.f +++ b/src/Integrals_Monoelec/read_write.irp.f @@ -1,5 +1,6 @@ BEGIN_PROVIDER [ logical, read_ao_one_integrals ] &BEGIN_PROVIDER [ logical, read_mo_one_integrals ] +&BEGIN_PROVIDER [ logical, read_only_mo_one_integrals ] &BEGIN_PROVIDER [ logical, write_ao_one_integrals ] &BEGIN_PROVIDER [ logical, write_mo_one_integrals ] @@ -21,10 +22,14 @@ write_ao_one_integrals = .False. else - print *, 'bielec_integrals/disk_access_ao_integrals has a wrong type' + print *, 'monoelec_integrals/disk_access_ao_integrals has a wrong type' stop 1 endif + + if (disk_access_only_mo_one_integrals.EQ.'Read')then + read_only_mo_one_integrals = .True. + endif if (disk_access_mo_one_integrals.EQ.'Read') then read_mo_one_integrals = .True. @@ -39,7 +44,7 @@ write_mo_one_integrals = .False. else - print *, 'bielec_integrals/disk_access_mo_integrals has a wrong type' + print *, 'monoelec_integrals/disk_access_mo_integrals has a wrong type' stop 1 endif