diff --git a/src/tc_bi_ortho/test_tc_two_rdm.irp.f b/src/tc_bi_ortho/test_tc_two_rdm.irp.f
new file mode 100644
index 00000000..ecdeef43
--- /dev/null
+++ b/src/tc_bi_ortho/test_tc_two_rdm.irp.f
@@ -0,0 +1,60 @@
+program test_tc_rdm
+
+  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. 
+  !
+  END_DOC
+
+  my_grid_becke = .True.
+  PROVIDE tc_grid1_a tc_grid1_r
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
+  read_wf = .True.
+  touch read_wf
+
+  print*, ' nb of states = ', N_states
+  print*, ' nb of det    = ', N_det
+
+  call test()
+
+end
+
+subroutine  test
+ implicit none
+ integer :: h1,p1,h2,p2,i,j,istate
+ double precision :: rdm, integral, accu,ref
+ double precision :: hmono, htwoe, hthree, htot
+ accu = 0.d0
+ do h1 = 1, mo_num
+  do p1 = 1, mo_num
+   do h2 = 1, mo_num
+    do p2 = 1, mo_num
+     integral = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
+     rdm = tc_two_rdm(p1,h1,p2,h2)
+     accu += integral * rdm
+    enddo
+   enddo
+  enddo
+ enddo
+ accu *= 0.5d0
+ print*,'accu = ',accu
+! print*,tc_two_rdm(1,1,1,1),mo_bi_ortho_tc_two_e(1,1,1,1)
+ ref = 0.d0
+ do i = 1, N_det
+  do j = 1, N_det
+!   if(i.ne.j)cycle
+   call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot)
+   do istate = 1,N_states
+!  print*,'i,j',i,j,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * state_average_weight(istate) * htwoe 
+!  print*,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) ,  htwoe 
+    ref += psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * state_average_weight(istate) * htwoe 
+   enddo
+  enddo
+ enddo
+ print*,' ref = ',ref 
+
+end
diff --git a/src/tc_bi_ortho/two_rdm_naive.irp.f b/src/tc_bi_ortho/two_rdm_naive.irp.f
new file mode 100644
index 00000000..9694c653
--- /dev/null
+++ b/src/tc_bi_ortho/two_rdm_naive.irp.f
@@ -0,0 +1,132 @@
+BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+ ! tc_two_rdm(p,s,q,r) = <Phi| a^dager_p
+ END_DOC
+ integer :: i,j,istate,m,mm,nn
+ integer                        :: exc(0:2,2,2)
+ double precision               :: phase
+ double precision               :: contrib
+ integer                        :: h1,p1,s1,h2,p2,s2,degree
+ integer, allocatable           :: occ(:,:)
+ integer                        :: n_occ_ab(2),other_spin(2)
+ other_spin(1) = 2
+ other_spin(2) = 1
+ allocate(occ(N_int*bit_kind_size,2))
+ tc_two_rdm = 0.d0
+
+ do i = 1, N_det ! psi_left 
+  do j = 1, N_det ! psi_right 
+   call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
+   if(degree.gt.2)cycle
+   if(degree.gt.0)then
+    ! get excitation operators: from psi_det(j) --> psi_det(i)
+    call get_excitation(psi_det(1,1,j),psi_det(1,1,i),exc,degree,phase,N_int)
+    call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+    contrib = psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * phase * state_average_weight(1)
+    do istate = 2, N_states
+     contrib += psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * phase * state_average_weight(istate)
+    enddo
+    if(degree == 2)then
+     call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+    else if(degree==1)then
+!     cycle
+     ! occupation of the determinant psi_det(j)
+     call bitstring_to_list_ab(psi_det(1,1,j), occ, n_occ_ab, N_int) 
+ 
+     ! run over the electrons of opposite spin than the excitation
+     s2 = other_spin(s1)
+     do mm = 1, n_occ_ab(s2) 
+      m = occ(mm,s2)
+      h2 = m 
+      p2 = m 
+      call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+     enddo
+     ! run over the electrons of same spin than the excitation
+     s2 = s1
+     do mm = 1, n_occ_ab(s2) 
+      m = occ(mm,s2)
+      h2 = m 
+      p2 = m 
+      if(h2.le.h1)cycle
+      call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+     enddo
+    endif
+   else if(degree == 0)then
+    contrib = psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) *  state_average_weight(1)
+!    print*,'contrib',contrib
+    do istate = 2, N_states
+     contrib += psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) *  state_average_weight(istate)
+    enddo
+    ! occupation of the determinant psi_det(j)
+    call bitstring_to_list_ab(psi_det(1,1,j), occ, n_occ_ab, N_int) 
+    s1 = 1 ! alpha electrons
+     do nn = 1, n_occ_ab(s1)
+      h1 = occ(nn,s1)
+      p1 = occ(nn,s1)
+      ! run over the couple of alpha-beta electrons 
+      s2 = other_spin(s1)
+      do mm = 1, n_occ_ab(s2) 
+       m = occ(mm,s2)
+       h2 = m 
+       p2 = m 
+       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+      enddo
+      ! run over the couple of alpha-alpha electrons 
+      s2 = s1
+      do mm = 1, n_occ_ab(s2) 
+       m = occ(mm,s2)
+       h2 = m 
+       p2 = m 
+       if(h2.le.h1)cycle
+       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+      enddo
+    enddo
+    s1 = 2
+     do nn = 1, n_occ_ab(s1)
+      h1 = occ(nn,s1)
+      p1 = occ(nn,s1)
+      ! run over the couple of beta-beta electrons 
+      s2 = s1
+      do mm = 1, n_occ_ab(s2) 
+       m = occ(mm,s2)
+       h2 = m 
+       p2 = m 
+       if(h2.le.h1)cycle
+       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+      enddo
+    enddo
+   endif
+  enddo
+ enddo
+
+END_PROVIDER 
+
+subroutine update_tc_rdm(h1,p1,h2,p2,s1,s2,array,sze,contrib)
+ implicit none
+ integer, intent(in) :: h1,p1,h2,p2,s1,s2,sze
+ double precision, intent(in) :: contrib
+ double precision, intent(inout) :: array(sze, sze, sze, sze)
+ integer :: istate
+ if(s1.ne.s2)then
+   array(p1,h1,p2,h2) += contrib
+   ! permutation for particle symmetry
+   array(p2,h2,p1,h1) += contrib
+ else ! same spin double excitation 
+   array(p1,h1,p2,h2) += contrib
+   ! exchange 
+   ! exchanging the holes 
+   array(p2,h1,p1,h2) -= contrib
+   ! exchanging the particles 
+   array(p1,h2,p2,h1) -= contrib
+
+   ! permutation for particle symmetry
+   array(p2,h2,p1,h1) += contrib
+   ! exchange 
+   ! exchanging the holes 
+   array(p1,h2,p2,h1) -= contrib
+   ! exchanging the particles 
+   array(p2,h1,p1,h2) -= contrib
+ endif
+
+end