From 4e35f9dbf622eb9d2ccfb04653f53c90874533ff Mon Sep 17 00:00:00 2001
From: eginer <giner.emmanuel@gmail.com>
Date: Wed, 15 Mar 2023 11:55:03 +0100
Subject: [PATCH 1/3] does not work

---
 src/cipsi/selection.irp.f                    |  45 +--
 src/cipsi_tc_bi_ortho/cipsi.irp.f            |   2 +-
 src/cipsi_tc_bi_ortho/pt2.irp.f              |   2 +-
 src/cipsi_tc_bi_ortho/selection.irp.f        | 365 +++++++++----------
 src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f |   2 +-
 src/fci_tc_bi/diagonalize_ci.irp.f           |   5 +-
 src/fci_tc_bi/generators.irp.f               |   7 +-
 src/fci_tc_bi/selectors.irp.f                |   9 +-
 src/tc_bi_ortho/psi_det_tc_sorted.irp.f      |  29 +-
 src/tc_bi_ortho/tc_h_eigvectors.irp.f        |   3 +
 10 files changed, 203 insertions(+), 266 deletions(-)

diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f
index 62d7c52c..6f40a809 100644
--- a/src/cipsi/selection.irp.f
+++ b/src/cipsi/selection.irp.f
@@ -312,9 +312,6 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
   end do
   deallocate(indices)
 
-!  !$OMP CRITICAL
-!  print *,  'Step1: ', i_generator, preinteresting(0)
-!  !$OMP END CRITICAL
 
   allocate(banned(mo_num, mo_num,2), bannedOrb(mo_num, 2))
   allocate (mat(N_states, mo_num, mo_num))
@@ -466,17 +463,8 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
           fullinteresting(sze+1) = i
         end if
       end do
-
       allocate (fullminilist (N_int, 2, fullinteresting(0)), &
                     minilist (N_int, 2,     interesting(0)) )
-!      if(pert_2rdm)then
-!        allocate(coef_fullminilist_rev(N_states,fullinteresting(0)))
-!        do i=1,fullinteresting(0)
-!          do j = 1, N_states
-!            coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j)
-!          enddo
-!        enddo
-!      endif
 
       do i=1,fullinteresting(0)
         fullminilist(:,:,i) = psi_det_sorted(:,:,fullinteresting(i))
@@ -524,33 +512,19 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
 
             call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting)
             if(fullMatch) cycle
-! !$OMP CRITICAL
-!  print *,  'Step3: ', i_generator, h1, interesting(0)
-! !$OMP END CRITICAL
 
             call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting)
-
-!            if(.not.pert_2rdm)then
-             call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf)
-!            else
-!             call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0))
-!            endif
+            call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf)
           end if
         enddo
         if(s1 /= s2) monoBdo = .false.
       enddo
       deallocate(fullminilist,minilist)
-!      if(pert_2rdm)then
-!        deallocate(coef_fullminilist_rev)
-!      endif
     enddo
   enddo
   deallocate(preinteresting, prefullinteresting, interesting, fullinteresting)
   deallocate(banned, bannedOrb,mat)
 end subroutine
-
-
-
 subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf)
   use bitmasks
   use selection_types
@@ -606,18 +580,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
 ! to a determinant of the future. In that case, the determinant will be
 ! detected as already generated when generating in the future with a
 ! double excitation.
-!
-!      if (.not.do_singles) then
-!        if ((h1 == p1) .or. (h2 == p2)) then
-!          cycle
-!        endif
-!      endif
-!
-!      if (.not.do_doubles) then
-!        if ((h1 /= p1).and.(h2 /= p2)) then
-!          cycle
-!        endif
-!      endif
 ! -----
 
       if(bannedOrb(p2, s2)) cycle
@@ -974,13 +936,10 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
 
         call get_mask_phase(psi_det_sorted(1,1,interesting(i)), phasemask,N_int)
         if(nt == 4) then
