diff --git a/src/tc_bi_ortho/print_tc_energy.irp.f b/src/tc_bi_ortho/print_tc_energy.irp.f
new file mode 100644
index 00000000..c9f8cecb
--- /dev/null
+++ b/src/tc_bi_ortho/print_tc_energy.irp.f
@@ -0,0 +1,47 @@
+program print_tc_energy
+  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 write_tc_energy
+end
+
+subroutine write_tc_energy()
+
+  implicit none
+  integer          :: i, j, k
+  double precision :: hmono, htwoe, htot
+  double precision :: E_TC, O_TC
+
+  do k = 1, n_states
+
+    E_TC = 0.d0
+    do i = 1, N_det
+      do j = 1, N_det
+        call hmat_bi_ortho(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, htot)
+        E_TC = E_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * htot
+        !E_TC = E_TC + leigvec_tc_bi_orth(i,k) * reigvec_tc_bi_orth(j,k) * htot
+      enddo
+    enddo
+
+    O_TC = 0.d0
+    do i = 1, N_det
+      O_TC = O_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(i,k)
+      !O_TC = O_TC + leigvec_tc_bi_orth(i,k) * reigvec_tc_bi_orth(i,k)
+    enddo
+
+    print *, ' state :', k
+    print *, " E_TC = ", E_TC
+    print *, " O_TC = ", O_TC
+
+  enddo
+
+end
+
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 521acff5..05d132d5 100644
--- a/src/tc_bi_ortho/psi_r_l_prov.irp.f
+++ b/src/tc_bi_ortho/psi_r_l_prov.irp.f
@@ -136,7 +136,7 @@ BEGIN_PROVIDER [ double precision, psi_r_coef_bi_ortho, (psi_det_size,N_states)
 END_PROVIDER
 
 
-subroutine save_tc_wavefunction_general(ndet,nstates,psidet,sze,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
@@ -195,9 +195,16 @@ end
 subroutine save_tc_bi_ortho_wavefunction
  implicit none
  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 save_tc_wavefunction_general( N_det, N_states, psi_det_sorted_tc, N_det &
+                                   , size(psi_l_coef_sorted_bi_ortho, 1), psi_l_coef_sorted_bi_ortho &
+                                   , psi_r_coef_sorted_bi_ortho )
+  !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)
+  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
diff --git a/src/tc_bi_ortho/tc_bi_ortho.irp.f b/src/tc_bi_ortho/tc_bi_ortho.irp.f
index c11164b0..ef2e5659 100644
--- a/src/tc_bi_ortho/tc_bi_ortho.irp.f
+++ b/src/tc_bi_ortho/tc_bi_ortho.irp.f
@@ -1,8 +1,14 @@
 program tc_bi_ortho
-  implicit none
+
   BEGIN_DOC
-! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together with the energy. Saves the left-right wave functions at the end. 
+  !
+  ! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together 
+  !        with the energy. Saves the left-right wave functions at the end. 
+  !
   END_DOC
+
+  implicit none
+
   print *, 'Hello world'
   my_grid_becke = .True.
   my_n_pt_r_grid = 30
@@ -15,6 +21,7 @@ program tc_bi_ortho
   print*, ' nb of det    = ', N_det
 
   call routine_diag()
+  call write_tc_energy()
   call save_tc_bi_ortho_wavefunction()
 end
 
@@ -82,3 +89,36 @@ end
 
 
 
+subroutine write_tc_energy()
+
+  implicit none
+  integer          :: i, j, k
+  double precision :: hmono, htwoe, htot
+  double precision :: E_TC, O_TC
+
+  do k = 1, n_states
+
+    E_TC = 0.d0
+    do i = 1, N_det
+      do j = 1, N_det
+        htot = htilde_matrix_elmt_bi_ortho(j,i)
+        !call hmat_bi_ortho(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, htot)
+        E_TC = E_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * htot
+        !E_TC = E_TC + leigvec_tc_bi_orth(i,k) * reigvec_tc_bi_orth(j,k) * htot
+      enddo
+    enddo
+
+    O_TC = 0.d0
+    do i = 1, N_det
+      O_TC = O_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(i,k)
+      !O_TC = O_TC + leigvec_tc_bi_orth(i,k) * reigvec_tc_bi_orth(i,k)
+    enddo
+
+    print *, ' state :', k
+    print *, " E_TC = ", E_TC
+    print *, " O_TC = ", O_TC
+
+  enddo
+
+end
+
diff --git a/src/tc_bi_ortho/tc_h_eigvectors.irp.f b/src/tc_bi_ortho/tc_h_eigvectors.irp.f
index dc010701..fa9d9f5c 100644
--- a/src/tc_bi_ortho/tc_h_eigvectors.irp.f
+++ b/src/tc_bi_ortho/tc_h_eigvectors.irp.f
@@ -44,7 +44,7 @@ end
   END_DOC
 
   implicit none
-  integer                       :: i, idx_dress, j, istate
+  integer                       :: i, idx_dress, j, istate, k
   logical                       :: converged, dagger
   integer                       :: n_real_tc_bi_orth_eigval_right,igood_r,igood_l
   integer,          allocatable :: iorder(:)
@@ -168,13 +168,39 @@ end
 
     deallocate(H_jj)
    endif
+
   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
-   do j = 1, N_det
-    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
+  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)
+  do i = 1, N_states
+    norm_ground_left_right_bi_orth = 0.d0
+    do j = 1, N_det
+      norm_ground_left_right_bi_orth += leigvec_tc_bi_orth(j,i) * reigvec_tc_bi_orth(j,i)
+    enddo
+    print*,'norm l/r = ',norm_ground_left_right_bi_orth
+  enddo
+
+  ! ---
+
+  double precision, allocatable :: buffer(:,:)
+  allocate(buffer(N_det,N_states))
+
+  do k = 1, N_states
+    do i = 1, N_det
+      buffer(i,k) = leigvec_tc_bi_orth(i,k)
+    enddo
+  enddo
+  call ezfio_set_tc_bi_ortho_psi_l_coef_bi_ortho(buffer)
+
+  do k = 1, N_states
+    do i = 1, N_det
+      buffer(i,k) = reigvec_tc_bi_orth(i,k)
+    enddo
+  enddo
+  call ezfio_set_tc_bi_ortho_psi_r_coef_bi_ortho(buffer)
+
+  deallocate(buffer)
+
+  ! ---
 
 END_PROVIDER