Looks like everything is working.

src/csf/configuration_CI_sigma_helpers.irp.f

@@ -231,7 +231,7 @@
 
            ! SOMO
            NalphaIcfg += 1
-           !print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg
+           !print *,i,j,"|",NalphaIcfg
            alphasIcfg(1,1,NalphaIcfg) = Jsomo
            alphasIcfg(1,2,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1)
         endif

src/csf/sigma_vector.irp.f

@@ -908,6 +908,10 @@ subroutine calculate_sigma_vector_cfg_nst(psi_out, psi_in, n_st, sze, istart, ie
   real*8, external :: mo_two_e_integral
   real*8, external :: get_two_e_integral
   real*8          :: diag_energies(n_CSF)
+  !PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals
+
+
+  print *," sze = ",sze
   call calculate_preconditioner_cfg(diag_energies)
 
   MS = 0
@@ -1043,6 +1047,7 @@ subroutine calculate_sigma_vector_cfg_nst(psi_out, psi_in, n_st, sze, istart, ie
      endi   = psi_config_data(i,2)
 
      ! Returns all unique (checking the past) singly excited cfgs connected to I
+     Nalphas_Icfg = 0
      call obtain_associated_alphaI(i, Icfg, alphas_Icfg, Nalphas_Icfg)
      ! TODO : remove doubly excited for return
      ! Here we do 2x the loop. One to count for the size of the matrix, then we compute.

src/davidson/diagonalization_hcsf_dressed.irp.f

@@ -303,6 +303,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
         ! -----------------------------------
 
         !call convertWFfromCSFtoDET(N_st_diag,U_csf(1,shift+1),U)
+        PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals
         if ((sze > 100000).and.distributed_davidson) then
 
             !call convertWFfromCSFtoDET(N_st_diag,U_csf(1,shift+1),U)
@@ -311,6 +312,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
             !call convertWFfromDETtoCSF(N_st_diag,U,U_csf(1,shift+1))
             !call convertWFfromDETtoCSF(N_st_diag,W,W_csf(1,shift+1))
             !call calculate_sigma_vector_cfg_nst(W_csf(1,shift+1),U_csf(1,shift+1),N_st_diag,sze_csf,1,sze_csf,0,1)
+            ! TODO : psi_det_size ? for psi_det
             do kk=1,N_st_diag
                 call calculate_sigma_vector_cfg_nst(W_csf(1,shift+kk),U_csf(1,shift+kk),1,sze_csf,1,sze_csf,0,1)
             enddo

src/determinants/slater_rules.irp.f

@@ -977,12 +977,12 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
   call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx)
   if (Nstate == 1) then
 
-    do ii=1,idx(0)
-      i = idx(ii)
-      !DIR$ FORCEINLINE
-      call i_H_j(keys(1,1,i),key,Nint,hij)
-      i_H_psi_array(1) = i_H_psi_array(1) + coef(i,1)*hij
-    enddo
+   do ii=1,idx(0)
+     i = idx(ii)
+     !DIR$ FORCEINLINE
+     call i_H_j(keys(1,1,i),key,Nint,hij)
+     i_H_psi_array(1) = i_H_psi_array(1) + coef(i,1)*hij
+   enddo
 
   else