-!          call get_d2_reference(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
           call get_d2(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
         else if(nt == 3) then
-!          call get_d1_reference(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
           call get_d1(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
         else
-!          call get_d0_reference(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
           call get_d0(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
         end if
     else if(nt == 4) then
@@ -1540,8 +1499,6 @@ subroutine past_d2(banned, p, sp)
   end if
 end
 
-
-
 subroutine spot_isinwf(mask, det, i_gen, N, banned, fullMatch, interesting)
   use bitmasks
   implicit none
diff --git a/src/cipsi_tc_bi_ortho/cipsi.irp.f b/src/cipsi_tc_bi_ortho/cipsi.irp.f
index b1941068..fb907cb3 100644
--- a/src/cipsi_tc_bi_ortho/cipsi.irp.f
+++ b/src/cipsi_tc_bi_ortho/cipsi.irp.f
@@ -64,7 +64,7 @@ subroutine run_cipsi
   endif
 
   if (N_det > N_det_max) then
-    psi_det(1:N_int,1:2,1:N_det) = psi_det_sorted_tc_gen(1:N_int,1:2,1:N_det)
+    psi_det(1:N_int,1:2,1:N_det) = psi_det_generators(1:N_int,1:2,1:N_det)
     psi_coef(1:N_det,1:N_states) = psi_coef_sorted_tc_gen(1:N_det,1:N_states)
     N_det = N_det_max
     soft_touch N_det psi_det psi_coef
diff --git a/src/cipsi_tc_bi_ortho/pt2.irp.f b/src/cipsi_tc_bi_ortho/pt2.irp.f
index e7dca456..13b4dff4 100644
--- a/src/cipsi_tc_bi_ortho/pt2.irp.f
+++ b/src/cipsi_tc_bi_ortho/pt2.irp.f
@@ -52,7 +52,7 @@ subroutine pt2_tc_bi_ortho
 !  call routine_save_right
 
   if (N_det > N_det_max) then
-    psi_det(1:N_int,1:2,1:N_det) = psi_det_sorted_tc_gen(1:N_int,1:2,1:N_det)
+    psi_det(1:N_int,1:2,1:N_det) = psi_det_generators(1:N_int,1:2,1:N_det)
     psi_coef(1:N_det,1:N_states) = psi_coef_sorted_tc_gen(1:N_det,1:N_states)
     N_det = N_det_max
     soft_touch N_det psi_det psi_coef
diff --git a/src/cipsi_tc_bi_ortho/selection.irp.f b/src/cipsi_tc_bi_ortho/selection.irp.f
index 633ca815..393023f2 100644
--- a/src/cipsi_tc_bi_ortho/selection.irp.f
+++ b/src/cipsi_tc_bi_ortho/selection.irp.f
@@ -91,7 +91,6 @@ subroutine select_connected(i_generator,E0,pt2_data,b,subset,csubset)
 end subroutine select_connected
 
 
-
 double precision function get_phase_bi(phasemask, s1, s2, h1, p1, h2, p2, Nint)
   use bitmasks
   implicit none
@@ -136,7 +135,7 @@ double precision function get_phase_bi(phasemask, s1, s2, h1, p1, h2, p2, Nint)
 end
 
 
-subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock_diag_tmp, E0, pt2_data, buf, subset, csubset)
+subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2_data,buf,subset,csubset)
   use bitmasks
   use selection_types
   implicit none
@@ -266,7 +265,6 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock
   enddo
 
   do k = 1, nmax
-
     i = indices(k)
     mobMask(1,1) = iand(negMask(1,1), psi_det_sorted_tc(1,1,i))
     mobMask(1,2) = iand(negMask(1,2), psi_det_sorted_tc(1,2,i))
@@ -304,10 +302,10 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock
         prefullinteresting(sze+1) = i
       endif
     endif
-
   enddo
   deallocate(indices)
 
+
   allocate( banned(mo_num, mo_num,2), bannedOrb(mo_num, 2) )
   allocate( mat(N_states, mo_num, mo_num) )
   allocate( mat_l(N_states, mo_num, mo_num), mat_r(N_states, mo_num, mo_num) )
@@ -463,17 +461,11 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock
                     minilist (N_int, 2,     interesting(0))  )
 
       do i = 1, fullinteresting(0)
-        do k = 1, N_int
-          fullminilist(k,1,i) = psi_det_sorted_tc(k,1,fullinteresting(i))
-          fullminilist(k,2,i) = psi_det_sorted_tc(k,2,fullinteresting(i))
-        enddo
+        fullminilist(:,:,i) = psi_det_sorted_tc(:,:,fullinteresting(i))
       enddo
 
       do i = 1, interesting(0)
-        do k = 1, N_int
-          minilist(k,1,i) = psi_det_sorted_tc(k,1,interesting(i))
-          minilist(k,2,i) = psi_det_sorted_tc(k,2,interesting(i))
-        enddo
+        minilist(:,:,i) = psi_det_sorted_tc(:,:,interesting(i))
       enddo
 
       do s2 = s1, 2
@@ -516,196 +508,19 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock
             if(fullMatch) cycle
 
             call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting, mat_l, mat_r)
-
             call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf, mat_l, mat_r)
           endif
-
         enddo
-
         if(s1 /= s2) monoBdo = .false.
       enddo
-
       deallocate(fullminilist, minilist)
-
     enddo
   enddo
-
   deallocate(preinteresting, prefullinteresting, interesting, fullinteresting)
   deallocate(banned, bannedOrb,mat)
   deallocate(mat_l, mat_r)
-
 end subroutine select_singles_and_doubles
-
-! ---
-
-subroutine spot_isinwf(mask, det, i_gen, N, banned, fullMatch, interesting)
-
-  use bitmasks
-  implicit none
-
-  BEGIN_DOC
-  ! Identify the determinants in det which are in the internal space. These are
-  ! the determinants that can be produced by creating two particles on the mask.
-  END_DOC
-
-  integer, intent(in) :: i_gen, N
-  integer, intent(in) :: interesting(0:N)
-  integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N)
-  logical, intent(inout) :: banned(mo_num, mo_num)
-  logical, intent(out) :: fullMatch
-
-  integer :: i, j, na, nb, list(3)
-  integer(bit_kind) :: myMask(N_int, 2), negMask(N_int, 2)
-
-  fullMatch = .false.
-
-  do i=1,N_int
-    negMask(i,1) = not(mask(i,1))
-    negMask(i,2) = not(mask(i,2))
-  end do
-
-  genl : do i=1, N
-    ! If det(i) can't be generated by the mask, cycle
-    do j=1, N_int
-      if(iand(det(j,1,i), mask(j,1)) /= mask(j, 1)) cycle genl
-      if(iand(det(j,2,i), mask(j,2)) /= mask(j, 2)) cycle genl
-    end do
-
-    ! If det(i) < det(i_gen), it hs already been considered
-    if(interesting(i) < i_gen) then
-      fullMatch = .true.
-      return
-    end if
-
-    ! Identify the particles
-    do j=1, N_int
-      myMask(j, 1) = iand(det(j, 1, i), negMask(j, 1))
-      myMask(j, 2) = iand(det(j, 2, i), negMask(j, 2))
-    end do
-
-    call bitstring_to_list_in_selection(myMask(1,1), list(1), na, N_int)
-    call bitstring_to_list_in_selection(myMask(1,2), list(na+1), nb, N_int)
-    banned(list(1), list(2)) = .true.
-  end do genl
-
-end subroutine spot_isinwf
-
-! ---
-
-subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting, mat_l, mat_r)
-
-  BEGIN_DOC
-  ! Computes the contributions A(r,s) by
-  ! comparing the external determinant to all the internal determinants det(i).
-  ! an applying two particles (r,s) to the mask.
-  END_DOC
-
-  use bitmasks
-  implicit none
-
-  integer, intent(in)             :: sp, i_gen, N_sel
-  integer, intent(in)             :: interesting(0:N_sel)
-  integer(bit_kind),intent(in)    :: mask(N_int, 2), det(N_int, 2, N_sel)
-  logical, intent(inout)          :: bannedOrb(mo_num, 2), banned(mo_num, mo_num, 2)
-  double precision, intent(inout) :: mat(N_states, mo_num, mo_num)
-  double precision, intent(inout) :: mat_l(N_states, mo_num, mo_num), mat_r(N_states, mo_num, mo_num)
-
-  integer                         :: i, ii, j, k, l, h(0:2,2), p(0:4,2), nt
-  integer(bit_kind)               :: perMask(N_int, 2), mobMask(N_int, 2), negMask(N_int, 2)
-  integer(bit_kind)               :: phasemask(N_int,2)
-
-
-  PROVIDE psi_selectors_coef_transp_tc psi_det_sorted_tc
-  PROVIDE psi_selectors_rcoef_bi_orth_transp psi_selectors_lcoef_bi_orth_transp
-
-
-  mat   = 0d0
-  mat_l = 0d0
-  mat_r = 0d0
-
-  do i = 1, N_int
-    negMask(i,1) = not(mask(i,1))
-    negMask(i,2) = not(mask(i,2))
-  end do
-
-  do i = 1, N_sel
-    if(interesting(i) < 0) then
-      stop 'prefetch interesting(i) and det(i)'
-    endif
-
-    mobMask(1,1) = iand(negMask(1,1), det(1,1,i))
-    mobMask(1,2) = iand(negMask(1,2), det(1,2,i))
-    nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2))
-
-    if(nt > 4) cycle
-
-    do j = 2, N_int
-      mobMask(j,1) = iand(negMask(j,1), det(j,1,i))
-      mobMask(j,2) = iand(negMask(j,2), det(j,2,i))
-      nt = nt + popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2))
-    enddo
-
-    if(nt > 4) cycle
-
-    if (interesting(i) == i_gen) then
-      if(sp == 3) then
-        do k = 1, mo_num
-          do j = 1, mo_num
-            banned(j,k,2) = banned(k,j,1)
-          enddo
-        enddo
-      else
-        do k = 1, mo_num
-          do l = k+1, mo_num
-            banned(l,k,1) = banned(k,l,1)
-          enddo
-        enddo
-      endif
-    endif
-
-    if (interesting(i) >= i_gen) then
-
-        call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int)
-        call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int)
-
-        perMask(1,1) = iand(mask(1,1), not(det(1,1,i)))
-        perMask(1,2) = iand(mask(1,2), not(det(1,2,i)))
-        do j=2,N_int
-          perMask(j,1) = iand(mask(j,1), not(det(j,1,i)))
-          perMask(j,2) = iand(mask(j,2), not(det(j,2,i)))
-        end do
-!        call get_d3_h  ( det(1,1,i), bannedOrb, banned, mat         , mask, p, sp, psi_selectors_coef_transp_tc (1, interesting(i)) )
-!        call get_d3_htc( det(1,1,i), bannedOrb, banned, mat_r, mat_l, mask, p, sp, psi_selectors_rcoef_bi_orth_transp(1, interesting(i)) &
-!                       , psi_selectors_lcoef_bi_orth_transp(1, interesting(i)) )
-
-        call bitstring_to_list_in_selection(perMask(1,1), h(1,1), h(0,1), N_int)
-        call bitstring_to_list_in_selection(perMask(1,2), h(1,2), h(0,2), N_int)
-
-        call get_mask_phase(psi_det_sorted_tc(1,1,interesting(i)), phasemask,N_int)
-        if(nt == 4) then
-          call get_d2_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
-        elseif(nt == 3) then
-          call get_d1_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
-        else
-          call get_d0_new (det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
-        endif
-    elseif(nt == 4) then
-        call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int)
-        call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int)
-        call past_d2(banned, p, sp)
-    elseif(nt == 3) then
-        call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int)
-        call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int)
-        call past_d1(bannedOrb, p)
-    endif
-  enddo
-
-end subroutine splash_pq
-
-! ---
-
 subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf, mat_l, mat_r)
-
   use bitmasks
   use selection_types
   implicit none
@@ -740,7 +555,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
 
 
   PROVIDE dominant_dets_of_cfgs N_dominant_dets_of_cfgs
-
   do jstate = 1, N_states
     do istate = 1, N_states
       s_weight(istate,jstate) = dsqrt(selection_weight(istate)*selection_weight(jstate))
@@ -780,11 +594,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
       !  endif
       !endif
 
-      ! MANU: ERREUR dans les calculs puisque < I | H | J > = 0
-      ! n'implique pas < I | H_TC | J > = 0 ??
-      !val = maxval(abs(mat(1:N_states, p1, p2)))
-      !if( val == 0d0) cycle
-
       call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
 
       if(do_only_cas) then
@@ -811,7 +620,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
       if(excitation_max >= 0) then
         do_cycle = .True.
         if(excitation_ref == 1) then
-          call get_excitation_degree(HF_bitmask, det(1,1), degree, N_int)
+          call get_excitation_degree(HF_bitmask,det(1,1),degree,N_int)
           do_cycle = do_cycle .and. (degree > excitation_max)
         elseif(excitation_ref == 2) then
           do k = 1, N_dominant_dets_of_cfgs
@@ -995,12 +804,118 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
 
     enddo ! end do p2
   enddo   ! end do p1
-
 end subroutine fill_buffer_double
-
 ! ---
+subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting, mat_l, mat_r)
+
+  BEGIN_DOC
+  ! Computes the contributions A(r,s) by
+  ! comparing the external determinant to all the internal determinants det(i).
+  ! an applying two particles (r,s) to the mask.
+  END_DOC
+
+  use bitmasks
+  implicit none
+
+  integer, intent(in)             :: sp, i_gen, N_sel
+  integer, intent(in)             :: interesting(0:N_sel)
+  integer(bit_kind),intent(in)    :: mask(N_int, 2), det(N_int, 2, N_sel)
+  logical, intent(inout)          :: bannedOrb(mo_num, 2), banned(mo_num, mo_num, 2)
+  double precision, intent(inout) :: mat(N_states, mo_num, mo_num)
+  double precision, intent(inout) :: mat_l(N_states, mo_num, mo_num), mat_r(N_states, mo_num, mo_num)
+
+  integer                         :: i, ii, j, k, l, h(0:2,2), p(0:4,2), nt
+  integer(bit_kind)               :: perMask(N_int, 2), mobMask(N_int, 2), negMask(N_int, 2)
+  integer(bit_kind)               :: phasemask(N_int,2)
 
 
+  PROVIDE psi_selectors_coef_transp_tc psi_det_sorted_tc
+  PROVIDE psi_selectors_rcoef_bi_orth_transp psi_selectors_lcoef_bi_orth_transp
+
+
+  mat   = 0d0
+  mat_l = 0d0
+  mat_r = 0d0
+
+  do i = 1, N_int
+    negMask(i,1) = not(mask(i,1))
+    negMask(i,2) = not(mask(i,2))
+  end do
+
+  do i = 1, N_sel
+    if(interesting(i) < 0) then
+      stop 'prefetch interesting(i) and det(i)'
+    endif
+
+    mobMask(1,1) = iand(negMask(1,1), det(1,1,i))
+    mobMask(1,2) = iand(negMask(1,2), det(1,2,i))
+    nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2))
+
+    if(nt > 4) cycle
+
+    do j = 2, N_int
+      mobMask(j,1) = iand(negMask(j,1), det(j,1,i))
+      mobMask(j,2) = iand(negMask(j,2), det(j,2,i))
+      nt = nt + popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2))
+    enddo
+
+    if(nt > 4) cycle
+
+    if (interesting(i) == i_gen) then
+      if(sp == 3) then
+        do k = 1, mo_num
+          do j = 1, mo_num
+            banned(j,k,2) = banned(k,j,1)
+          enddo
+        enddo
+      else
+        do k = 1, mo_num
+          do l = k+1, mo_num
+            banned(l,k,1) = banned(k,l,1)
+          enddo
+        enddo
+      endif
+    endif
+
+    if (interesting(i) >= i_gen) then
+
+        call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int)
+        call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int)
+
+        perMask(1,1) = iand(mask(1,1), not(det(1,1,i)))
+        perMask(1,2) = iand(mask(1,2), not(det(1,2,i)))
+        do j=2,N_int
+          perMask(j,1) = iand(mask(j,1), not(det(j,1,i)))
+          perMask(j,2) = iand(mask(j,2), not(det(j,2,i)))
+        end do
+!        call get_d3_h  ( det(1,1,i), bannedOrb, banned, mat         , mask, p, sp, psi_selectors_coef_transp_tc (1, interesting(i)) )
+!        call get_d3_htc( det(1,1,i), bannedOrb, banned, mat_r, mat_l, mask, p, sp, psi_selectors_rcoef_bi_orth_transp(1, interesting(i)) &
+!                       , psi_selectors_lcoef_bi_orth_transp(1, interesting(i)) )
+
+        call bitstring_to_list_in_selection(perMask(1,1), h(1,1), h(0,1), N_int)
+        call bitstring_to_list_in_selection(perMask(1,2), h(1,2), h(0,2), N_int)
+
+        call get_mask_phase(psi_det_sorted_tc(1,1,interesting(i)), phasemask,N_int)
+        if(nt == 4) then
+          call get_d2_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
+        elseif(nt == 3) then
+          call get_d1_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
+        else
+          call get_d0_new (det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
+        endif
+    elseif(nt == 4) then
+        call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int)
+        call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int)
+        call past_d2(banned, p, sp)
+    elseif(nt == 3) then
+        call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int)
+        call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int)
+        call past_d1(bannedOrb, p)
+    endif
+  enddo
+
+end subroutine splash_pq
+! ---
 subroutine past_d1(bannedOrb, p)
 
   use bitmasks
@@ -1043,9 +958,61 @@ subroutine past_d2(banned, p, sp)
       end do
     end do
   end if
-
 end subroutine past_d2
 
+subroutine spot_isinwf(mask, det, i_gen, N, banned, fullMatch, interesting)
+
+  use bitmasks
+  implicit none
+
+  BEGIN_DOC
+  ! Identify the determinants in det which are in the internal space. These are
+  ! the determinants that can be produced by creating two particles on the mask.
+  END_DOC
+
+  integer, intent(in) :: i_gen, N
+  integer, intent(in) :: interesting(0:N)
+  integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N)
+  logical, intent(inout) :: banned(mo_num, mo_num)
+  logical, intent(out) :: fullMatch
+
+  integer :: i, j, na, nb, list(3)
+  integer(bit_kind) :: myMask(N_int, 2), negMask(N_int, 2)
+
+  fullMatch = .false.
+
+  do i=1,N_int
+    negMask(i,1) = not(mask(i,1))
+    negMask(i,2) = not(mask(i,2))
+  end do
+
+  genl : do i=1, N
+    ! If det(i) can't be generated by the mask, cycle
+    do j=1, N_int
+      if(iand(det(j,1,i), mask(j,1)) /= mask(j, 1)) cycle genl
+      if(iand(det(j,2,i), mask(j,2)) /= mask(j, 2)) cycle genl
+    end do
+
+    ! If det(i) < det(i_gen), it hs already been considered
+    if(interesting(i) < i_gen) then
+      fullMatch = .true.
+      return
+    end if
+
+    ! Identify the particles
+    do j=1, N_int
+      myMask(j, 1) = iand(det(j, 1, i), negMask(j, 1))
+      myMask(j, 2) = iand(det(j, 2, i), negMask(j, 2))
+    end do
+
+!    call debug_det(myMask, N_int)
+    call bitstring_to_list_in_selection(myMask(1,1), list(1), na, N_int)
+    call bitstring_to_list_in_selection(myMask(1,2), list(na+1), nb, N_int)
+    banned(list(1), list(2)) = .true.
+  end do genl
+
+end subroutine spot_isinwf
+
 ! ---
 
 subroutine bitstring_to_list_in_selection( string, list, n_elements, Nint)
diff --git a/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f b/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f
index e67287d3..64e7e6ba 100644
--- a/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f
+++ b/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f
@@ -54,7 +54,7 @@ subroutine run_stochastic_cipsi
 
 
 !  if (N_det > N_det_max) then
-!    psi_det(1:N_int,1:2,1:N_det) = psi_det_sorted_tc_gen(1:N_int,1:2,1:N_det)
+!    psi_det(1:N_int,1:2,1:N_det) = psi_det_generators(1:N_int,1:2,1:N_det)
 !    psi_coef(1:N_det,1:N_states) = psi_coef_sorted_tc_gen(1:N_det,1:N_states)
 !    N_det = N_det_max
 !    soft_touch N_det psi_det psi_coef
diff --git a/src/fci_tc_bi/diagonalize_ci.irp.f b/src/fci_tc_bi/diagonalize_ci.irp.f
index 56c561ac..c8369e93 100644
--- a/src/fci_tc_bi/diagonalize_ci.irp.f
+++ b/src/fci_tc_bi/diagonalize_ci.irp.f
@@ -49,9 +49,8 @@ subroutine diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2)
       psi_coef(i,j) = dabs(psi_l_coef_bi_ortho(i,j) * psi_r_coef_bi_ortho(i,j))   
     enddo
   enddo
-  SOFT_TOUCH  eigval_left_tc_bi_orth  eigval_right_tc_bi_orth  leigvec_tc_bi_orth  reigvec_tc_bi_orth norm_ground_left_right_bi_orth psi_coef  psi_l_coef_bi_ortho psi_r_coef_bi_ortho
-
-
+  SOFT_TOUCH eigval_left_tc_bi_orth  eigval_right_tc_bi_orth leigvec_tc_bi_orth reigvec_tc_bi_orth norm_ground_left_right_bi_orth 
+  SOFT_TOUCH psi_l_coef_bi_ortho psi_r_coef_bi_ortho psi_coef  
 
   call save_tc_bi_ortho_wavefunction
 end
diff --git a/src/fci_tc_bi/generators.irp.f b/src/fci_tc_bi/generators.irp.f
index 55c0cbb9..bf972423 100644
--- a/src/fci_tc_bi/generators.irp.f
+++ b/src/fci_tc_bi/generators.irp.f
@@ -43,9 +43,14 @@ END_PROVIDER
  ! For Single reference wave functions, the generator is the
  ! Hartree-Fock determinant
  END_DOC
- psi_det_sorted_tc_gen = psi_det_sorted_tc
+ psi_det_sorted_tc_gen =  psi_det_sorted_tc
  psi_coef_sorted_tc_gen = psi_coef_sorted_tc
  psi_det_sorted_tc_gen_order = psi_det_sorted_tc_order
+ integer :: i
+! do i = 1,N_det
+!  print*,'i = ',i
+!  call debug_det(psi_det_sorted_tc(1,1,i),N_int)
+! enddo
 END_PROVIDER
 
 
diff --git a/src/fci_tc_bi/selectors.irp.f b/src/fci_tc_bi/selectors.irp.f
index 3830927b..94aa4b01 100644
--- a/src/fci_tc_bi/selectors.irp.f
+++ b/src/fci_tc_bi/selectors.irp.f
@@ -47,13 +47,20 @@ END_PROVIDER
   enddo
   do k=1,N_states
     do i=1,N_det_selectors
-      psi_selectors_coef(i,k) = psi_coef_sorted_tc_gen(i,k)
+      psi_selectors_coef(i,k) =  dsqrt(dabs(psi_l_coef_sorted_bi_ortho(i,k) * psi_r_coef_sorted_bi_ortho(i,k)))
       psi_selectors_coef_tc(i,1,k) = psi_l_coef_sorted_bi_ortho(i,k)
       psi_selectors_coef_tc(i,2,k) = psi_r_coef_sorted_bi_ortho(i,k)
+!      psi_selectors_coef_tc(i,1,k) = psi_l_coef_bi_ortho(i,k)
+!      psi_selectors_coef_tc(i,2,k) = psi_r_coef_bi_ortho(i,k)
 !      psi_selectors_coef_tc(i,1,k) = 1.d0
 !      psi_selectors_coef_tc(i,2,k) = 1.d0
     enddo
   enddo
+  print*,'selectors '
+  do i = 1, N_det_selectors
+   print*,i,dabs(psi_selectors_coef_tc(i,1,1)*psi_selectors_coef_tc(i,2,1))
+   call debug_det(psi_selectors(1,1,i),N_int)
+  enddo
 
 END_PROVIDER
 
diff --git a/src/tc_bi_ortho/psi_det_tc_sorted.irp.f b/src/tc_bi_ortho/psi_det_tc_sorted.irp.f
index 42617557..35c78468 100644
--- a/src/tc_bi_ortho/psi_det_tc_sorted.irp.f
+++ b/src/tc_bi_ortho/psi_det_tc_sorted.irp.f
@@ -10,7 +10,10 @@ BEGIN_PROVIDER [ double precision, psi_average_norm_contrib_tc, (psi_det_size) ]
     
   psi_average_norm_contrib_tc(:) = 0.d0
   do k=1,N_states
