Accelerated selection : no more O(N2)

plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f

@@ -375,7 +375,7 @@ BEGIN_PROVIDER [ integer, size_tbc ]
   BEGIN_DOC
 ! Size of the tbc array
   END_DOC
-  size_tbc = 2*N_det_generators + fragment_count*fragment_first
+  size_tbc = (comb_teeth+1)*N_det_generators + fragment_count*fragment_first
 END_PROVIDER
 
 subroutine get_carlo_workbatch(computed, comb, Ncomb, tbc)

plugins/Full_CI_ZMQ/selection.irp.f

@@ -275,6 +275,9 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
   use bitmasks
   use selection_types
   implicit none
+  BEGIN_DOC
+!            WARNING /!\ : It is assumed that the generators and selectors are psi_det_sorted
+  END_DOC
   
   integer, intent(in)            :: i_generator, subset
   integer(bit_kind), intent(in)  :: hole_mask(N_int,2), particle_mask(N_int,2)
@@ -332,8 +335,35 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
 !  ! ======
 
   
+  integer, allocatable :: indices(:), exc_degree(:), iorder(:)
   integer(bit_kind), allocatable:: preinteresting_det(:,:,:)
-  allocate (preinteresting_det(N_int,2,N_det))
+  integer :: l_a, nmax
+  allocate (preinteresting_det(N_int,2,N_det), indices(N_det),  &
+            exc_degree(N_det_alpha_unique))
+  do i=1,N_det_alpha_unique
+    call get_excitation_degree_spin(psi_det_alpha_unique(1,i), &
+      psi_det_generators(1,1,i_generator), exc_degree(i), N_int)
+  enddo
+
+  k=1
+  do j=1,N_det_beta_unique
+    call get_excitation_degree_spin(psi_det_beta_unique(1,j), &
+      psi_det_generators(1,2,i_generator), nt, N_int)
+    if (nt > 4) cycle
+    do l_a=psi_bilinear_matrix_columns_loc(j), psi_bilinear_matrix_columns_loc(j+1)-1
+      i = psi_bilinear_matrix_rows(l_a)
+      if (nt + exc_degree(i) <= 4) then
+        indices(k) = psi_det_sorted_order(psi_bilinear_matrix_order(l_a))
+        k=k+1
+      endif
+    enddo
+  enddo
+  nmax=k-1
+  allocate(iorder(nmax))
+  do i=1,nmax
+    iorder(i) = i
+  enddo
+  call isort(indices,iorder,nmax)
 
   preinteresting(0) = 0
   prefullinteresting(0) = 0
@@ -343,7 +373,8 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
     negMask(i,2) = not(psi_det_generators(i,2,i_generator))
   end do
   
-  do i=1,N_det
+  do k=1,nmax
+    i = indices(k)
     mobMask(1,1) = iand(negMask(1,1), psi_det_sorted(1,1,i))
     mobMask(1,2) = iand(negMask(1,2), psi_det_sorted(1,2,i))
     nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2)) 

src/Determinants/determinants.irp.f

@@ -324,9 +324,12 @@ END_PROVIDER
  BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted, (N_int,2,psi_det_size) ]
 &BEGIN_PROVIDER [ double precision, psi_coef_sorted, (psi_det_size,N_states) ]
 &BEGIN_PROVIDER [ double precision, psi_average_norm_contrib_sorted, (psi_det_size) ]
+&BEGIN_PROVIDER [ integer, psi_det_sorted_order, (psi_det_size) ]
  implicit none
  BEGIN_DOC
  ! Wave function sorted by determinants contribution to the norm (state-averaged)
+ !
+ ! psi_det_sorted_order(i) -> k : index in psi_det
  END_DOC
  integer :: i,j,k
  integer, allocatable ::  iorder(:)
@@ -345,7 +348,9 @@ END_PROVIDER
     psi_coef_sorted(i,k) = psi_coef(iorder(i),k)
   enddo
   psi_average_norm_contrib_sorted(i) = -psi_average_norm_contrib_sorted(i)
+  psi_det_sorted_order(i) = i
  enddo
+ call iset_order(psi_det_sorted_order,iorder,N_det)
 
  deallocate(iorder)