Fixed selection

plugins/Full_CI_ZMQ/pt2_stoch.irp.f

@@ -24,7 +24,7 @@ subroutine run
   E_CI_before = pt2_E0_denominator(1) + nuclear_repulsion
   threshold_selectors = 1.d0
   threshold_generators = 1d0 
-  relative_error = 1.d-3
+  relative_error = -1.d-3
 !  relative_error = 1.d-8
   call ZMQ_pt2(E_CI_before, pt2, relative_error)
   print *,  'Final step'

plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f

@@ -27,7 +27,7 @@ subroutine ZMQ_pt2(E, pt2,relative_error)
   double precision, external :: omp_get_wtime 
   double precision :: time
 
-  allocate(pt2_detail(N_states, N_det_generators), comb(N_det_generators), computed(N_det_generators), tbc(0:size_tbc))
+  allocate(pt2_detail(N_states, N_det_generators+1), comb(N_det_generators), computed(N_det_generators), tbc(0:size_tbc))
   sumabove = 0d0
   sum2above = 0d0
   Nabove = 0d0
@@ -256,11 +256,15 @@ subroutine pt2_collector(E, b, tbc, comb, Ncomb, computed, pt2_detail, sumabove,
       end do
 
       E0 = sum(pt2_detail(1,:first_det_of_teeth(tooth)-1))
-      prop = ((1d0 - dfloat(comb_teeth - tooth + 1) * comb_step) - pt2_cweight(first_det_of_teeth(tooth)-1))
+      if (tooth <= comb_teeth) then
-      prop = prop * pt2_weight_inv(first_det_of_teeth(tooth))
+        prop = ((1d0 - dfloat(comb_teeth - tooth + 1) * comb_step) - pt2_cweight(first_det_of_teeth(tooth)-1))
-      E0 += pt2_detail(1,first_det_of_teeth(tooth)) * prop
+        prop = prop * pt2_weight_inv(first_det_of_teeth(tooth))
-      avg = E0 + (sumabove(tooth) / Nabove(tooth))
+        E0 += pt2_detail(1,first_det_of_teeth(tooth)) * prop
-      eqt = sqrt(1d0 / (Nabove(tooth)-1) * abs(sum2above(tooth) / Nabove(tooth) - (sumabove(tooth)/Nabove(tooth))**2))
+        avg = E0 + (sumabove(tooth) / Nabove(tooth))
+        eqt = sqrt(1d0 / (Nabove(tooth)-1) * abs(sum2above(tooth) / Nabove(tooth) - (sumabove(tooth)/Nabove(tooth))**2))
+      else
+        eqt = 0.d0
+      endif
       call wall_time(time)
       if (dabs(eqt/avg) < relative_error) then
         ! Termination
@@ -387,16 +391,16 @@ subroutine get_carlo_workbatch(computed, comb, Ncomb, tbc)
   integer                        :: i, j, last_full, dets(comb_teeth)
   integer                        :: icount, n
   integer                        :: k, l
-  l=1
+  l=first_det_of_comb
   call RANDOM_NUMBER(comb)
   do i=1,size(comb)
     comb(i) = comb(i) * comb_step
     !DIR$ FORCEINLINE
     call add_comb(comb(i), computed, tbc, size_tbc, comb_teeth)
     Ncomb = i
+    if (tbc(0) == N_det_generators) return
     do while (computed(l))
       l=l+1
-      if (l == N_det_generators+1) return
     enddo
     k=tbc(0)+1
     tbc(k) = l
@@ -511,15 +515,25 @@ end subroutine
   pt2_weight(1) = psi_coef_generators(1,1)**2
   pt2_cweight(1) = psi_coef_generators(1,1)**2
   
-  do i=2,N_det_generators
+  do i=1,N_det_generators
     pt2_weight(i) = psi_coef_generators(i,1)**2
-    pt2_cweight(i) = pt2_cweight(i-1) + psi_coef_generators(i,1)**2
+  enddo
+
+  ! Important to loop backwards for numerical precision
+  pt2_cweight(N_det_generators) = pt2_weight(N_det_generators)
+  do i=N_det_generators-1,1,-1
+    pt2_cweight(i) = pt2_weight(i) + pt2_cweight(i+1) 
   end do
   
   do i=1,N_det_generators
-    pt2_weight(i)  = pt2_weight(i) / pt2_cweight(N_det_generators)
+    pt2_weight(i)  = pt2_weight(i) / pt2_cweight(1)
-    pt2_cweight(i) = pt2_cweight(i) / pt2_cweight(N_det_generators)
+    pt2_cweight(i) = pt2_cweight(i) / pt2_cweight(1)
   enddo
+
+  do i=1,N_det_generators-1
+    pt2_cweight(i) = 1.d0 - pt2_cweight(i+1) 
+  end do
+  pt2_cweight(N_det_generators) = 1.d0
   
   norm_left = 1d0
   

plugins/Full_CI_ZMQ/selection.irp.f

@@ -320,36 +320,21 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
   call bitstring_to_list_ab(hole    , hole_list    , N_holes    , N_int)
   call bitstring_to_list_ab(particle, particle_list, N_particles, N_int)
 
-!  ! ======
-!  ! If the subset doesn't exist, return
-!  logical :: will_compute
-!  will_compute = subset == 0
-!
-!  if (.not.will_compute) then
-!    maskInd = N_holes(1)*N_holes(2) + N_holes(2)*((N_holes(2)-1)/2) + N_holes(1)*((N_holes(1)-1)/2)
-!    will_compute = (maskInd >= subset)
-!    if (.not.will_compute) then
-!      return
-!    endif
-!  endif
-!  ! ======
-
-  
   integer, allocatable :: indices(:), exc_degree(:), iorder(:)
   integer(bit_kind), allocatable:: preinteresting_det(:,:,:)
   integer :: l_a, nmax
   allocate (preinteresting_det(N_int,2,N_det), indices(N_det),  &
-            exc_degree(N_det_alpha_unique))
+            exc_degree(max(N_det_alpha_unique,N_det_beta_unique)))
+  k=1
   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