+    print*,'in psi_average_norm_contrib_tc'
     do i=1,N_det
+      print*,i,dabs(psi_l_coef_bi_ortho(i,k)*psi_r_coef_bi_ortho(i,k))
+      call debug_det(psi_det(1,1,i),N_int)
       psi_average_norm_contrib_tc(i) = psi_average_norm_contrib_tc(i) +    &
           dabs(psi_l_coef_bi_ortho(i,k)*psi_r_coef_bi_ortho(i,k))*state_average_weight(k)
     enddo
@@ -26,11 +29,18 @@ END_PROVIDER
 &BEGIN_PROVIDER [ double precision, psi_coef_sorted_tc, (psi_det_size,N_states) ]
 &BEGIN_PROVIDER [ double precision, psi_average_norm_contrib_sorted_tc, (psi_det_size) ]
 &BEGIN_PROVIDER [ integer, psi_det_sorted_tc_order, (psi_det_size) ]
+&BEGIN_PROVIDER [double precision, psi_r_coef_sorted_bi_ortho, (psi_det_size, N_states)]
+&BEGIN_PROVIDER [double precision, psi_l_coef_sorted_bi_ortho, (psi_det_size, N_states)]
    implicit none
    BEGIN_DOC
    ! Wave function sorted by determinants contribution to the norm (state-averaged)
    !
    ! psi_det_sorted_tc_order(i) -> k : index in psi_det
+   !
+   ! psi_r_coef_sorted_bi_ortho : right coefficients corresponding to psi_det_sorted_tc
+   !
+   ! psi_l_coef_sorted_bi_ortho : left  coefficients corresponding to psi_det_sorted_tc
+ END_DOC
    END_DOC
    integer                        :: i,j,k
    integer, allocatable           :: iorder(:)
@@ -39,7 +49,7 @@ END_PROVIDER
      psi_average_norm_contrib_sorted_tc(i) = -psi_average_norm_contrib_tc(i)
      iorder(i) = i
    enddo
-!   call dsort(psi_average_norm_contrib_sorted_tc,iorder,N_det)
+   call dsort(psi_average_norm_contrib_sorted_tc,iorder,N_det)
    do i=1,N_det
      do j=1,N_int
        psi_det_sorted_tc(j,1,i) = psi_det(j,1,iorder(i))
@@ -66,24 +76,13 @@ END_PROVIDER
    psi_average_norm_contrib_sorted_tc(N_det+1:psi_det_size) = 0.d0
    psi_det_sorted_tc_order(N_det+1:psi_det_size) = 0
 
-   deallocate(iorder)
-
-END_PROVIDER
-
- BEGIN_PROVIDER [double precision, psi_r_coef_sorted_bi_ortho, (psi_det_size, N_states)]
-&BEGIN_PROVIDER [double precision, psi_l_coef_sorted_bi_ortho, (psi_det_size, N_states)]
- BEGIN_DOC
- ! psi_r_coef_sorted_bi_ortho : right coefficients corresponding to psi_det_sorted_tc
- ! psi_l_coef_sorted_bi_ortho : left  coefficients corresponding to psi_det_sorted_tc
- END_DOC
-   implicit none
-   integer                       :: i, j, k
    psi_r_coef_sorted_bi_ortho = 0.d0
    psi_l_coef_sorted_bi_ortho = 0.d0
    do i = 1, N_det
-    psi_r_coef_sorted_bi_ortho(i,1) = psi_r_coef_bi_ortho(psi_det_sorted_tc_order(i),1)
-    psi_l_coef_sorted_bi_ortho(i,1) = psi_l_coef_bi_ortho(psi_det_sorted_tc_order(i),1)
+    psi_r_coef_sorted_bi_ortho(i,1:N_states) = psi_r_coef_bi_ortho(psi_det_sorted_tc_order(i),1:N_states)
+    psi_l_coef_sorted_bi_ortho(i,1:N_states) = psi_l_coef_bi_ortho(psi_det_sorted_tc_order(i),1:N_states)
    enddo
+   deallocate(iorder)
 
 END_PROVIDER
 
diff --git a/src/tc_bi_ortho/tc_h_eigvectors.irp.f b/src/tc_bi_ortho/tc_h_eigvectors.irp.f
index d39b7a29..d12bbb4e 100644
--- a/src/tc_bi_ortho/tc_h_eigvectors.irp.f
+++ b/src/tc_bi_ortho/tc_h_eigvectors.irp.f
@@ -133,7 +133,10 @@
   call bi_normalize(leigvec_tc_bi_orth,reigvec_tc_bi_orth,size(reigvec_tc_bi_orth,1),N_det,N_states)
    print*,'leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1) = ',leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1)
    norm_ground_left_right_bi_orth = 0.d0
+   print*,'In diago'
    do j = 1, N_det
+    print*,j,dabs(leigvec_tc_bi_orth(j,1) * reigvec_tc_bi_orth(j,1))
+    call debug_det(psi_det(1,1,j),N_int)
     norm_ground_left_right_bi_orth += leigvec_tc_bi_orth(j,1) * reigvec_tc_bi_orth(j,1)
    enddo
    print*,'norm l/r = ',norm_ground_left_right_bi_orth

From d1068047e8df56adbdd3006a114b88e9995fb28f Mon Sep 17 00:00:00 2001
From: eginer <giner.emmanuel@gmail.com>
Date: Wed, 15 Mar 2023 14:33:10 +0100
Subject: [PATCH 2/3] trying to debug some psi_det_generators stuffs in
 fci_tc_bi

---
 src/cipsi_tc_bi_ortho/selection.irp.f   | 11 ++++---
 src/fci_tc_bi/generators.irp.f          |  8 +++++
 src/fci_tc_bi/selectors.irp.f           | 15 ++++++++++
 src/tc_bi_ortho/psi_det_tc_sorted.irp.f | 40 ++++++++++++++++++++++++-
 src/tc_bi_ortho/tc_h_eigvectors.irp.f   |  3 ++
 5 files changed, 72 insertions(+), 5 deletions(-)

diff --git a/src/cipsi_tc_bi_ortho/selection.irp.f b/src/cipsi_tc_bi_ortho/selection.irp.f
index 633ca815..0c3f0451 100644
--- a/src/cipsi_tc_bi_ortho/selection.irp.f
+++ b/src/cipsi_tc_bi_ortho/selection.irp.f
@@ -464,15 +464,15 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock
 
       do i = 1, fullinteresting(0)
         do k = 1, N_int
-          fullminilist(k,1,i) = psi_det_sorted_tc(k,1,fullinteresting(i))
-          fullminilist(k,2,i) = psi_det_sorted_tc(k,2,fullinteresting(i))
+          fullminilist(k,1,i) = psi_selectors(k,1,fullinteresting(i))
+          fullminilist(k,2,i) = psi_selectors(k,2,fullinteresting(i))
         enddo
       enddo
 
       do i = 1, interesting(0)
         do k = 1, N_int
-          minilist(k,1,i) = psi_det_sorted_tc(k,1,interesting(i))
-          minilist(k,2,i) = psi_det_sorted_tc(k,2,interesting(i))
+          minilist(k,1,i) = psi_selectors(k,1,interesting(i))
+          minilist(k,2,i) = psi_selectors(k,2,interesting(i))
         enddo
       enddo
 
@@ -628,7 +628,10 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
     negMask(i,2) = not(mask(i,2))
   end do
 
+  print*,'in selection '
   do i = 1, N_sel
+!    call debug_det(det(1,1,i),N_int)
+!    print*,i,dabs(psi_selectors_coef_transp_tc(1,2,i) * psi_selectors_coef_transp_tc(1,1,i))
     if(interesting(i) < 0) then
       stop 'prefetch interesting(i) and det(i)'
     endif
diff --git a/src/fci_tc_bi/generators.irp.f b/src/fci_tc_bi/generators.irp.f
index 55c0cbb9..250a1f71 100644
--- a/src/fci_tc_bi/generators.irp.f
+++ b/src/fci_tc_bi/generators.irp.f
@@ -31,6 +31,14 @@ END_PROVIDER
  END_DOC
  psi_det_generators(1:N_int,1:2,1:N_det) = psi_det_sorted_tc(1:N_int,1:2,1:N_det)
  psi_coef_generators(1:N_det,1:N_states) = psi_coef_sorted_tc(1:N_det,1:N_states)
+ integer :: i
+! print*,'generators '
+ do i = 1, N_det
+      if(N_det.ne.1)then
+      print*,'writing generators'
+      write(33,*) psi_det_generators(1,1,i), psi_det_generators(1,2,i)
+      endif
+ enddo
 
 END_PROVIDER
 
diff --git a/src/fci_tc_bi/selectors.irp.f b/src/fci_tc_bi/selectors.irp.f
index 3830927b..3c12bb07 100644
--- a/src/fci_tc_bi/selectors.irp.f
+++ b/src/fci_tc_bi/selectors.irp.f
@@ -43,15 +43,27 @@ END_PROVIDER
     do k=1,N_int
       psi_selectors(k,1,i) = psi_det_sorted_tc(k,1,i)
       psi_selectors(k,2,i) = psi_det_sorted_tc(k,2,i)
+!      psi_selectors(k,2,i) = psi_det(k,2,i)
+!      psi_selectors(k,2,i) = psi_det(k,2,i)
     enddo
   enddo
+ print*,'selectors '
   do k=1,N_states
     do i=1,N_det_selectors
       psi_selectors_coef(i,k) = psi_coef_sorted_tc_gen(i,k)
