Fixed make_s2_eigenfunction

config/gfortran_debug.cfg

@@ -13,7 +13,7 @@
 FC           : gfortran -g -ffree-line-length-none -I . -static-libgcc 
 LAPACK_LIB   : -llapack -lblas
 IRPF90       : irpf90
-IRPF90_FLAGS : --ninja --assert --align=32 
+IRPF90_FLAGS : --ninja --align=32 
 
 # Global options
 ################

src/Determinants/H_apply.irp.f

@@ -222,8 +222,12 @@ subroutine remove_duplicates_in_psi_det(found_duplicates)
     do while (bit_tmp(j)==bit_tmp(i))
       if (duplicate(j)) then
         j += 1
+        if (j > N_det) then
+          exit
+        else
           cycle
         endif
+      endif
       duplicate(j) = .True.
       do k=1,N_int
         if ( (psi_det_sorted_bit(k,1,i) /= psi_det_sorted_bit(k,1,j) ) &

src/Determinants/occ_pattern.irp.f

@@ -52,8 +52,8 @@ subroutine occ_pattern_to_dets(o,d,sze,n_alpha,Nint)
   integer(bit_kind),intent(out)  :: d(Nint,2,sze)
   
   integer                        :: i, k, nt, na, nd, amax
-  integer                        :: list_todo(n_alpha)
+  integer                        :: list_todo(2*n_alpha)
-  integer                        :: list_a(n_alpha)
+  integer                        :: list_a(2*n_alpha)
 
   amax = n_alpha
   do k=1,Nint
@@ -69,13 +69,24 @@ subroutine occ_pattern_to_dets(o,d,sze,n_alpha,Nint)
 
   sze = nd
   
+  integer :: ne(2), l
+  l=0
   do i=1,nd
+    ne(1) = 0
+    ne(2) = 0
+    l=l+1
     ! Doubly occupied orbitals
     do k=1,Nint
-      d(k,1,i) = ior(d(k,1,i),o(k,2))
+      d(k,1,l) = ior(d(k,1,i),o(k,2))
-      d(k,2,i) = ior(d(k,2,i),o(k,2))
+      d(k,2,l) = ior(d(k,2,i),o(k,2))
+      ne(1) += popcnt(d(k,1,l))
+      ne(2) += popcnt(d(k,2,l))
     enddo
+    if ( (ne(1) /= elec_alpha_num).or.(ne(2) /= elec_beta_num) ) then
+      l = l-1
+    endif
   enddo
+  sze = l
 
 end
 
@@ -169,11 +180,11 @@ end
   endif
   j = i+1
   do while (bit_tmp(j)==bit_tmp(i))
+    if (duplicate(j)) then
+      j+=1
       if (j>N_det) then
         exit
       endif
-    if (duplicate(j)) then
-      j+=1
       cycle
     endif
     duplicate(j) = .True.
@@ -185,6 +196,9 @@ end
       endif
     enddo
     j+=1
+    if (j>N_det) then
+      exit
+    endif
   enddo
  enddo
 
@@ -256,6 +270,20 @@ subroutine make_s2_eigenfunction
           det_buffer(k,1,N_det_new) = d(k,1,j)
           det_buffer(k,2,N_det_new) = d(k,2,j)
         enddo
+!  integer :: ne(2)
+!  ne(:) = 0
+!  do k=1,N_int
+!    ne(1) += popcnt(d(k,1,j))
+!    ne(2) += popcnt(d(k,2,j))
+!  enddo
+!  if (ne(1) /= elec_alpha_num) then
+!    call debug_det(d(1,1,j),N_int)
+!    stop "ALPHA"
+!  endif
+!  if (ne(2) /= elec_beta_num) then
+!    call debug_det(d(1,1,j),N_int)
+!    stop "BETA"
+!  endif
         if (N_det_new == bufsze) then
           call fill_H_apply_buffer_no_selection(bufsze,det_buffer,N_int,0)
           N_det_new = 0
@@ -274,7 +302,8 @@ subroutine make_s2_eigenfunction
   call copy_H_apply_buffer_to_wf
   SOFT_TOUCH N_det psi_coef psi_det
   print *,  'Added determinants for S^2'
-!  call remove_duplicates_in_psi_det
+!  logical :: found
+!  call remove_duplicates_in_psi_det(found)
 
 end