+    if (nt > 2) 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
@@ -358,6 +343,27 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
       endif
     enddo
   enddo
+  
+  do i=1,N_det_beta_unique
+    call get_excitation_degree_spin(psi_det_beta_unique(1,i), &
+      psi_det_generators(1,2,i_generator), exc_degree(i), N_int)
+  enddo
+
+  do j=1,N_det_alpha_unique
+    call get_excitation_degree_spin(psi_det_alpha_unique(1,j), &
+      psi_det_generators(1,1,i_generator), nt, N_int)
+    if (nt > 1) cycle
+    do l_a=psi_bilinear_matrix_transp_rows_loc(j), psi_bilinear_matrix_transp_rows_loc(j+1)-1
+      i = psi_bilinear_matrix_transp_columns(l_a)
+      if (exc_degree(i) < 3) cycle
+      if (nt + exc_degree(i) <= 4) then
+        indices(k) = psi_det_sorted_order(                   &
+                        psi_bilinear_matrix_order(           &
+                          psi_bilinear_matrix_transp_order(l_a)))
+        k=k+1
+      endif
+    enddo
+  enddo
   nmax=k-1
   allocate(iorder(nmax))
   do i=1,nmax
@@ -365,6 +371,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
   enddo
   call isort(indices,iorder,nmax)
 
+
   preinteresting(0) = 0
   prefullinteresting(0) = 0
   

plugins/Full_CI_ZMQ/selection_buffer.irp.f

@@ -53,7 +53,7 @@ subroutine merge_selection_buffers(b1, b2)
   BEGIN_DOC
 ! Merges the selection buffers b1 and b2 into b2
   END_DOC
-  type(selection_buffer), intent(in) :: b1
+  type(selection_buffer), intent(inout) :: b1
   type(selection_buffer), intent(inout) :: b2
   integer(bit_kind), pointer     :: detmp(:,:,:)
   double precision, pointer      :: val(:)
@@ -62,6 +62,9 @@ subroutine merge_selection_buffers(b1, b2)
   allocate( val(size(b1%val)), detmp(N_int, 2, size(b1%det,3)) )
   i1=1
   i2=1
+  do while (b1%val(b1%cur) > b2%mini)
+    b1%cur = b1%cur-1
+  enddo
   do i=1,nmwen
     if ( (i1 > b1%cur).and.(i2 > b2%cur) ) then
       exit 

plugins/Full_CI_ZMQ/zmq_selection.irp.f

@@ -98,7 +98,7 @@ subroutine selection_collector(b, N, pt2)
 
   zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
   zmq_socket_pull = new_zmq_pull_socket()
-  call create_selection_buffer(N, N*2, b2)
+  call create_selection_buffer(N, N*8, b2)
   allocate(task_id(N_det_generators))
   more = 1
   pt2(:) = 0d0
@@ -107,7 +107,11 @@ subroutine selection_collector(b, N, pt2)
     call pull_selection_results(zmq_socket_pull, pt2_mwen, b2%val(1), b2%det(1,1,1), b2%cur, task_id, ntask)
 
     pt2 += pt2_mwen
-    call merge_selection_buffers(b2,b)
+    do i=1, b2%cur
+      call add_to_selection_buffer(b, b2%det(1,1,i), b2%val(i))
+      if (b2%val(i) > b%mini) exit
+    end do
+
     do i=1, ntask
       if(task_id(i) == 0) then
           print *,  "Error in collector"

src/Determinants/determinants.irp.f

@@ -348,9 +348,11 @@ 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)
+ do i=1,N_det
+   psi_det_sorted_order(iorder(i)) = i
+ enddo
+
 
  deallocate(iorder)
 

src/Determinants/spindeterminants.irp.f

@@ -438,7 +438,7 @@ BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order_reverse , (N_det) ]
   use bitmasks
   implicit none
   BEGIN_DOC
-! Order which allors to go from psi_bilinear_matrix to psi_det
+! Order which allows to go from psi_bilinear_matrix to psi_det
   END_DOC
   integer                        :: k
   !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(k)