+!      psi_selectors_coef_tc(i,1,k) = psi_l_coef_bi_ortho(i,k)
+!      psi_selectors_coef_tc(i,2,k) = psi_r_coef_bi_ortho(i,k)
       psi_selectors_coef_tc(i,1,k) = psi_l_coef_sorted_bi_ortho(i,k)
       psi_selectors_coef_tc(i,2,k) = psi_r_coef_sorted_bi_ortho(i,k)
+!      call debug_det(psi_selectors(1,1,i),N_int)
+      if(N_det.ne.1)then
+       print*,'writing selectors'
+       write(34,*)psi_selectors(1,1,i),psi_selectors(1,2,i)
+       write(40,'(F10.7)')dabs(psi_selectors_coef_tc(i,1,1) * psi_selectors_coef_tc(i,2,1))
+      endif
 !      psi_selectors_coef_tc(i,1,k) = 1.d0
 !      psi_selectors_coef_tc(i,2,k) = 1.d0
+       
     enddo
   enddo
 
@@ -71,6 +83,9 @@ END_PROVIDER
       psi_selectors_coef_transp_tc(k,1,i) = psi_selectors_coef_tc(i,1,k)
       psi_selectors_coef_transp_tc(k,2,i) = psi_selectors_coef_tc(i,2,k)
     enddo
+    if(N_det.ne.1)then
+    write(41,'(F10.7)')dabs(psi_selectors_coef_transp_tc(1,1,i)*psi_selectors_coef_transp_tc(1,2,i))
+    endif
   enddo
 END_PROVIDER
 
diff --git a/src/tc_bi_ortho/psi_det_tc_sorted.irp.f b/src/tc_bi_ortho/psi_det_tc_sorted.irp.f
index 42617557..e8477dec 100644
--- a/src/tc_bi_ortho/psi_det_tc_sorted.irp.f
+++ b/src/tc_bi_ortho/psi_det_tc_sorted.irp.f
@@ -34,13 +34,19 @@ END_PROVIDER
    END_DOC
    integer                        :: i,j,k
    integer, allocatable           :: iorder(:)
+   print *, 'providing psi_det_sorted_tc'
    allocate ( iorder(N_det) )
+   print*,'before '
    do i=1,N_det
      psi_average_norm_contrib_sorted_tc(i) = -psi_average_norm_contrib_tc(i)
      iorder(i) = i
+     print*,i,iorder(i),psi_average_norm_contrib_sorted_tc(i)
    enddo
-!   call dsort(psi_average_norm_contrib_sorted_tc,iorder,N_det)
+   call dsort(psi_average_norm_contrib_sorted_tc,iorder,N_det)
+   print*,'after '
    do i=1,N_det
+!     iorder(i) = i
+     print*,i,iorder(i),psi_average_norm_contrib_sorted_tc(i)
      do j=1,N_int
        psi_det_sorted_tc(j,1,i) = psi_det(j,1,iorder(i))
        psi_det_sorted_tc(j,2,i) = psi_det(j,2,iorder(i))
@@ -67,6 +73,23 @@ END_PROVIDER
    psi_det_sorted_tc_order(N_det+1:psi_det_size) = 0
 
    deallocate(iorder)
+  logical :: pouet
+ pouet  = .true.
+ do i = 1, N_det
+  if(psi_average_norm_contrib_sorted_tc(i) == 0.d0)then
+   pouet = .False.
+   exit
+  endif
+ enddo
+ 
+ if(pouet.and.N_det.ne.1)then
+ print*,'writing sorted'
+ do i = 1, N_det
+!   call debug_det(psi_det_sorted_tc(1,1,i),N_int)
+   print*,i,psi_average_norm_contrib_sorted_tc(i)
+   write(35,*)psi_det_sorted_tc(1,1,i),psi_det_sorted_tc(1,2,i)
+ enddo
+ endif
 
 END_PROVIDER
 
@@ -84,6 +107,21 @@ END_PROVIDER
     psi_r_coef_sorted_bi_ortho(i,1) = psi_r_coef_bi_ortho(psi_det_sorted_tc_order(i),1)
     psi_l_coef_sorted_bi_ortho(i,1) = psi_l_coef_bi_ortho(psi_det_sorted_tc_order(i),1)
    enddo
+   logical :: pouet
+   pouet  = .true.
+   do i = 1, N_det
+    if(psi_l_coef_sorted_bi_ortho(i,1) == 0.d0)then
+     pouet = .False.
+     exit
+    endif
+   enddo
+   if(pouet.and.N_det.ne.1)then
+   print*,'psi_r_coef_sorted_bi_ortho'
+    do i = 1, N_det
+     print*,psi_r_coef_bi_ortho(psi_det_sorted_tc_order(i),1)
+     write(42,'(F10.7)')dabs(psi_r_coef_sorted_bi_ortho(i,1)*psi_l_coef_sorted_bi_ortho(i,1))
+    enddo
+   endif
 
 END_PROVIDER
 
diff --git a/src/tc_bi_ortho/tc_h_eigvectors.irp.f b/src/tc_bi_ortho/tc_h_eigvectors.irp.f
index d39b7a29..c66ff036 100644
--- a/src/tc_bi_ortho/tc_h_eigvectors.irp.f
+++ b/src/tc_bi_ortho/tc_h_eigvectors.irp.f
@@ -133,7 +133,10 @@
   call bi_normalize(leigvec_tc_bi_orth,reigvec_tc_bi_orth,size(reigvec_tc_bi_orth,1),N_det,N_states)
    print*,'leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1) = ',leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1)
    norm_ground_left_right_bi_orth = 0.d0
+   print*,'after diago'
    do j = 1, N_det
+    call debug_det(psi_det(1,1,j),N_int)
+    print*,j,dabs(leigvec_tc_bi_orth(j,1) * reigvec_tc_bi_orth(j,1))
     norm_ground_left_right_bi_orth += leigvec_tc_bi_orth(j,1) * reigvec_tc_bi_orth(j,1)
    enddo
    print*,'norm l/r = ',norm_ground_left_right_bi_orth

From 22fb8c17e21f986c7e5b18faf01ed206fa6ea73d Mon Sep 17 00:00:00 2001
From: eginer <giner.emmanuel@gmail.com>
Date: Thu, 16 Mar 2023 14:00:21 +0100
Subject: [PATCH 3/3] fixed the bug of misalignement between coefs and
 determinants in fci_tc_bi_ortho

---
 .../pt2_stoch_routines.irp.f                  |   1 -
 src/cipsi_tc_bi_ortho/selection.irp.f         |  29 +++--
 src/cipsi_tc_bi_ortho/selection_buffer.irp.f  |  12 +-
 src/cipsi_tc_bi_ortho/slave_cipsi.irp.f       |   2 -
 src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f  |   1 +
 src/fci_tc_bi/generators.irp.f                |   8 --
 src/fci_tc_bi/selectors.irp.f                 |  47 +-------
 src/tc_bi_ortho/psi_det_tc_sorted.irp.f       | 104 ++++++++----------
 src/tc_bi_ortho/psi_r_l_prov.irp.f            |  26 ++---
 src/tc_bi_ortho/tc_bi_ortho.irp.f             |  26 ++---
 src/tc_bi_ortho/tc_h_eigvectors.irp.f         |  25 ++++-
 src/tc_keywords/EZFIO.cfg                     |   8 ++
 12 files changed, 124 insertions(+), 165 deletions(-)

diff --git a/src/cipsi_tc_bi_ortho/pt2_stoch_routines.irp.f b/src/cipsi_tc_bi_ortho/pt2_stoch_routines.irp.f
index 027b74c5..284b2bc8 100644
--- a/src/cipsi_tc_bi_ortho/pt2_stoch_routines.irp.f
+++ b/src/cipsi_tc_bi_ortho/pt2_stoch_routines.irp.f
@@ -134,7 +134,6 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
   PROVIDE psi_det_hii selection_weight pseudo_sym
   PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max
   PROVIDE excitation_beta_max  excitation_alpha_max excitation_max
-  PROVIDE psi_selectors_rcoef_bi_orth_transp psi_selectors_lcoef_bi_orth_transp
 
   if (h0_type == 'CFG') then
     PROVIDE psi_configuration_hii det_to_configuration
diff --git a/src/cipsi_tc_bi_ortho/selection.irp.f b/src/cipsi_tc_bi_ortho/selection.irp.f
index 0c3f0451..4c271a4b 100644
--- a/src/cipsi_tc_bi_ortho/selection.irp.f
+++ b/src/cipsi_tc_bi_ortho/selection.irp.f
@@ -181,7 +181,6 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock
   PROVIDE psi_bilinear_matrix_rows psi_det_sorted_tc_order psi_bilinear_matrix_order
   PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
   PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp_tc
-  PROVIDE psi_selectors_rcoef_bi_orth_transp psi_selectors_lcoef_bi_orth_transp
 
   PROVIDE banned_excitation
 
@@ -616,7 +615,6 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
 
 
   PROVIDE psi_selectors_coef_transp_tc psi_det_sorted_tc
-  PROVIDE psi_selectors_rcoef_bi_orth_transp psi_selectors_lcoef_bi_orth_transp
 
 
   mat   = 0d0
@@ -628,7 +626,7 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
     negMask(i,2) = not(mask(i,2))
   end do
 
-  print*,'in selection '
+!  print*,'in selection '
   do i = 1, N_sel
 !    call debug_det(det(1,1,i),N_int)
 !    print*,i,dabs(psi_selectors_coef_transp_tc(1,2,i) * psi_selectors_coef_transp_tc(1,1,i))
