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 ee8abcec..b28c417f 100644
--- a/src/tc_bi_ortho/psi_r_l_prov.irp.f
+++ b/src/tc_bi_ortho/psi_r_l_prov.irp.f
@@ -192,33 +192,51 @@ subroutine save_tc_wavefunction_general(ndet, nstates, psidet, sze, dim_psicoef,
   endif
 end
 
+! ---
+
 subroutine save_tc_bi_ortho_wavefunction()
+
   implicit none
-  if(save_sorted_tc_wf)then
+
+  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)
+    call routine_save_right_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 )
+    call routine_save_right_bi_ortho()
+
   endif
-  call routine_save_right_bi_ortho()
+
 end
 
-subroutine routine_save_right_bi_ortho
+! ---
+
+subroutine routine_save_right_sorted_bi_ortho()
+
   implicit none
+  integer                       :: i
   double precision, allocatable :: coef_tmp(:,:)
-  integer :: i
+
   allocate(coef_tmp(N_det, N_states))
   do i = 1, N_det
     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_sorted_tc, size(coef_tmp, 1), coef_tmp(1,1))
-end                     
+  deallocate(coef_tmp)
+
+end
+
+subroutine routine_save_left_right_sorted_bi_ortho()
 
-subroutine routine_save_left_right_bi_ortho
   implicit none
+  integer                       :: i, n_states_tmp
   double precision, allocatable :: coef_tmp(:,:)
-  integer :: i,n_states_tmp
+
   n_states_tmp = 2
   allocate(coef_tmp(N_det, n_states_tmp))
   do i = 1, N_det
@@ -226,5 +244,26 @@ subroutine routine_save_left_right_bi_ortho
     coef_tmp(i,2) = psi_l_coef_bi_ortho(i,1)
   enddo
   call save_wavefunction_general_unormalized(N_det, n_states_tmp, psi_det, size(coef_tmp, 1), coef_tmp(1,1))
+  deallocate(coef_tmp)
 end
 
+! ---
+
+subroutine routine_save_right_bi_ortho()
+
+  implicit none
+  integer                       :: i
+  double precision, allocatable :: coef_tmp(:,:)
+
+  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)
+  enddo
+  call save_wavefunction_general_unormalized(N_det, N_states, psi_det, size(coef_tmp, 1), coef_tmp(1,1))
+  deallocate(coef_tmp)
+
+end
+
+! ---
+
+
diff --git a/src/tc_bi_ortho/save_lr_bi_ortho_states.irp.f b/src/tc_bi_ortho/save_lr_bi_ortho_states.irp.f
deleted file mode 100644
index 5eb3c069..00000000
--- a/src/tc_bi_ortho/save_lr_bi_ortho_states.irp.f
+++ /dev/null
@@ -1,15 +0,0 @@
-program tc_bi_ortho
-  implicit none
-  BEGIN_DOC
-! TODO : Put the documentation of the program here
-  END_DOC
-  print *, 'Hello world'
-  my_grid_becke = .True.
-  my_n_pt_r_grid = 30
-  my_n_pt_a_grid = 50
-  read_wf = .True.
-  touch read_wf
-  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
-  call routine_save_left_right_bi_ortho
-! call test
-end
diff --git a/src/tc_bi_ortho/tc_h_eigvectors.irp.f b/src/tc_bi_ortho/tc_h_eigvectors.irp.f
index 3d452b68..a288810b 100644
--- a/src/tc_bi_ortho/tc_h_eigvectors.irp.f
+++ b/src/tc_bi_ortho/tc_h_eigvectors.irp.f
@@ -266,7 +266,7 @@ end
     converged = .False.
     i_it = 0
     do while (.not. converged)
-      call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_dagger_u_0_opt)
+      call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_u_0_opt)
       i_it += 1
       if(i_it .gt. 5) exit
     enddo
@@ -324,42 +324,54 @@ END_PROVIDER
 
 
 
-subroutine bi_normalize(u_l,u_r,n,ld,nstates)
+subroutine bi_normalize(u_l, u_r, n, ld, nstates)
+
+  BEGIN_DOC
   !!!! Normalization of the scalar product of the left/right eigenvectors
+  END_DOC
+
+  implicit none
+  integer,          intent(in)    :: n, ld, nstates
   double precision, intent(inout) :: u_l(ld,nstates), u_r(ld,nstates)
-  integer, intent(in) :: n,ld,nstates
-  integer :: i
-  double precision  :: accu, tmp 
+  integer                         :: i, j
+  double precision                :: accu, tmp
+
   do i = 1, nstates
-   !!!! Normalization of right eigenvectors |Phi>
-   accu = 0.d0
-   do j = 1, n
-    accu += u_r(j,i) * u_r(j,i)
-   enddo
-   accu = 1.d0/dsqrt(accu)
-   print*,'accu_r = ',accu
-   do j = 1, n
-    u_r(j,i) *= accu 
-   enddo
-   tmp = u_r(1,i) / dabs(u_r(1,i))
-   do j = 1, n
-    u_r(j,i) *= tmp
-   enddo
-   !!!! Adaptation of the norm of the left eigenvector such that <chi|Phi> = 1
-   accu = 0.d0
-   do j = 1, n
-    accu += u_l(j,i) * u_r(j,i)
-!    print*,j, u_l(j,i) , u_r(j,i)
-   enddo
-   if(accu.gt.0.d0)then
+
+    !!!! Normalization of right eigenvectors |Phi>
+    accu = 0.d0
+    do j = 1, n
+      accu += u_r(j,i) * u_r(j,i)
+    enddo
     accu = 1.d0/dsqrt(accu)
-   else
-    accu = 1.d0/dsqrt(-accu)
-   endif
-   tmp = (u_l(1,i) * u_r(1,i) )/dabs(u_l(1,i) * u_r(1,i))
-   do j = 1, n
-    u_l(j,i) *= accu * tmp
-    u_r(j,i) *= accu 
-   enddo
+    print*,'accu_r = ',accu
+    do j = 1, n
+      u_r(j,i) *= accu
+    enddo
+    tmp = u_r(1,i) / dabs(u_r(1,i))
+    do j = 1, n
+      u_r(j,i) *= tmp
+    enddo
+
+    !!!! Adaptation of the norm of the left eigenvector such that <chi|Phi> = 1
+    accu = 0.d0
+    do j = 1, n
+      accu += u_l(j,i) * u_r(j,i)
+      !print*,j, u_l(j,i) , u_r(j,i)
+    enddo
+    print*,'accu_lr = ', accu
+    if(accu.gt.0.d0)then
+      accu = 1.d0/dsqrt(accu)
+    else
+      accu = 1.d0/dsqrt(-accu)
+    endif
+    tmp = (u_l(1,i) * u_r(1,i) )/dabs(u_l(1,i) * u_r(1,i))
+    do j = 1, n
+      u_l(j,i) *= accu * tmp
+      u_r(j,i) *= accu
+    enddo
+
   enddo
+
 end
+