diff --git a/src/tc_bi_ortho/test_tc_two_rdm.irp.f b/src/tc_bi_ortho/test_tc_two_rdm.irp.f
index ecdeef43..3e556312 100644
--- a/src/tc_bi_ortho/test_tc_two_rdm.irp.f
+++ b/src/tc_bi_ortho/test_tc_two_rdm.irp.f
@@ -35,6 +35,10 @@ subroutine  test
     do p2 = 1, mo_num
      integral = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
      rdm = tc_two_rdm(p1,h1,p2,h2)
+!     if(dabs(rdm).gt.1.d-10)then
+!      print*,h1,p1,h2,p2
+!      print*,rdm,integral,rdm*integral
+!     endif
      accu += integral * rdm
     enddo
    enddo
@@ -42,19 +46,28 @@ subroutine  test
  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)
+! print*,mo_bi_ortho_tc_two_e(2,15,2,1)
+! print*,mo_bi_ortho_tc_two_e(15,2,2,1)
+! print*,mo_bi_ortho_tc_two_e(2,1,2,15)
+! print*,mo_bi_ortho_tc_two_e(2,1,15,2)
  ref = 0.d0
  do i = 1, N_det
   do j = 1, N_det
-!   if(i.ne.j)cycle
+!   if(i.eq.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*,'i,j',i,j
+!  print*,psi_l_coef_bi_ortho(i,istate) , psi_r_coef_bi_ortho(j,istate) , htwoe 
+!  print*,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * htwoe 
+!   if(i.ne.j)then
 !  print*,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) ,  htwoe 
+!  print*,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * htwoe 
+!   endif
     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 
+ print*,'delta= ',ref-accu
 
 end
diff --git a/src/tc_bi_ortho/two_rdm_naive.irp.f b/src/tc_bi_ortho/two_rdm_naive.irp.f
index 9694c653..8fd34975 100644
--- a/src/tc_bi_ortho/two_rdm_naive.irp.f
+++ b/src/tc_bi_ortho/two_rdm_naive.irp.f
@@ -30,7 +30,6 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
     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) 
  
@@ -48,13 +47,12 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
       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
+!    cycle
     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
@@ -115,18 +113,12 @@ subroutine update_tc_rdm(h1,p1,h2,p2,s1,s2,array,sze,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(p2,h1,p1,h2) -= contrib
+   ! exchanging the 
    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