@@ -677,9 +675,6 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
           perMask(j,1) = iand(mask(j,1), not(det(j,1,i)))
           perMask(j,2) = iand(mask(j,2), not(det(j,2,i)))
         end do
-!        call get_d3_h  ( det(1,1,i), bannedOrb, banned, mat         , mask, p, sp, psi_selectors_coef_transp_tc (1, interesting(i)) )
-!        call get_d3_htc( det(1,1,i), bannedOrb, banned, mat_r, mat_l, mask, p, sp, psi_selectors_rcoef_bi_orth_transp(1, interesting(i)) &
-!                       , psi_selectors_lcoef_bi_orth_transp(1, interesting(i)) )
 
         call bitstring_to_list_in_selection(perMask(1,1), h(1,1), h(0,1), N_int)
         call bitstring_to_list_in_selection(perMask(1,2), h(1,2), h(0,2), N_int)
@@ -921,15 +916,26 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
           endif
           val = 4.d0 * psi_h_alpha * alpha_h_psi 
           tmp = dsqrt(delta_E * delta_E + val)
-          if (delta_E < 0.d0) then
-              tmp = -tmp
-          endif
-          e_pert(istate) = 0.5d0 * (tmp - delta_E)
+!          if (delta_E < 0.d0) then
+!              tmp = -tmp
+!          endif
+          e_pert(istate) = 0.25 * val / delta_E
+!          e_pert(istate) = 0.5d0 * (tmp - delta_E)
           if(dsqrt(dabs(tmp)).gt.1.d-4.and.dabs(alpha_h_psi).gt.1.d-4)then
-           coef(istate)   = e_pert(istate) / alpha_h_psi
+           coef(istate)   = e_pert(istate) / psi_h_alpha 
           else
            coef(istate)   = alpha_h_psi / delta_E 
           endif
+
+          if(selection_tc == 1)then
+           if(e_pert(istate).lt.0.d0)then 
+            e_pert(istate)=0.d0
+           else 
+            e_pert(istate)=-e_pert(istate)
+           endif
+          else if(selection_tc == -1)then
+           if(e_pert(istate).gt.0.d0)e_pert(istate)=0.d0
+          endif
           
 !         if(selection_tc     ==  1 )then
 !          if(e_pert(istate).lt.0.d0)then
@@ -943,7 +949,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
         enddo
 
 
-      do_diag = sum(dabs(coef)) > 0.001d0 .and. N_states > 1
 
       do istate = 1, N_states
 
diff --git a/src/cipsi_tc_bi_ortho/selection_buffer.irp.f b/src/cipsi_tc_bi_ortho/selection_buffer.irp.f
index 10132086..0bd51464 100644
--- a/src/cipsi_tc_bi_ortho/selection_buffer.irp.f
+++ b/src/cipsi_tc_bi_ortho/selection_buffer.irp.f
@@ -125,7 +125,11 @@ subroutine merge_selection_buffers(b1, b2)
   enddo
   b2%det => detmp
   b2%val => val
-  b2%mini = min(b2%mini,b2%val(b2%N))
+!  if(selection_tc == 1)then
+!   b2%mini = max(b2%mini,b2%val(b2%N))
+!  else
+   b2%mini = min(b2%mini,b2%val(b2%N))
+!  endif
   b2%cur = nmwen
 end
 
@@ -157,7 +161,11 @@ subroutine sort_selection_buffer(b)
   end do
   deallocate(b%det,iorder)
   b%det => detmp
-  b%mini = min(b%mini,b%val(b%N))
+!  if(selection_tc == 1)then
+!   b%mini = max(b%mini,b%val(b%N))
+!  else
+   b%mini = min(b%mini,b%val(b%N))
+!  endif
   b%cur = nmwen
 end subroutine
 
diff --git a/src/cipsi_tc_bi_ortho/slave_cipsi.irp.f b/src/cipsi_tc_bi_ortho/slave_cipsi.irp.f
index c3a49280..6343bf8b 100644
--- a/src/cipsi_tc_bi_ortho/slave_cipsi.irp.f
+++ b/src/cipsi_tc_bi_ortho/slave_cipsi.irp.f
@@ -17,7 +17,6 @@ end
 
 subroutine provide_everything
   PROVIDE H_apply_buffer_allocated mo_two_e_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context N_states_diag
-  PROVIDE psi_selectors_rcoef_bi_orth_transp psi_selectors_lcoef_bi_orth_transp
 
   PROVIDE pt2_e0_denominator mo_num N_int ci_energy mpi_master zmq_state zmq_context
   PROVIDE psi_det psi_coef threshold_generators state_average_weight
@@ -312,7 +311,6 @@ subroutine run_slave_main
           PROVIDE psi_bilinear_matrix_rows psi_det_sorted_tc_order psi_bilinear_matrix_order
           PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
           PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted_tc
-          PROVIDE psi_selectors_rcoef_bi_orth_transp psi_selectors_lcoef_bi_orth_transp
 
           PROVIDE psi_det_hii selection_weight pseudo_sym pt2_min_parallel_tasks
 
diff --git a/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f b/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f
index e67287d3..e7ee4be9 100644
--- a/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f
+++ b/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f
@@ -108,6 +108,7 @@ subroutine run_stochastic_cipsi
     ept2(N_iter-1) = E_tc + nuclear_repulsion + (pt2_data % pt2(1))/norm
     pt1(N_iter-1) = dsqrt(pt2_data % overlap(1,1))
     call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2)
+!    stop
     if (qp_stop()) exit
   enddo
 !  print*,'data to extrapolate '
diff --git a/src/fci_tc_bi/generators.irp.f b/src/fci_tc_bi/generators.irp.f
index 250a1f71..55c0cbb9 100644
--- a/src/fci_tc_bi/generators.irp.f
+++ b/src/fci_tc_bi/generators.irp.f
@@ -31,14 +31,6 @@ END_PROVIDER
  END_DOC
  psi_det_generators(1:N_int,1:2,1:N_det) = psi_det_sorted_tc(1:N_int,1:2,1:N_det)
  psi_coef_generators(1:N_det,1:N_states) = psi_coef_sorted_tc(1:N_det,1:N_states)
- integer :: i
-! print*,'generators '
- do i = 1, N_det
-      if(N_det.ne.1)then
-      print*,'writing generators'
-      write(33,*) psi_det_generators(1,1,i), psi_det_generators(1,2,i)
-      endif
- enddo
 
 END_PROVIDER
 
diff --git a/src/fci_tc_bi/selectors.irp.f b/src/fci_tc_bi/selectors.irp.f
index 3c12bb07..4d3de7d0 100644
--- a/src/fci_tc_bi/selectors.irp.f
+++ b/src/fci_tc_bi/selectors.irp.f
@@ -18,15 +18,6 @@ BEGIN_PROVIDER [ integer, N_det_selectors]
   double precision               :: norm, norm_max
   call write_time(6)
   N_det_selectors = N_det
-!  norm = 1.d0
-!  do i=1,N_det
-!    norm = norm - psi_average_norm_contrib_tc(i)
-!    if (norm - 1.d-10 < 1.d0 - threshold_selectors) then
-!      N_det_selectors = i
-!      exit
-!    endif
-!  enddo
-  N_det_selectors = max(N_det_selectors,N_det_generators)
   call write_int(6,N_det_selectors,'Number of selectors')
 END_PROVIDER
 
@@ -43,27 +34,13 @@ END_PROVIDER
     do k=1,N_int
       psi_selectors(k,1,i) = psi_det_sorted_tc(k,1,i)
       psi_selectors(k,2,i) = psi_det_sorted_tc(k,2,i)
-!      psi_selectors(k,2,i) = psi_det(k,2,i)
-!      psi_selectors(k,2,i) = psi_det(k,2,i)
     enddo
   enddo
- print*,'selectors '
   do k=1,N_states
     do i=1,N_det_selectors
-      psi_selectors_coef(i,k) = psi_coef_sorted_tc_gen(i,k)
-!      psi_selectors_coef_tc(i,1,k) = psi_l_coef_bi_ortho(i,k)
-!      psi_selectors_coef_tc(i,2,k) = psi_r_coef_bi_ortho(i,k)
+      psi_selectors_coef(i,k)      = psi_coef_sorted_tc_gen(i,k)
       psi_selectors_coef_tc(i,1,k) = psi_l_coef_sorted_bi_ortho(i,k)
       psi_selectors_coef_tc(i,2,k) = psi_r_coef_sorted_bi_ortho(i,k)
-!      call debug_det(psi_selectors(1,1,i),N_int)
-      if(N_det.ne.1)then
-       print*,'writing selectors'
-       write(34,*)psi_selectors(1,1,i),psi_selectors(1,2,i)
-       write(40,'(F10.7)')dabs(psi_selectors_coef_tc(i,1,1) * psi_selectors_coef_tc(i,2,1))
-      endif
-!      psi_selectors_coef_tc(i,1,k) = 1.d0
-!      psi_selectors_coef_tc(i,2,k) = 1.d0
-       
     enddo
   enddo
 
@@ -83,31 +60,9 @@ END_PROVIDER
       psi_selectors_coef_transp_tc(k,1,i) = psi_selectors_coef_tc(i,1,k)
       psi_selectors_coef_transp_tc(k,2,i) = psi_selectors_coef_tc(i,2,k)
     enddo
-    if(N_det.ne.1)then
-    write(41,'(F10.7)')dabs(psi_selectors_coef_transp_tc(1,1,i)*psi_selectors_coef_transp_tc(1,2,i))
-    endif
   enddo
 END_PROVIDER
 
- BEGIN_PROVIDER [ double precision, psi_selectors_rcoef_bi_orth_transp, (N_states, psi_det_size) ]
-&BEGIN_PROVIDER [ double precision, psi_selectors_lcoef_bi_orth_transp, (N_states, psi_det_size) ]
-
-  implicit none
-  integer :: i, k
-
-  psi_selectors_rcoef_bi_orth_transp = 0.d0
-  psi_selectors_lcoef_bi_orth_transp = 0.d0
-
-  print*,'N_det,N_det_selectors',N_det,N_det_selectors
-  do i = 1, N_det_selectors
-    do k = 1, N_states
-      psi_selectors_rcoef_bi_orth_transp(k,i) = psi_r_coef_sorted_bi_ortho(i,k)
-      psi_selectors_lcoef_bi_orth_transp(k,i) = psi_l_coef_sorted_bi_ortho(i,k)
-    enddo
-  enddo
-
-END_PROVIDER
-
 BEGIN_PROVIDER [ integer, psi_selectors_size ]
  implicit none
  psi_selectors_size = psi_det_size
diff --git a/src/tc_bi_ortho/psi_det_tc_sorted.irp.f b/src/tc_bi_ortho/psi_det_tc_sorted.irp.f
index e8477dec..2d2111d6 100644
--- a/src/tc_bi_ortho/psi_det_tc_sorted.irp.f
+++ b/src/tc_bi_ortho/psi_det_tc_sorted.irp.f
@@ -11,7 +11,7 @@ BEGIN_PROVIDER [ double precision, psi_average_norm_contrib_tc, (psi_det_size) ]
   psi_average_norm_contrib_tc(:) = 0.d0
   do k=1,N_states
     do i=1,N_det
-      psi_average_norm_contrib_tc(i) = psi_average_norm_contrib_tc(i) +    &
+      psi_average_norm_contrib_tc(i) = & 
           dabs(psi_l_coef_bi_ortho(i,k)*psi_r_coef_bi_ortho(i,k))*state_average_weight(k)
     enddo
   enddo
@@ -26,39 +26,54 @@ END_PROVIDER
 &BEGIN_PROVIDER [ double precision, psi_coef_sorted_tc, (psi_det_size,N_states) ]
 &BEGIN_PROVIDER [ double precision, psi_average_norm_contrib_sorted_tc, (psi_det_size) ]
 &BEGIN_PROVIDER [ integer, psi_det_sorted_tc_order, (psi_det_size) ]
+&BEGIN_PROVIDER [double precision, psi_r_coef_sorted_bi_ortho, (psi_det_size, N_states)]
+&BEGIN_PROVIDER [double precision, psi_l_coef_sorted_bi_ortho, (psi_det_size, N_states)]
    implicit none
    BEGIN_DOC
    ! Wave function sorted by determinants contribution to the norm (state-averaged)
    !
    ! psi_det_sorted_tc_order(i) -> k : index in psi_det
+   !
+   ! psi_r_coef_sorted_bi_ortho : right coefficients corresponding to psi_det_sorted_tc
+   !
+   ! psi_l_coef_sorted_bi_ortho : left  coefficients corresponding to psi_det_sorted_tc
    END_DOC
    integer                        :: i,j,k
    integer, allocatable           :: iorder(:)
-   print *, 'providing psi_det_sorted_tc'
    allocate ( iorder(N_det) )
-   print*,'before '
+!   print*,'before = '
    do i=1,N_det
-     psi_average_norm_contrib_sorted_tc(i) = -psi_average_norm_contrib_tc(i)
      iorder(i) = i
-     print*,i,iorder(i),psi_average_norm_contrib_sorted_tc(i)
+     psi_average_norm_contrib_sorted_tc(i) = -psi_average_norm_contrib_tc(i)
+!    print*,'------------'
+!     call debug_det(psi_det(1,1,i),N_int)
+!     print*,i,psi_average_norm_contrib_tc(i)
+!    print*,i,psi_l_coef_bi_ortho(iorder(i),1:N_states),psi_r_coef_bi_ortho(iorder(i),1:N_states)
+!    print*,'------------'
    enddo
    call dsort(psi_average_norm_contrib_sorted_tc,iorder,N_det)
-   print*,'after '
+!   print*,'after = '
    do i=1,N_det
-!     iorder(i) = i
-     print*,i,iorder(i),psi_average_norm_contrib_sorted_tc(i)
+     psi_average_norm_contrib_sorted_tc(i) = -psi_average_norm_contrib_sorted_tc(i)
      do j=1,N_int
        psi_det_sorted_tc(j,1,i) = psi_det(j,1,iorder(i))
        psi_det_sorted_tc(j,2,i) = psi_det(j,2,iorder(i))
      enddo
-    psi_average_norm_contrib_sorted_tc(i) = -psi_average_norm_contrib_sorted_tc(i)
-    psi_det_sorted_tc_order(iorder(i)) = i
+     psi_det_sorted_tc_order(iorder(i)) = i
+!     if(iorder(i).ne.i)then
+!      print*,'changed the order for ',i,iorder(i)
+!     endif
+!      print*,'------------'
+!      call debug_det(psi_det_sorted_tc(1,1,i),N_int)
+!      print*,i,psi_average_norm_contrib_tc(i)
+!      print*,i,psi_l_coef_bi_ortho(iorder(i),1:N_states),psi_r_coef_bi_ortho(iorder(i),1:N_states)
+!      print*,'------------'
    enddo
    double precision :: accu
    do k=1,N_states
     accu = 0.d0
     do i=1,N_det
-     psi_coef_sorted_tc(i,k) = dsqrt(dabs(psi_l_coef_bi_ortho(iorder(i),k)*psi_r_coef_bi_ortho(iorder(i),k)))
+     psi_coef_sorted_tc(i,k) = dsqrt(psi_average_norm_contrib_sorted_tc(i))
      accu += psi_coef_sorted_tc(i,k)**2
     enddo
     accu = 1.d0/dsqrt(accu)
@@ -72,60 +87,33 @@ END_PROVIDER
    psi_average_norm_contrib_sorted_tc(N_det+1:psi_det_size) = 0.d0
    psi_det_sorted_tc_order(N_det+1:psi_det_size) = 0
 
-   deallocate(iorder)
-  logical :: pouet
- pouet  = .true.
- do i = 1, N_det
-  if(psi_average_norm_contrib_sorted_tc(i) == 0.d0)then
-   pouet = .False.
-   exit
-  endif
- enddo
- 
- if(pouet.and.N_det.ne.1)then
- print*,'writing sorted'
- do i = 1, N_det
-!   call debug_det(psi_det_sorted_tc(1,1,i),N_int)
-   print*,i,psi_average_norm_contrib_sorted_tc(i)
-   write(35,*)psi_det_sorted_tc(1,1,i),psi_det_sorted_tc(1,2,i)
- enddo
- endif
-
-END_PROVIDER
-
- BEGIN_PROVIDER [double precision, psi_r_coef_sorted_bi_ortho, (psi_det_size, N_states)]
-&BEGIN_PROVIDER [double precision, psi_l_coef_sorted_bi_ortho, (psi_det_size, N_states)]
- BEGIN_DOC
- ! psi_r_coef_sorted_bi_ortho : right coefficients corresponding to psi_det_sorted_tc
- ! psi_l_coef_sorted_bi_ortho : left  coefficients corresponding to psi_det_sorted_tc
- END_DOC
-   implicit none
-   integer                       :: i, j, k
    psi_r_coef_sorted_bi_ortho = 0.d0
    psi_l_coef_sorted_bi_ortho = 0.d0
    do i = 1, N_det
-    psi_r_coef_sorted_bi_ortho(i,1) = psi_r_coef_bi_ortho(psi_det_sorted_tc_order(i),1)
-    psi_l_coef_sorted_bi_ortho(i,1) = psi_l_coef_bi_ortho(psi_det_sorted_tc_order(i),1)
+    psi_r_coef_sorted_bi_ortho(i,1:N_states) = psi_r_coef_bi_ortho(iorder(i),1:N_states)
+    psi_l_coef_sorted_bi_ortho(i,1:N_states) = psi_l_coef_bi_ortho(iorder(i),1:N_states)
    enddo
-   logical :: pouet
-   pouet  = .true.
-   do i = 1, N_det
-    if(psi_l_coef_sorted_bi_ortho(i,1) == 0.d0)then
-     pouet = .False.
-     exit
-    endif
-   enddo
-   if(pouet.and.N_det.ne.1)then
-   print*,'psi_r_coef_sorted_bi_ortho'
-    do i = 1, N_det
-     print*,psi_r_coef_bi_ortho(psi_det_sorted_tc_order(i),1)
-     write(42,'(F10.7)')dabs(psi_r_coef_sorted_bi_ortho(i,1)*psi_l_coef_sorted_bi_ortho(i,1))
-    enddo
-   endif
+
+   deallocate(iorder)
+!  logical :: pouet
+! pouet  = .true.
+! do i = 1, N_det
+!  if(psi_average_norm_contrib_sorted_tc(i) == 0.d0)then
+!   pouet = .False.
+!   exit
+!  endif
+! enddo
+! 
+! if(pouet.and.N_det.ne.1)then
+! print*,'writing sorted'
+! do i = 1, N_det
+!   print*,i,psi_average_norm_contrib_sorted_tc(i)
+!   write(35,*)psi_det_sorted_tc(1,1,i),psi_det_sorted_tc(1,2,i)
+! enddo
+! endif
 
 END_PROVIDER
 
-
  BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_tc_bit, (N_int,2,psi_det_size) ]
 &BEGIN_PROVIDER [ double precision, psi_coef_sorted_tc_bit, (psi_det_size,N_states) ]
    implicit none
diff --git a/src/tc_bi_ortho/psi_r_l_prov.irp.f b/src/tc_bi_ortho/psi_r_l_prov.irp.f
index ac9b0e74..521acff5 100644
--- a/src/tc_bi_ortho/psi_r_l_prov.irp.f
+++ b/src/tc_bi_ortho/psi_r_l_prov.irp.f
@@ -136,15 +136,15 @@ BEGIN_PROVIDER [ double precision, psi_r_coef_bi_ortho, (psi_det_size,N_states)
 END_PROVIDER
 
 
-subroutine save_tc_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psilcoef,psircoef)
+subroutine save_tc_wavefunction_general(ndet,nstates,psidet,sze,dim_psicoef,psilcoef,psircoef)
   implicit none
   BEGIN_DOC
   !  Save the wave function into the |EZFIO| file
   END_DOC
   use bitmasks
   include 'constants.include.F'
-  integer, intent(in)            :: ndet,nstates,dim_psicoef
-  integer(bit_kind), intent(in)  :: psidet(N_int,2,ndet)
+  integer, intent(in)            :: ndet,nstates,dim_psicoef,sze
+  integer(bit_kind), intent(in)  :: psidet(N_int,2,sze)
   double precision, intent(in)   :: psilcoef(dim_psicoef,nstates)
   double precision, intent(in)   :: psircoef(dim_psicoef,nstates)
   integer*8, allocatable         :: psi_det_save(:,:,:)
@@ -188,23 +188,17 @@ subroutine save_tc_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psilcoef
     call ezfio_set_tc_bi_ortho_psi_r_coef_bi_ortho(psir_coef_save)
     deallocate (psil_coef_save,psir_coef_save)
 
-!    allocate (psi_coef_save(ndet_qp_edit,nstates))
-!    do k=1,nstates
-!      do i=1,ndet_qp_edit
-!        psi_coef_save(i,k) = psicoef(i,k)
-!      enddo
-!    enddo
-!
-!    call ezfio_set_determinants_psi_coef_qp_edit(psi_coef_save)
-!    deallocate (psi_coef_save)
-
     call write_int(6,ndet,'Saved determinantsi and psi_r/psi_l coef')
   endif
 end
 
 subroutine save_tc_bi_ortho_wavefunction
  implicit none
- call save_tc_wavefunction_general(N_det,N_states,psi_det,size(psi_l_coef_bi_ortho, 1),psi_l_coef_bi_ortho,psi_r_coef_bi_ortho)
+ if(save_sorted_tc_wf)then
+  call save_tc_wavefunction_general(N_det,N_states,psi_det_sorted_tc,size(psi_det_sorted_tc, 3),size(psi_l_coef_sorted_bi_ortho, 1),psi_l_coef_sorted_bi_ortho,psi_r_coef_sorted_bi_ortho)
+ else
+  call save_tc_wavefunction_general(N_det,N_states,psi_det,size(psi_det, 3), size(psi_l_coef_bi_ortho, 1),psi_l_coef_bi_ortho,psi_r_coef_bi_ortho)
+ endif
  call routine_save_right_bi_ortho
 end
 
@@ -214,9 +208,9 @@ subroutine routine_save_right_bi_ortho
  integer :: i
  allocate(coef_tmp(N_det, N_states))
  do i = 1, N_det
-  coef_tmp(i,1:N_states) = psi_r_coef_bi_ortho(i,1:N_states)
+  coef_tmp(i,1:N_states) = psi_r_coef_sorted_bi_ortho(i,1:N_states)
  enddo
- call save_wavefunction_general_unormalized(N_det,N_states,psi_det,size(coef_tmp,1),coef_tmp(1,1))
+ call save_wavefunction_general_unormalized(N_det,N_states,psi_det_sorted_tc,size(coef_tmp,1),coef_tmp(1,1))
 end                     
 
 subroutine routine_save_left_right_bi_ortho
diff --git a/src/tc_bi_ortho/tc_bi_ortho.irp.f b/src/tc_bi_ortho/tc_bi_ortho.irp.f
index cfa24f3b..2d51f6f0 100644
--- a/src/tc_bi_ortho/tc_bi_ortho.irp.f
+++ b/src/tc_bi_ortho/tc_bi_ortho.irp.f
@@ -11,7 +11,7 @@ program tc_bi_ortho
   touch read_wf
   touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
   call routine_diag
-! call test
+  call save_tc_bi_ortho_wavefunction
 end
 
 subroutine test
@@ -19,18 +19,19 @@ subroutine test
  integer :: i,j 
  double precision :: hmono,htwoe,hthree,htot
  use bitmasks
-
- print*,'test'
-! call htilde_mu_mat_bi_ortho(psi_det(1,1,1), psi_det(1,1,2), N_int, hmono, htwoe, hthree, htot)
- call double_htilde_mu_mat_bi_ortho(N_int,psi_det(1,1,1), psi_det(1,1,2), hmono, htwoe, htot)
- print*,hmono, htwoe,  htot
+ print*,'reading the wave function '
+ do i = 1, N_det
+  call debug_det(psi_det(1,1,i),N_int)
+  print*,i,psi_l_coef_bi_ortho(i,1)*psi_r_coef_bi_ortho(i,1)
+  print*,i,psi_l_coef_bi_ortho(i,1),psi_r_coef_bi_ortho(i,1)
+ enddo
 
 end
 
 subroutine routine_diag
  implicit none
 ! provide eigval_right_tc_bi_orth
-  provide overlap_bi_ortho
+!  provide overlap_bi_ortho
 !  provide htilde_matrix_elmt_bi_ortho
  integer ::i,j
  print*,'eigval_right_tc_bi_orth = ',eigval_right_tc_bi_orth(1)
@@ -46,16 +47,7 @@ subroutine routine_diag
  print*,'e_corr_double_bi_orth   = ',e_corr_double_bi_orth
  print*,'Left/right eigenvectors'
  do i = 1,N_det
-  write(*,'(I5,X,(100(F12.7,X)))')i,leigvec_tc_bi_orth(i,1),reigvec_tc_bi_orth(i,1)
+  write(*,'(I5,X,(100(F12.7,X)))')i,leigvec_tc_bi_orth(i,1),reigvec_tc_bi_orth(i,1),leigvec_tc_bi_orth(i,1)*reigvec_tc_bi_orth(i,1)
  enddo
- do j=1,N_states
-   do i=1,N_det
-     psi_l_coef_bi_ortho(i,j) = leigvec_tc_bi_orth(i,j)
-     psi_r_coef_bi_ortho(i,j) = reigvec_tc_bi_orth(i,j)
-   enddo
- enddo
- SOFT_TOUCH psi_l_coef_bi_ortho psi_r_coef_bi_ortho
- call save_tc_bi_ortho_wavefunction
-! call routine_save_left_right_bi_ortho
 end
 
diff --git a/src/tc_bi_ortho/tc_h_eigvectors.irp.f b/src/tc_bi_ortho/tc_h_eigvectors.irp.f
index c66ff036..1ccda822 100644
--- a/src/tc_bi_ortho/tc_h_eigvectors.irp.f
+++ b/src/tc_bi_ortho/tc_h_eigvectors.irp.f
@@ -12,6 +12,25 @@
  enddo
  END_PROVIDER
 
+subroutine diagonalize_CI_tc
+  implicit none
+  BEGIN_DOC
+!  Replace the coefficients of the |CI| states by the coefficients of the
+!  eigenstates of the |CI| matrix.
+  END_DOC
+  integer                        :: i,j
+  do j=1,N_states
+    do i=1,N_det
+      psi_l_coef_bi_ortho(i,j) = leigvec_tc_bi_orth(i,j)
+      psi_r_coef_bi_ortho(i,j) = reigvec_tc_bi_orth(i,j)
+    enddo
+  enddo
+!  psi_energy(1:N_states) = CI_electronic_energy(1:N_states)
+!  psi_s2(1:N_states) = CI_s2(1:N_states)
+
+  SOFT_TOUCH psi_l_coef_bi_ortho psi_r_coef_bi_ortho
+end
+
 
 
  BEGIN_PROVIDER [double precision, eigval_right_tc_bi_orth, (N_states)]
@@ -133,10 +152,10 @@
   call bi_normalize(leigvec_tc_bi_orth,reigvec_tc_bi_orth,size(reigvec_tc_bi_orth,1),N_det,N_states)
    print*,'leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1) = ',leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1)
    norm_ground_left_right_bi_orth = 0.d0
-   print*,'after diago'
+!   print*,'after diago'
    do j = 1, N_det
-    call debug_det(psi_det(1,1,j),N_int)
-    print*,j,dabs(leigvec_tc_bi_orth(j,1) * reigvec_tc_bi_orth(j,1))
+!    call debug_det(psi_det(1,1,j),N_int)
+!    print*,j,dabs(leigvec_tc_bi_orth(j,1) * reigvec_tc_bi_orth(j,1))
     norm_ground_left_right_bi_orth += leigvec_tc_bi_orth(j,1) * reigvec_tc_bi_orth(j,1)
    enddo
    print*,'norm l/r = ',norm_ground_left_right_bi_orth
diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg
index 8765cd6e..e397e700 100644
--- a/src/tc_keywords/EZFIO.cfg
+++ b/src/tc_keywords/EZFIO.cfg
@@ -207,3 +207,11 @@ type: logical
 doc: If |true|, only the right part of WF is used to compute spin dens
 interface: ezfio,provider,ocaml
 default: False
+
+[save_sorted_tc_wf]
+type: logical
+doc: If |true|, save the bi-ortho wave functions in a sorted way
+interface: ezfio,provider,ocaml
+default: True
+
+