diff --git a/src/bi_ortho_mos/overlap.irp.f b/src/bi_ortho_mos/overlap.irp.f
index 09cec15f..1c88d16f 100644
--- a/src/bi_ortho_mos/overlap.irp.f
+++ b/src/bi_ortho_mos/overlap.irp.f
@@ -76,67 +76,85 @@ END_PROVIDER
 
 ! ---
 
-
  BEGIN_PROVIDER [ double precision, overlap_mo_r, (mo_num, mo_num)]
 &BEGIN_PROVIDER [ double precision, overlap_mo_l, (mo_num, mo_num)]
- implicit none
- BEGIN_DOC
-! overlap_mo_r_mo(j,i) = <MO_i|MO_R_j>
- END_DOC
- integer :: i,j,p,q
- overlap_mo_r= 0.d0
- overlap_mo_l= 0.d0
- do i = 1, mo_num
-   do j = 1, mo_num
-    do p = 1, ao_num
-     do q = 1, ao_num
-      overlap_mo_r(j,i) += mo_r_coef(q,i) * mo_r_coef(p,j) * ao_overlap(q,p) 
-      overlap_mo_l(j,i) += mo_l_coef(q,i) * mo_l_coef(p,j) * ao_overlap(q,p)
-     enddo
+
+  BEGIN_DOC
+  ! overlap_mo_r_mo(j,i) = <MO_i|MO_R_j>
+  END_DOC
+
+  implicit none
+  integer :: i, j, p, q
+
+  overlap_mo_r = 0.d0
+  overlap_mo_l = 0.d0
+  do i = 1, mo_num
+    do j = 1, mo_num
+      do p = 1, ao_num
+        do q = 1, ao_num
+          overlap_mo_r(j,i) += mo_r_coef(q,i) * mo_r_coef(p,j) * ao_overlap(q,p) 
+          overlap_mo_l(j,i) += mo_l_coef(q,i) * mo_l_coef(p,j) * ao_overlap(q,p)
+        enddo
+      enddo
     enddo
-   enddo
- enddo
+  enddo
+
 END_PROVIDER 
 
+! ---
+
  BEGIN_PROVIDER [ double precision, overlap_mo_r_mo, (mo_num, mo_num)]
 &BEGIN_PROVIDER [ double precision, overlap_mo_l_mo, (mo_num, mo_num)]
- implicit none
- BEGIN_DOC
-! overlap_mo_r_mo(j,i) = <MO_j|MO_R_i>
- END_DOC
- integer :: i,j,p,q
- overlap_mo_r_mo = 0.d0
- overlap_mo_l_mo = 0.d0
- do i = 1, mo_num
-   do j = 1, mo_num
-    do p = 1, ao_num
-     do q = 1, ao_num
-      overlap_mo_r_mo(j,i) += mo_coef(p,j) * mo_r_coef(q,i) * ao_overlap(q,p)
-      overlap_mo_l_mo(j,i) += mo_coef(p,j) * mo_l_coef(q,i) * ao_overlap(q,p)
-     enddo
+
+  BEGIN_DOC
+  ! overlap_mo_r_mo(j,i) = <MO_j|MO_R_i>
+  END_DOC
+
+  implicit none
+  integer :: i, j, p, q
+
+  overlap_mo_r_mo = 0.d0
+  overlap_mo_l_mo = 0.d0
+  do i = 1, mo_num
+    do j = 1, mo_num
+      do p = 1, ao_num
+        do q = 1, ao_num
+          overlap_mo_r_mo(j,i) += mo_coef(p,j) * mo_r_coef(q,i) * ao_overlap(q,p)
+          overlap_mo_l_mo(j,i) += mo_coef(p,j) * mo_l_coef(q,i) * ao_overlap(q,p)
+        enddo
+      enddo
     enddo
-   enddo
- enddo
+  enddo
+
 END_PROVIDER 
 
+! ---
+
  BEGIN_PROVIDER [ double precision, angle_left_right, (mo_num)]
 &BEGIN_PROVIDER [ double precision, max_angle_left_right]
- implicit none
- BEGIN_DOC
+
+  BEGIN_DOC
   ! angle_left_right(i) = angle between the left-eigenvector chi_i and the right-eigenvector phi_i
- END_DOC
- integer :: i,j
- double precision :: left,right,arg
- do i = 1, mo_num
-  left = overlap_mo_l(i,i)
-  right = overlap_mo_r(i,i)
-  arg = min(overlap_bi_ortho(i,i)/(left*right),1.d0)
-  arg = max(arg,-1.d0)
-  angle_left_right(i) = dacos(arg) * 180.d0/dacos(-1.d0)
- enddo
- double precision :: angle(mo_num)
- angle(1:mo_num) = dabs(angle_left_right(1:mo_num))
- max_angle_left_right = maxval(angle)
+  END_DOC
+
+  implicit none
+  integer          :: i, j
+  double precision :: left, right, arg
+  double precision :: angle(mo_num)
+
+  do i = 1, mo_num
+    left  = overlap_mo_l(i,i)
+    right = overlap_mo_r(i,i)
+    arg = min(overlap_bi_ortho(i,i)/(left*right),1.d0)
+    arg = max(arg, -1.d0)
+    angle_left_right(i) = dacos(arg) * 180.d0/dacos(-1.d0)
+  enddo
+
+  angle(1:mo_num) = dabs(angle_left_right(1:mo_num))
+  max_angle_left_right = maxval(angle)
+
 END_PROVIDER 
 
+! ---
+
 
diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f
index eda9642c..00c9dc4d 100644
--- a/src/cipsi/selection.irp.f
+++ b/src/cipsi/selection.irp.f
@@ -732,7 +732,10 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
         if (delta_E < 0.d0) then
             tmp = -tmp
         endif
+
+        !e_pert(istate) = alpha_h_psi * alpha_h_psi / (E0(istate) - Hii)
         e_pert(istate) = 0.5d0 * (tmp - delta_E)
+
         if (dabs(alpha_h_psi) > 1.d-4) then
           coef(istate) = e_pert(istate) / alpha_h_psi
         else
diff --git a/src/non_hermit_dav/biorthog.irp.f b/src/non_hermit_dav/biorthog.irp.f
index 633eff0d..0b8c78c2 100644
--- a/src/non_hermit_dav/biorthog.irp.f
+++ b/src/non_hermit_dav/biorthog.irp.f
@@ -252,7 +252,7 @@ end subroutine non_hrmt_real_diag_new
 
 ! ---
 
-subroutine non_hrmt_bieig(n, A, leigvec, reigvec, n_real_eigv, eigval)
+subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, eigval)
 
   BEGIN_DOC
   ! 
@@ -266,13 +266,14 @@ subroutine non_hrmt_bieig(n, A, leigvec, reigvec, n_real_eigv, eigval)
   implicit none
   integer,          intent(in)  :: n
   double precision, intent(in)  :: A(n,n)
+  double precision, intent(in)  :: thr_d, thr_nd
   integer,          intent(out) :: n_real_eigv
   double precision, intent(out) :: reigvec(n,n), leigvec(n,n), eigval(n)
 
   integer                       :: i, j
   integer                       :: n_good
   double precision              :: thr, thr_cut, thr_diag, thr_norm
-  double precision              :: accu_d, accu_nd, thr_d, thr_nd
+  double precision              :: accu_d, accu_nd
 
   integer,          allocatable :: list_good(:), iorder(:)
   double precision, allocatable :: WR(:), WI(:), VL(:,:), VR(:,:)
@@ -329,7 +330,7 @@ subroutine non_hrmt_bieig(n, A, leigvec, reigvec, n_real_eigv, eigval)
   n_good = 0
   thr    = 1.d-5
   do i = 1, n
-      print*, 'Re(i) + Im(i)', WR(i), WI(i)
+    print*, 'Re(i) + Im(i)', WR(i), WI(i)
     if(dabs(WI(i)) .lt. thr) then
       n_good += 1
     else
@@ -393,22 +394,22 @@ subroutine non_hrmt_bieig(n, A, leigvec, reigvec, n_real_eigv, eigval)
   ! -------------------------------------------------------------------------------------
   !                               check bi-orthogonality
 
-  thr_d  = 1d-10 ! -7
-  thr_nd = 1d-10 ! -7
+  thr_diag = 10.d0
+  thr_norm = 1d+10
 
   allocate( S(n_real_eigv,n_real_eigv) )
-  call check_biorthog(n, n_real_eigv, leigvec, reigvec, accu_d, accu_nd, S, .false.)
+  call check_biorthog(n, n_real_eigv, leigvec, reigvec, accu_d, accu_nd, S, thr_d, thr_nd, .false.)
 
-  if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv)) .lt. thr_d) ) then
+  if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv))/dble(n_real_eigv) .lt. thr_d) ) then
 
     print *, ' lapack vectors are normalized and bi-orthogonalized'
     deallocate(S)
     return
 
-  elseif( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv)) .gt. thr_d) ) then
+  elseif( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv))/dble(n_real_eigv) .gt. thr_d) ) then
 
     print *, ' lapack vectors are not normalized but bi-orthogonalized'
-    call check_biorthog_binormalize(n, n_real_eigv, leigvec, reigvec, .true.)
+    call check_biorthog_binormalize(n, n_real_eigv, leigvec, reigvec, thr_d, thr_nd, .true.)
 
     call check_EIGVEC(n, n, A, eigval, leigvec, reigvec, thr_diag, thr_norm, .true.)
 
@@ -427,17 +428,17 @@ subroutine non_hrmt_bieig(n, A, leigvec, reigvec, n_real_eigv, eigval)
     call impose_biorthog_degen_eigvec(n, eigval, leigvec, reigvec)
 
 
-    !call impose_orthog_biorthog_degen_eigvec(n, eigval, leigvec, reigvec)
+    !call impose_orthog_biorthog_degen_eigvec(n, thr_d, thr_nd, eigval, leigvec, reigvec)
 
     !call impose_unique_biorthog_degen_eigvec(n, eigval, mo_coef, ao_overlap, leigvec, reigvec)
 
     ! ---
 
-    call check_biorthog(n, n_real_eigv, leigvec, reigvec, accu_d, accu_nd, S, .false.)
+    call check_biorthog(n, n_real_eigv, leigvec, reigvec, accu_d, accu_nd, S, thr_d, thr_nd, .false.)
     if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv)) .gt. thr_d) ) then
-      call check_biorthog_binormalize(n, n_real_eigv, leigvec, reigvec, .true.)
+      call check_biorthog_binormalize(n, n_real_eigv, leigvec, reigvec, thr_d, thr_nd, .true.)
     endif
-    call check_biorthog(n, n_real_eigv, leigvec, reigvec, accu_d, accu_nd, S, .true.)
+    call check_biorthog(n, n_real_eigv, leigvec, reigvec, accu_d, accu_nd, S, thr_d, thr_nd, .true.)
 
     !call impose_biorthog_qr(n, n_real_eigv, leigvec, reigvec)
     !call impose_biorthog_lu(n, n_real_eigv, leigvec, reigvec)
diff --git a/src/non_hermit_dav/lapack_diag_non_hermit.irp.f b/src/non_hermit_dav/lapack_diag_non_hermit.irp.f
index d79e25ba..8d47f124 100644
--- a/src/non_hermit_dav/lapack_diag_non_hermit.irp.f
+++ b/src/non_hermit_dav/lapack_diag_non_hermit.irp.f
@@ -573,21 +573,22 @@ end subroutine non_hrmt_general_real_diag
 
 ! ---
 
-subroutine impose_biorthog_qr(m, n, Vl, Vr)
+subroutine impose_biorthog_qr(m, n, thr_d, thr_nd, Vl, Vr)
 
   implicit none 
-  integer, intent(in)             :: m, n
+  integer,          intent(in)    :: m, n
+  double precision, intent(in)    :: thr_d, thr_nd
   double precision, intent(inout) :: Vl(m,n), Vr(m,n)
 
   integer                         :: i, j
   integer                         :: LWORK, INFO
-  double precision                :: accu_nd, accu_d, thr_nd, thr_d
+  double precision                :: accu_nd, accu_d
   double precision, allocatable   :: TAU(:), WORK(:)
   double precision, allocatable   :: S(:,:), R(:,:), tmp(:,:)
 
   ! ---
 
-  call check_biorthog_binormalize(m, n, Vl, Vr, .false.)
+  call check_biorthog_binormalize(m, n, Vl, Vr, thr_d, thr_nd, .false.)
   
   ! ---
   
@@ -609,9 +610,7 @@ subroutine impose_biorthog_qr(m, n, Vl, Vr)
   enddo
   accu_nd = dsqrt(accu_nd)
 
-  thr_d  = 1d-10
-  thr_nd = 1d-08
-  if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n)) .lt. thr_d)) then
+  if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n))/dble(n) .lt. thr_d)) then
     print *, ' bi-orthogonal vectors without QR !'
     deallocate(S)
     return
@@ -1011,7 +1010,7 @@ subroutine check_degen(n, m, eigval, leigvec, reigvec)
   double precision                :: ei, ej, de, de_thr, accu_nd
   double precision, allocatable   :: S(:,:)
 
-  de_thr = 1d-7
+  de_thr = 1d-6
 
   do i = 1, m-1
     ei = eigval(i)
@@ -1082,7 +1081,7 @@ subroutine impose_weighted_orthog_svd(n, m, W, C)
   double precision, allocatable   :: S(:,:), tmp(:,:)
   double precision, allocatable   :: U(:,:), Vt(:,:), D(:)
 
-  print *, ' apply SVD to orthogonalize vectors'
+  print *, ' apply SVD to orthogonalize & normalize weighted vectors'
 
   ! ---
 
@@ -1097,7 +1096,7 @@ subroutine impose_weighted_orthog_svd(n, m, W, C)
             , 0.d0, S, size(S, 1) )
   deallocate(tmp)
 
-  print *, ' eigenvec overlap bef SVD: '
+  print *, ' overlap bef SVD: '
   do i = 1, m
     write(*, '(1000(F16.10,X))') S(i,:)
   enddo
@@ -1160,7 +1159,7 @@ subroutine impose_weighted_orthog_svd(n, m, W, C)
             , 0.d0, S, size(S, 1) )
   deallocate(tmp)
 
-  print *, ' eigenvec overlap aft SVD: '
+  print *, ' overlap aft SVD: '
   do i = 1, m
     write(*, '(1000(F16.10,X))') S(i,:)
   enddo
@@ -1185,7 +1184,7 @@ subroutine impose_orthog_svd(n, m, C)
   double precision, allocatable   :: S(:,:), tmp(:,:)
   double precision, allocatable   :: U(:,:), Vt(:,:), D(:)
 
-  print *, ' apply SVD to orthogonalize vectors'
+  print *, ' apply SVD to orthogonalize & normalize vectors'
 
   ! ---
 
@@ -1379,7 +1378,7 @@ subroutine impose_orthog_GramSchmidt(n, m, C)
   double precision, allocatable   :: S(:,:)
 
   print *, ''
-  print *, ' apply Gram-Schmidt to orthogonalize vectors'
+  print *, ' apply Gram-Schmidt to orthogonalize & normalize vectors'
   print *, ''
 
   ! ---
@@ -1663,22 +1662,19 @@ end subroutine get_halfinv_svd
 
 ! ---
 
-subroutine check_biorthog_binormalize(n, m, Vl, Vr, stop_ifnot)
+subroutine check_biorthog_binormalize(n, m, Vl, Vr, thr_d, thr_nd, stop_ifnot)
 
   implicit none
   
   integer,          intent(in)    :: n, m
   logical,          intent(in)    :: stop_ifnot
+  double precision, intent(in)    :: thr_d, thr_nd
   double precision, intent(inout) :: Vl(n,m), Vr(n,m)
 
   integer                         :: i, j
-  double precision                :: thr_d, thr_nd
   double precision                :: accu_d, accu_nd, s_tmp
   double precision, allocatable   :: S(:,:)
 
-  thr_d  = 1d-6
-  thr_nd = 1d-7
-
   print *, ' check bi-orthonormality'
 
   ! ---
@@ -1713,7 +1709,7 @@ subroutine check_biorthog_binormalize(n, m, Vl, Vr, stop_ifnot)
       endif
     enddo
   enddo
-  accu_nd = dsqrt(accu_nd)
+  accu_nd = dsqrt(accu_nd) / dble(m)
   print*, '    diag acc: ', accu_d
   print*, ' nondiag acc: ', accu_nd
 
@@ -1755,7 +1751,7 @@ subroutine check_biorthog_binormalize(n, m, Vl, Vr, stop_ifnot)
       endif
     enddo
   enddo
-  accu_nd = dsqrt(accu_nd)
+  accu_nd = dsqrt(accu_nd) / dble(m)
   print *, '    diag acc: ', accu_d
   print *, ' nondiag acc: ', accu_nd
 
@@ -1774,22 +1770,19 @@ end subroutine check_biorthog_binormalize
 
 ! ---
 
-subroutine check_weighted_biorthog(n, m, W, Vl, Vr, accu_d, accu_nd, S, stop_ifnot)
+subroutine check_weighted_biorthog(n, m, W, Vl, Vr, thr_d, thr_nd, accu_d, accu_nd, S, stop_ifnot)
 
   implicit none
   
   integer,          intent(in)  :: n, m
   double precision, intent(in)  :: Vl(n,m), Vr(n,m), W(n,n)
+  double precision, intent(in)  :: thr_d, thr_nd
   logical,          intent(in)  :: stop_ifnot
   double precision, intent(out) :: accu_d, accu_nd, S(m,m)
 
   integer                       :: i, j
-  double precision              :: thr_d, thr_nd
   double precision, allocatable :: SS(:,:), tmp(:,:)
 
-  thr_d  = 1d-6
-  thr_nd = 1d-08
-
   print *, ' check weighted bi-orthogonality'
 
   ! ---
@@ -1841,22 +1834,19 @@ end subroutine check_weighted_biorthog
 
 ! ---
 
-subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, stop_ifnot)
+subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_ifnot)
 
   implicit none
   
   integer,          intent(in)  :: n, m
   double precision, intent(in)  :: Vl(n,m), Vr(n,m)
   logical,          intent(in)  :: stop_ifnot
+  double precision, intent(in)  :: thr_d, thr_nd
   double precision, intent(out) :: accu_d, accu_nd, S(m,m)
 
   integer                       :: i, j
-  double precision              :: thr_d, thr_nd
   double precision, allocatable :: SS(:,:)
 
-  thr_d  = 1d-6
-  thr_nd = 1d-08
-
   print *, ' check bi-orthogonality'
 
   ! ---
@@ -1880,7 +1870,7 @@ subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, stop_ifnot)
       endif
     enddo
   enddo
-  accu_nd = dsqrt(accu_nd)
+  accu_nd = dsqrt(accu_nd) / dble(m)
 
   print *, ' accu_nd = ', accu_nd
   print *, ' accu_d  = ', dabs(accu_d-dble(m))/dble(m)
@@ -2029,7 +2019,7 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
 
       call impose_biorthog_svd(n, m, L, R)
 
-      !call impose_biorthog_qr(n, m, L, R)
+      !call impose_biorthog_qr(n, m, thr_d, thr_nd, L, R)
 
       ! ---
 
@@ -2047,11 +2037,12 @@ end subroutine impose_biorthog_degen_eigvec
 
 ! ---
 
-subroutine impose_orthog_biorthog_degen_eigvec(n, e0, L0, R0)
+subroutine impose_orthog_biorthog_degen_eigvec(n, thr_d, thr_nd, e0, L0, R0)
 
   implicit none
 
   integer,          intent(in)    :: n
+  double precision, intent(in)    :: thr_d, thr_nd
   double precision, intent(in)    :: e0(n)
   double precision, intent(inout) :: L0(n,n), R0(n,n)
 
@@ -2116,12 +2107,12 @@ subroutine impose_orthog_biorthog_degen_eigvec(n, e0, L0, R0)
 
       ! ---
   
-      call impose_biorthog_qr(n, m, L, R)
+      call impose_biorthog_qr(n, m, thr_d, thr_nd, L, R)
 
       allocate(S(m,m))
-      call check_biorthog(n, m, L, R, accu_d, accu_nd, S, .true.)
-      !call check_biorthog(n, m, L, L, accu_d, accu_nd, S, .true.)
-      !call check_biorthog(n, m, R, R, accu_d, accu_nd, S, .false.)
+      call check_biorthog(n, m, L, R, accu_d, accu_nd, S, thr_d, thr_nd, .true.)
+      !call check_biorthog(n, m, L, L, accu_d, accu_nd, S, thr_d, thr_nd, .true.)
+      !call check_biorthog(n, m, R, R, accu_d, accu_nd, S, thr_d, thr_nd, .false.)
       deallocate(S)
 
       ! ---
@@ -2140,11 +2131,12 @@ end subroutine impose_orthog_biorthog_degen_eigvec
 
 ! ---
 
-subroutine impose_unique_biorthog_degen_eigvec(n, e0, C0, W0, L0, R0)
+subroutine impose_unique_biorthog_degen_eigvec(n, thr_d, thr_nd, e0, C0, W0, L0, R0)
 
   implicit none
 
   integer,          intent(in)    :: n
+  double precision, intent(in)    :: thr_d, thr_nd
   double precision, intent(in)    :: e0(n), W0(n,n), C0(n,n)
   double precision, intent(inout) :: L0(n,n), R0(n,n)
 
@@ -2255,7 +2247,7 @@ subroutine impose_unique_biorthog_degen_eigvec(n, e0, C0, W0, L0, R0)
       call get_inv_half_nonsymmat_diago(S, m, S_inv_half, complex_root)
       if(complex_root)then
         call impose_biorthog_svd(n, m, L, R)
-        !call impose_biorthog_qr(n, m, L, R)
+        !call impose_biorthog_qr(n, m, thr_d, thr_nd, L, R)
       else
         call bi_ortho_s_inv_half(m, L, R, S_inv_half)
       endif
@@ -2502,8 +2494,288 @@ end subroutine impose_biorthog_svd
 
 ! ---
 
+subroutine impose_weighted_biorthog_qr(m, n, thr_d, thr_nd, Vl, W, Vr)
 
-subroutine impose_biorthog_svd_overlap(n, m, overlap, L, R)
+  implicit none 
+  integer,          intent(in)    :: m, n
+  double precision, intent(in)    :: thr_d, thr_nd
+  double precision, intent(inout) :: Vl(m,n), W(m,m), Vr(m,n)
+
+  integer                         :: i, j
+  integer                         :: LWORK, INFO
+  double precision                :: accu_nd, accu_d
+  double precision, allocatable   :: TAU(:), WORK(:)
+  double precision, allocatable   :: S(:,:), R(:,:), tmp(:,:), Stmp(:,:)
+
+
+  call check_weighted_biorthog_binormalize(m, n, Vl, W, Vr, thr_d, thr_nd, .false.)
+  
+  ! ---
+  
+  allocate(Stmp(n,m), S(n,n))
+  call dgemm( 'T', 'N', n, m, m, 1.d0        &
+            , Vl, size(Vl, 1), W, size(W, 1) &
+            , 0.d0, Stmp, size(Stmp, 1) )
+  call dgemm( 'N', 'N', n, n, m, 1.d0              &
+            , Stmp, size(Stmp, 1), Vr, size(Vr, 1) &
+            , 0.d0, S, size(S, 1) )
+  deallocate(Stmp)
+
+  accu_nd = 0.d0
+  accu_d  = 0.d0
+  do i = 1, n
+    do j = 1, n
+      if(i==j) then
+        accu_d += S(j,i)
+      else
+        accu_nd = accu_nd + S(j,i) * S(j,i)
+      endif
+    enddo
+  enddo
+  accu_nd = dsqrt(accu_nd)
+
+  if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n))/dble(n) .lt. thr_d)) then
+    print *, ' bi-orthogonal vectors without QR !'
+    deallocate(S)
+    return
+  endif
+
+  ! -------------------------------------------------------------------------------------
+  !                           QR factorization of S: S = Q x R
+
+
+  print *, ' apply QR decomposition ...'
+
+  allocate( TAU(n), WORK(1) )
+
+  LWORK = -1
+  call dgeqrf(n, n, S, n, TAU, WORK, LWORK, INFO)
+  if(INFO .ne. 0) then
+    print*,'dgeqrf failed !!', INFO
+    stop
+  endif
+
+  LWORK = max(n, int(WORK(1)))
+  deallocate(WORK)
+
+  allocate( WORK(LWORK) )
+  call dgeqrf(n, n, S, n, TAU, WORK, LWORK, INFO)
+  if(INFO .ne. 0) then
+    print*,'dgeqrf failed !!', INFO
+    stop
+  endif
+
+  ! save the upper triangular R
+  allocate( R(n,n) )
+  R(:,:) = S(:,:)
+
+  ! get Q
+  LWORK = -1
+  call dorgqr(n, n, n, S, n, TAU, WORK, LWORK, INFO)
+  if(INFO .ne. 0) then
+    print*,'dorgqr failed !!', INFO
+    stop
+  endif
+
+  LWORK = max(n, int(WORK(1)))
+  deallocate(WORK)
+
+  allocate( WORK(LWORK) )
+  call dorgqr(n, n, n, S, n, TAU, WORK, LWORK, INFO)
+  if(INFO .ne. 0) then
+    print*,'dorgqr failed !!', INFO
+    stop
+  endif
+
+  deallocate( WORK, TAU )
+
+  !
+  ! -------------------------------------------------------------------------------------
+
+  ! ---
+
+  ! -------------------------------------------------------------------------------------
+  !                               get bi-orhtog left & right vectors:
+  !                                           Vr' = Vr x inv(R) 
+  !                                           Vl' = inv(Q) x Vl =  Q.T   x Vl 
+
+  ! Q.T x Vl, where Q = S
+
+  allocate( tmp(n,m) )
+  call dgemm( 'T', 'T', n, m, n, 1.d0        &
+            , S, size(S, 1), Vl, size(Vl, 1) &
+            , 0.d0, tmp, size(tmp, 1) )
+
+  do i = 1, n
+    do j = 1, m
+      Vl(j,i) = tmp(i,j)
+    enddo
+  enddo
+  deallocate(tmp)
+
+  ! ---
+
+  ! inv(R) 
+  !print *, ' inversing upper triangular matrix ...'
+  call dtrtri("U", "N", n, R, n, INFO)
+  if(INFO .ne. 0) then
+    print*,'dtrtri failed !!', INFO
+    stop
+  endif
+  !print *, ' inversing upper triangular matrix OK' 
+
+  do i = 1, n-1
+    do j = i+1, n
+      R(j,i) = 0.d0
+    enddo
+  enddo
+
+  !print *, ' inv(R):'
+  !do i = 1, n
+  !  write(*, '(1000(F16.10,X))') R(i,:)
+  !enddo
+
+  ! Vr x inv(R) 
+  allocate( tmp(m,n) )
+  call dgemm( 'N', 'N', m, n, n, 1.d0        &
+            , Vr, size(Vr, 1), R, size(R, 1) &
+            , 0.d0, tmp, size(tmp, 1) )
+  deallocate( R )
+
+  do i = 1, n
+    do j = 1, m
+      Vr(j,i) = tmp(j,i)
+    enddo
+  enddo
+  deallocate(tmp)
+
+  call check_weighted_biorthog_binormalize(m, n, Vl, W, Vr, thr_d, thr_nd, .false.)
+
+  return
+end subroutine impose_weighted_biorthog_qr
+
+! ---
+
+subroutine check_weighted_biorthog_binormalize(n, m, Vl, W, Vr, thr_d, thr_nd, stop_ifnot)
+
+  implicit none
+  
+  integer,          intent(in)    :: n, m
+  logical,          intent(in)    :: stop_ifnot
+  double precision, intent(in)    :: thr_d, thr_nd
+  double precision, intent(inout) :: Vl(n,m), W(n,n), Vr(n,m)
+
+  integer                         :: i, j
+  double precision                :: accu_d, accu_nd, s_tmp
+  double precision, allocatable   :: S(:,:), Stmp(:,:)
+
+  print *, ' check weighted bi-orthonormality'
+
+  ! ---
+
+  allocate(Stmp(m,n), S(m,m))
+  call dgemm( 'T', 'N', m, n, n, 1.d0        &
+            , Vl, size(Vl, 1), W, size(W, 1) &
+            , 0.d0, Stmp, size(Stmp, 1) )
+  call dgemm( 'N', 'N', m, m, n, 1.d0              &
+            , Stmp, size(Stmp, 1), Vr, size(Vr, 1) &
+            , 0.d0, S, size(S, 1) )
+  deallocate(Stmp)
+  !print *, ' overlap matrix before:'
+  !do i = 1, m
+  !  write(*,'(1000(F16.10,X))') S(i,:)
+  !enddo
+
+  ! S(i,i) = -1
+  do i = 1, m
+    if( (S(i,i) + 1.d0) .lt. thr_d ) then
+      do j = 1, n
+        Vl(j,i) = -1.d0 * Vl(j,i)
+      enddo
+      S(i,i) = 1.d0
+    endif
+  enddo
+
+  accu_d  = 0.d0
+  accu_nd = 0.d0
+  do i = 1, m
+    do j = 1, m
+      if(i==j) then
+        accu_d = accu_d + S(i,i)
+      else
+        accu_nd = accu_nd + S(j,i) * S(j,i)
+      endif
+    enddo
+  enddo
+  accu_nd = dsqrt(accu_nd) / dble(m)
+  print*, '    diag acc: ', accu_d
+  print*, ' nondiag acc: ', accu_nd
+
+  ! ---
+
+  if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(m))/dble(m) .gt. thr_d) ) then
+
+    do i = 1, m
+      print *, i, S(i,i)
+      if(dabs(S(i,i) - 1.d0) .gt. thr_d) then
+        s_tmp = 1.d0 / dsqrt(S(i,i))
+        do j = 1, n
+          Vl(j,i) = Vl(j,i) * s_tmp 
+          Vr(j,i) = Vr(j,i) * s_tmp 
+        enddo
+      endif
+    enddo
+
+  endif
+
+  ! ---
+
+  allocate(Stmp(m,n))
+  call dgemm( 'T', 'N', m, n, n, 1.d0        &
+            , Vl, size(Vl, 1), W, size(W, 1) &
+            , 0.d0, Stmp, size(Stmp, 1) )
+  call dgemm( 'N', 'N', m, m, n, 1.d0              &
+            , Stmp, size(Stmp, 1), Vr, size(Vr, 1) &
+            , 0.d0, S, size(S, 1) )
+  deallocate(Stmp)
+  !print *, ' overlap matrix after:'
+  !do i = 1, m
+  !  write(*,'(1000(F16.10,X))') S(i,:)
+  !enddo
+
+  accu_d  = 0.d0
+  accu_nd = 0.d0
+  do i = 1, m
+    do j = 1, m
+      if(i==j) then
+        accu_d = accu_d + S(i,i)
+      else
+        accu_nd = accu_nd + S(j,i) * S(j,i)
+      endif
+    enddo
+  enddo
+  accu_nd = dsqrt(accu_nd) / dble(m)
+  print *, '    diag acc: ', accu_d
+  print *, ' nondiag acc: ', accu_nd
+
+  deallocate(S)
+
+  ! ---
+
+  if( stop_ifnot .and. ((accu_nd .gt. thr_nd) .or. (dabs(accu_d-dble(m))/dble(m) .gt. thr_d)) ) then
+    print *, accu_nd, thr_nd 
+    print *, dabs(accu_d-dble(m))/dble(m), thr_d
+    print *, ' weighted biorthog_binormalize failed !'
+    stop
+  endif
+
+end subroutine check_weighted_biorthog_binormalize
+
+! ---
+
+
+
+subroutine impose_weighted_biorthog_svd(n, m, overlap, L, R)
 
   implicit none
 
@@ -2527,6 +2799,7 @@ subroutine impose_biorthog_svd_overlap(n, m, overlap, L, R)
   call dgemm( 'T', 'N', m, m, n, 1.d0      &
             , L, size(L, 1), Stmp, size(Stmp, 1) &
             , 0.d0, S, size(S, 1) )
+  deallocate(Stmp)
 
   print *, ' overlap bef SVD: '
   do i = 1, m
@@ -2598,10 +2871,7 @@ subroutine impose_biorthog_svd_overlap(n, m, overlap, L, R)
 
   ! ---
 
-  allocate(S(m,m))
-!  call dgemm( 'T', 'N', m, m, n, 1.d0      &
-!            , L, size(L, 1), R, size(R, 1) &
-!            , 0.d0, S, size(S, 1) )
+  allocate(S(m,m),Stmp(n,m))
   ! S = C.T x overlap x C
   call dgemm( 'N', 'N', n, m, n, 1.d0      &
             , overlap, size(overlap, 1), R, size(R, 1) &
@@ -2609,6 +2879,7 @@ subroutine impose_biorthog_svd_overlap(n, m, overlap, L, R)
   call dgemm( 'T', 'N', m, m, n, 1.d0      &
             , L, size(L, 1), Stmp, size(Stmp, 1) &
             , 0.d0, S, size(S, 1) )
+  deallocate(Stmp)
 
   print *, ' overlap aft SVD with overlap: '
   do i = 1, m
@@ -2618,7 +2889,7 @@ subroutine impose_biorthog_svd_overlap(n, m, overlap, L, R)
 
   ! ---
 
-end subroutine impose_biorthog_svd
+end subroutine impose_weighted_biorthog_svd
 
 ! ---
 
diff --git a/src/tc_bi_ortho/tc_natorb.irp.f b/src/tc_bi_ortho/tc_natorb.irp.f
index 4144fcad..b2ccb091 100644
--- a/src/tc_bi_ortho/tc_natorb.irp.f
+++ b/src/tc_bi_ortho/tc_natorb.irp.f
@@ -1,125 +1,150 @@
 
-  BEGIN_PROVIDER [ double precision, natorb_tc_reigvec_mo, (mo_num, mo_num)]
- &BEGIN_PROVIDER [ double precision, natorb_tc_leigvec_mo, (mo_num, mo_num)]
- &BEGIN_PROVIDER [ double precision, natorb_tc_eigval, (mo_num)]
+! ---
+
+ BEGIN_PROVIDER [ double precision, natorb_tc_reigvec_mo, (mo_num, mo_num)]
+&BEGIN_PROVIDER [ double precision, natorb_tc_leigvec_mo, (mo_num, mo_num)]
+&BEGIN_PROVIDER [ double precision, natorb_tc_eigval, (mo_num)]
+
+  BEGIN_DOC
+  !
+  ! natorb_tc_reigvec_mo : RIGHT eigenvectors of the ground state transition matrix (equivalent of natural orbitals)
+  ! natorb_tc_leigvec_mo : LEFT  eigenvectors of the ground state transition matrix (equivalent of natural orbitals)
+  ! natorb_tc_eigval     : eigenvalues of the ground state transition matrix (equivalent of the occupation numbers). WARNINING :: can be negative !!
+  !
+  END_DOC
+
   implicit none
- BEGIN_DOC
- ! natorb_tc_reigvec_mo : RIGHT eigenvectors of the ground state transition matrix (equivalent of natural orbitals)
- ! natorb_tc_leigvec_mo : LEFT  eigenvectors of the ground state transition matrix (equivalent of natural orbitals)
- ! natorb_tc_eigval     : eigenvalues of the ground state transition matrix (equivalent of the occupation numbers). WARNINING :: can be negative !!
- END_DOC
-  double precision, allocatable :: dm_tmp(:,:),fock_diag(:)
-  double precision :: thr_deg
-  integer :: i,j,k,n_real
-  allocate( dm_tmp(mo_num,mo_num),fock_diag(mo_num))
+  integer                       :: i, j, k, n_real
+  double precision              :: thr_d, thr_nd, thr_deg, accu
+  double precision              :: accu_d, accu_nd
+  double precision, allocatable :: dm_tmp(:,:), fock_diag(:)
+
+  allocate(dm_tmp(mo_num,mo_num), fock_diag(mo_num))
+
   dm_tmp(:,:) = -tc_transition_matrix(:,:,1,1)
-  print*,'dm_tmp'
+
+  print *, ' dm_tmp'
   do i = 1, mo_num
-   fock_diag(i) = fock_matrix_tc_mo_tot(i,i)
-   write(*,'(100(F16.10,X))')-dm_tmp(:,i)
+    fock_diag(i) = fock_matrix_tc_mo_tot(i,i)
+    write(*, '(100(F16.10,X))') -dm_tmp(:,i)
   enddo
+
+  thr_d   = 1.d-6
+  thr_nd  = 1.d-6
   thr_deg = 1.d-3
-  call diag_mat_per_fock_degen(fock_diag,dm_tmp,mo_num,thr_deg,& 
-                               natorb_tc_leigvec_mo,natorb_tc_reigvec_mo,&
-                               natorb_tc_eigval)
+  call diag_mat_per_fock_degen( fock_diag, dm_tmp, mo_num, thr_d, thr_nd, thr_deg & 
+                              , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)
 !   call non_hrmt_bieig( mo_num, dm_tmp&
 !                      , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo& 
 !                      , n_real, natorb_tc_eigval )
- double precision :: accu
+
   accu = 0.d0
   do i = 1, n_real
-   print*,'natorb_tc_eigval(i) = ',-natorb_tc_eigval(i)
-   accu += -natorb_tc_eigval(i)
+    print*,'natorb_tc_eigval(i) = ',-natorb_tc_eigval(i)
+    accu += -natorb_tc_eigval(i)
   enddo
-  print*,'accu = ',accu
+  print *, ' accu = ', accu
+
   dm_tmp = 0.d0
   do i = 1, n_real
-   accu = 0.d0
-   do k = 1, mo_num
-    accu += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,i)
-   enddo
-   accu = 1.d0/dsqrt(dabs(accu))
-   natorb_tc_reigvec_mo(:,i) *= accu
-   natorb_tc_leigvec_mo(:,i) *= accu
-   do j = 1, n_real
+    accu = 0.d0
     do k = 1, mo_num
-     dm_tmp(j,i) += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,j)
+      accu += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,i)
+    enddo
+    accu = 1.d0/dsqrt(dabs(accu))
+    natorb_tc_reigvec_mo(:,i) *= accu
+    natorb_tc_leigvec_mo(:,i) *= accu
+    do j = 1, n_real
+      do k = 1, mo_num
+        dm_tmp(j,i) += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,j)
+      enddo
     enddo
-   enddo
   enddo
-  double precision :: accu_d, accu_nd
-  accu_d = 0.d0
+
+  accu_d  = 0.d0
   accu_nd = 0.d0
   do i = 1, mo_num
-   accu_d += dm_tmp(i,i)
- !  write(*,'(100(F16.10,X))')dm_tmp(:,i)
-   do j = 1, mo_num
-    if(i==j)cycle
-    accu_nd += dabs(dm_tmp(j,i))
-   enddo
+    accu_d += dm_tmp(i,i)
+    !write(*,'(100(F16.10,X))')dm_tmp(:,i)
+    do j = 1, mo_num
+      if(i==j)cycle
+      accu_nd += dabs(dm_tmp(j,i))
+    enddo
   enddo
-  print*,'Trace of the overlap between TC natural orbitals     ',accu_d
-  print*,'L1 norm of extra diagonal elements of overlap matrix ',accu_nd
+  print *, ' Trace of the overlap between TC natural orbitals     ', accu_d
+  print *, ' L1 norm of extra diagonal elements of overlap matrix ', accu_nd
+
+  deallocate(dm_tmp, fock_diag)
  
+END_PROVIDER 
+
+! ---
  
- END_PROVIDER 
- 
-  BEGIN_PROVIDER [ double precision, fock_diag_sorted_r_natorb, (mo_num, mo_num)]
- &BEGIN_PROVIDER [ double precision, fock_diag_sorted_l_natorb, (mo_num, mo_num)]
- &BEGIN_PROVIDER [ double precision, fock_diag_sorted_v_natorb, (mo_num)]
+ BEGIN_PROVIDER [ double precision, fock_diag_sorted_r_natorb, (mo_num, mo_num)]
+&BEGIN_PROVIDER [ double precision, fock_diag_sorted_l_natorb, (mo_num, mo_num)]
+&BEGIN_PROVIDER [ double precision, fock_diag_sorted_v_natorb, (mo_num)]
+
   implicit none
-  integer ::i,j,k
-  print*,'Diagonal elements of the Fock matrix before '
-  do i = 1, mo_num
-   write(*,*)i,Fock_matrix_tc_mo_tot(i,i)
-  enddo
+  integer                       :: i,j,k
+  integer,          allocatable :: iorder(:)
   double precision, allocatable :: fock_diag(:)
+
+  print *, ' Diagonal elements of the Fock matrix before '
+
+  do i = 1, mo_num
+   write(*,*) i, Fock_matrix_tc_mo_tot(i,i)
+  enddo
+
   allocate(fock_diag(mo_num))
   fock_diag = 0.d0
   do i = 1, mo_num
-   fock_diag(i) = 0.d0
-   do j = 1, mo_num
-    do k = 1, mo_num
-     fock_diag(i) += natorb_tc_leigvec_mo(k,i) * Fock_matrix_tc_mo_tot(k,j) * natorb_tc_reigvec_mo(j,i) 
+    fock_diag(i) = 0.d0
+    do j = 1, mo_num
+      do k = 1, mo_num
+        fock_diag(i) += natorb_tc_leigvec_mo(k,i) * Fock_matrix_tc_mo_tot(k,j) * natorb_tc_reigvec_mo(j,i) 
+      enddo
     enddo
-   enddo
   enddo
-  integer, allocatable :: iorder(:)
+
   allocate(iorder(mo_num))
   do i = 1, mo_num
    iorder(i) = i
   enddo 
-  call dsort(fock_diag,iorder,mo_num)
-  print*,'Diagonal elements of the Fock matrix after '
+  call dsort(fock_diag, iorder, mo_num)
+
+  print *, ' Diagonal elements of the Fock matrix after '
   do i = 1, mo_num
-   write(*,*)i,fock_diag(i)
+   write(*,*) i, fock_diag(i)
   enddo
+  deallocate(fock_diag)
+
   do i = 1, mo_num 
-   fock_diag_sorted_v_natorb(i) = natorb_tc_eigval(iorder(i))
-   do j = 1, mo_num
-    fock_diag_sorted_r_natorb(j,i) = natorb_tc_reigvec_mo(j,iorder(i))
-    fock_diag_sorted_l_natorb(j,i) = natorb_tc_leigvec_mo(j,iorder(i))
-   enddo
+    fock_diag_sorted_v_natorb(i) = natorb_tc_eigval(iorder(i))
+    do j = 1, mo_num
+      fock_diag_sorted_r_natorb(j,i) = natorb_tc_reigvec_mo(j,iorder(i))
+      fock_diag_sorted_l_natorb(j,i) = natorb_tc_leigvec_mo(j,iorder(i))
+    enddo
   enddo
+  deallocate(iorder)
  
- END_PROVIDER 
+END_PROVIDER 
  
+! --- 
  
+ BEGIN_PROVIDER [ double precision, natorb_tc_reigvec_ao, (ao_num, mo_num)]
+&BEGIN_PROVIDER [ double precision, natorb_tc_leigvec_ao, (ao_num, mo_num)]
+&BEGIN_PROVIDER [ double precision, overlap_natorb_tc_eigvec_ao, (mo_num, mo_num) ]
  
-  BEGIN_PROVIDER [ double precision, natorb_tc_reigvec_ao, (ao_num, mo_num)]
- &BEGIN_PROVIDER [ double precision, natorb_tc_leigvec_ao, (ao_num, mo_num)]
- &BEGIN_PROVIDER [ double precision, overlap_natorb_tc_eigvec_ao, (mo_num, mo_num) ]
+  BEGIN_DOC
+  ! EIGENVECTORS OF FOCK MATRIX ON THE AO BASIS and their OVERLAP
+  !
+  ! THE OVERLAP SHOULD BE THE SAME AS overlap_natorb_tc_eigvec_mo
+  END_DOC
  
-   BEGIN_DOC
-   ! EIGENVECTORS OF FOCK MATRIX ON THE AO BASIS and their OVERLAP
-   !
-   ! THE OVERLAP SHOULD BE THE SAME AS overlap_natorb_tc_eigvec_mo
-   END_DOC
- 
-   implicit none
-   integer                       :: i, j, k, q, p
-   double precision              :: accu, accu_d
-   double precision, allocatable :: tmp(:,:)
+  implicit none
+  integer                       :: i, j, k, q, p
+  double precision              :: accu, accu_d
+  double precision, allocatable :: tmp(:,:)
  
  
  !  ! MO_R x R
diff --git a/src/tc_scf/diago_bi_ort_tcfock.irp.f b/src/tc_scf/diago_bi_ort_tcfock.irp.f
index 00558de6..856b7382 100644
--- a/src/tc_scf/diago_bi_ort_tcfock.irp.f
+++ b/src/tc_scf/diago_bi_ort_tcfock.irp.f
@@ -11,20 +11,24 @@
   integer                       :: n_real_tc 
   integer                       :: i, k, l
   double precision              :: accu_d, accu_nd, accu_tmp
+  double precision              :: thr_d, thr_nd
   double precision              :: norm
   double precision, allocatable :: eigval_right_tmp(:)
 
+  thr_d  = 1d-6
+  thr_nd = 1d-6
+
   allocate( eigval_right_tmp(mo_num) )
 
   PROVIDE Fock_matrix_tc_mo_tot
 
- call non_hrmt_bieig( mo_num, Fock_matrix_tc_mo_tot          &
-                    , fock_tc_leigvec_mo, fock_tc_reigvec_mo & 
-                    , n_real_tc, eigval_right_tmp )
+  call non_hrmt_bieig( mo_num, Fock_matrix_tc_mo_tot, thr_d, thr_nd &
+                     , fock_tc_leigvec_mo, fock_tc_reigvec_mo       & 
+                     , n_real_tc, eigval_right_tmp )
   !if(max_ov_tc_scf)then
-!    call non_hrmt_fock_mat( mo_num, Fock_matrix_tc_mo_tot &
-!                       , fock_tc_leigvec_mo, fock_tc_reigvec_mo & 
-!                       , n_real_tc, eigval_right_tmp )
+  ! call non_hrmt_fock_mat( mo_num, Fock_matrix_tc_mo_tot, thr_d, thr_nd &
+  !                    , fock_tc_leigvec_mo, fock_tc_reigvec_mo & 
+  !                    , n_real_tc, eigval_right_tmp )
   !else 
   ! call non_hrmt_diag_split_degen_bi_orthog( mo_num, Fock_matrix_tc_mo_tot &
   !                    , fock_tc_leigvec_mo, fock_tc_reigvec_mo & 
@@ -59,17 +63,17 @@
       else
         accu_tmp = overlap_fock_tc_eigvec_mo(k,i)
         accu_nd += accu_tmp * accu_tmp
-        if(dabs(overlap_fock_tc_eigvec_mo(k,i)).gt.1.d-10)then
-         print*,'k,i',k,i,overlap_fock_tc_eigvec_mo(k,i)
+        if(dabs(overlap_fock_tc_eigvec_mo(k,i)) .gt. thr_nd)then
+         print *, 'k,i', k, i, overlap_fock_tc_eigvec_mo(k,i)
         endif
       endif
     enddo 
   enddo
   accu_nd = dsqrt(accu_nd)/accu_d
 
-  if( accu_nd .gt. 1d-8 ) then
+  if(accu_nd .gt. thr_nd) then
     print *, ' bi-orthog failed'
-    print*,'accu_nd MO = ', accu_nd
+    print*,'accu_nd MO = ', accu_nd, thr_nd
     print*,'overlap_fock_tc_eigvec_mo = '
     do i = 1, mo_num
       write(*,'(100(F16.10,X))') overlap_fock_tc_eigvec_mo(i,:)
@@ -77,13 +81,13 @@
    stop
   endif
 
-  if( dabs(accu_d - dble(mo_num)) .gt. 1e-7 ) then
+  if( dabs(accu_d - dble(mo_num))/dble(mo_num) .gt. thr_d ) then
     print *, 'mo_num     = ', mo_num 
-    print *, 'accu_d  MO = ', accu_d
+    print *, 'accu_d  MO = ', accu_d, thr_d
     print *, 'normalizing vectors ...'
     do i = 1, mo_num
       norm = dsqrt(dabs(overlap_fock_tc_eigvec_mo(i,i)))
-      if( norm.gt.1e-7 ) then
+      if(norm .gt. thr_d) then
         do k = 1, mo_num
           fock_tc_reigvec_mo(k,i) *= 1.d0/norm
           fock_tc_leigvec_mo(k,i) *= 1.d0/norm
diff --git a/src/tc_scf/fock_tc.irp.f b/src/tc_scf/fock_tc.irp.f
index 4907f2c8..0cb53295 100644
--- a/src/tc_scf/fock_tc.irp.f
+++ b/src/tc_scf/fock_tc.irp.f
@@ -116,17 +116,23 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ]
   endif
 END_PROVIDER
 
+! ---
 
 BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_tot, (mo_num, mo_num)]
-  implicit none
+
   BEGIN_DOC
- ! Total alpha+beta  TC Fock matrix : h_c + Two-e^TC terms on the MO basis
+  ! Total alpha+beta  TC Fock matrix : h_c + Two-e^TC terms on the MO basis
   END_DOC
+
+  implicit none
+
   Fock_matrix_tc_mo_tot = 0.5d0 * (Fock_matrix_tc_mo_alpha + Fock_matrix_tc_mo_beta)
   if(three_body_h_tc) then
     Fock_matrix_tc_mo_tot += fock_3_mat
   endif
+
   !call restore_symmetry(mo_num, mo_num, Fock_matrix_tc_mo_tot, mo_num, 1.d-10)
+
 END_PROVIDER 
 
 ! ---
diff --git a/src/tc_scf/rotate_tcscf_orbitals.irp.f b/src/tc_scf/rotate_tcscf_orbitals.irp.f
index d32d324d..d53991ed 100644
--- a/src/tc_scf/rotate_tcscf_orbitals.irp.f
+++ b/src/tc_scf/rotate_tcscf_orbitals.irp.f
@@ -73,7 +73,7 @@ subroutine maximize_overlap()
 
   ! ---
 
-  call rotate_degen_eigvec(n, m, e, C, W, L, R)
+  call rotate_degen_eigvec_to_maximize_overlap(n, m, e, C, W, L, R)
 
   ! ---
    
@@ -115,7 +115,7 @@ end subroutine maximize_overlap
 
 ! ---
 
-subroutine rotate_degen_eigvec(n, m, e0, C0, W0, L0, R0)
+subroutine rotate_degen_eigvec_to_maximize_overlap(n, m, e0, C0, W0, L0, R0)
 
   implicit none
 
@@ -124,7 +124,7 @@ subroutine rotate_degen_eigvec(n, m, e0, C0, W0, L0, R0)
   double precision, intent(inout) :: L0(n,m), R0(n,m)
 
 
-  integer                         :: i, j, k, kk, mm, id1
+  integer                         :: i, j, k, kk, mm, id1, tot_deg
   double precision                :: ei, ej, de, de_thr
   integer,          allocatable   :: deg_num(:)
   double precision, allocatable   :: L(:,:), R(:,:), C(:,:), Lnew(:,:), Rnew(:,:), tmp(:,:)
@@ -162,12 +162,19 @@ subroutine rotate_degen_eigvec(n, m, e0, C0, W0, L0, R0)
     enddo
   enddo
 
+  tot_deg = 0
   do i = 1, m
     if(deg_num(i).gt.1) then
       print *, ' degen on', i, deg_num(i)
+      tot_deg = tot_deg + 1
     endif
   enddo
 
+  if(tot_deg .eq. 0) then
+    print *, ' no degen'
+    return
+  endif
+
   ! ---
 
   do i = 1, m
@@ -243,6 +250,122 @@ subroutine rotate_degen_eigvec(n, m, e0, C0, W0, L0, R0)
 
   deallocate(S, Snew, T)
 
-end subroutine rotate_degen_eigvec
+end subroutine rotate_degen_eigvec_to_maximize_overlap
+
+! ---
+
+subroutine fix_right_to_one()
+
+  implicit none
+  integer                       :: i, j, m, n, mm, tot_deg
+  double precision              :: accu_d, accu_nd
+  double precision              :: de_thr, ei, ej, de
+  double precision              :: thr_d, thr_nd
+  integer,          allocatable :: deg_num(:)
+  double precision, allocatable :: R0(:,:), L0(:,:), W(:,:), e0(:)
+  double precision, allocatable :: R(:,:), L(:,:), S(:,:), Stmp(:,:), tmp(:,:)
+
+  thr_d  = 1d-7
+  thr_nd = 1d-7
+
+  n = ao_num
+  m = mo_num
+
+  allocate(L0(n,m), R0(n,m), W(n,n), e0(m))
+  L0 = mo_l_coef
+  R0 = mo_r_coef
+  W  = ao_overlap
+
+  print*, ' fock matrix diag elements'
+  do i = 1, m
+    e0(i) = Fock_matrix_tc_mo_tot(i,i)
+    print*, e0(i)
+  enddo
+
+  ! ---
+
+  allocate( deg_num(m) )
+  do i = 1, m
+    deg_num(i) = 1
+  enddo
+
+  de_thr = 1d-6
+
+  do i = 1, m-1
+    ei = e0(i)
+
+    ! already considered in degen vectors
+    if(deg_num(i).eq.0) cycle
+
+    do j = i+1, m
+      ej = e0(j)
+      de = dabs(ei - ej)
+
+      if(de .lt. de_thr) then
+        deg_num(i) = deg_num(i) + 1
+        deg_num(j) = 0
+      endif
+
+    enddo
+  enddo
+
+  deallocate(e0)
+
+  tot_deg = 0
+  do i = 1, m
+    if(deg_num(i).gt.1) then
+      print *, ' degen on', i, deg_num(i)
+      tot_deg = tot_deg + 1
+    endif
+  enddo
+
+  if(tot_deg .eq. 0) then
+    print *, ' no degen'
+    return
+  endif
+
+  ! ---
+
+  do i = 1, m
+    mm = deg_num(i)
+
+    if(mm .gt. 1) then
+
+      allocate(L(n,mm), R(n,mm))
+      do j = 1, mm
+        L(1:n,j) = L0(1:n,i+j-1)
+        R(1:n,j) = R0(1:n,i+j-1)
+      enddo
+
+      ! ---
+
+      call impose_weighted_orthog_svd(n, mm, W, R)
+      call impose_weighted_biorthog_qr(n, mm, thr_d, thr_nd, R, W, L)
+
+      ! ---
+
+      do j = 1, mm
+        L0(1:n,i+j-1) = L(1:n,j)
+        R0(1:n,i+j-1) = R(1:n,j)
+      enddo
+      deallocate(L, R)
+
+    endif
+  enddo
+
+  call check_weighted_biorthog_binormalize(n, m, L0, W, R0, thr_d, thr_nd, .true.)
+
+  deallocate(W, deg_num)
+
+  mo_l_coef = L0
+  mo_r_coef = R0
+  deallocate(L0, R0)
+
+  call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
+  call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
+  print *, ' orbitals are rotated '
+
+  return
+end subroutine fix_right_to_one
 
 ! ---
diff --git a/src/tc_scf/routines_rotates.irp.f b/src/tc_scf/routines_rotates.irp.f
index 51c218fa..1b03f338 100644
--- a/src/tc_scf/routines_rotates.irp.f
+++ b/src/tc_scf/routines_rotates.irp.f
@@ -1,203 +1,242 @@
-subroutine minimize_tc_orb_angles
- implicit none
- double precision :: thr_deg
- logical :: good_angles
- integer :: i
- good_angles = .False.
- thr_deg = thr_degen_tc
- call print_energy_and_mos
- i = 1
- print*,'Minimizing the angles between the TC orbitals'
- do while (.not. good_angles)
-  print*,'iteration = ',i
-  call routine_save_rotated_mos(thr_deg,good_angles)
-  thr_deg *= 10.d0
-  i+=1
-  if(i.gt.100)then
-   print*,'minimize_tc_orb_angles does not seem to converge ..'
-   print*,'Something is weird in the tc orbitals ...'
-   print*,'STOPPING'
-  endif
- enddo
- print*,'Converged ANGLES MINIMIZATION !!'
- call print_angles_tc
- call print_energy_and_mos
+
+! ---
+
+subroutine minimize_tc_orb_angles()
+
+  implicit none
+  logical          :: good_angles
+  integer          :: i
+  double precision :: thr_deg
+
+  good_angles = .False.
+  thr_deg = thr_degen_tc
+
+  call print_energy_and_mos()
+
+  print *, ' Minimizing the angles between the TC orbitals'
+  i = 1
+  do while (.not. good_angles)
+    print *, ' iteration = ', i
+    call routine_save_rotated_mos(thr_deg, good_angles)
+    thr_deg *= 10.d0
+    i += 1
+    if(i .gt. 100) then
+      print *, ' minimize_tc_orb_angles does not seem to converge ..'
+      print *, ' Something is weird in the tc orbitals ...'
+      print *, ' STOPPING'
+      stop
+    endif
+  enddo
+  print *, ' Converged ANGLES MINIMIZATION !!'
+
+  call print_angles_tc()
+  call print_energy_and_mos()
+
 end
 
-subroutine routine_save_rotated_mos(thr_deg,good_angles)
- implicit none
- double precision, intent(in) :: thr_deg
- logical, intent(out) :: good_angles
- good_angles = .False.
- integer :: i,j,k,n_degen_list,m,n,n_degen,ilast,ifirst
- double precision, allocatable :: mo_r_coef_good(:,:),mo_l_coef_good(:,:)
- allocate(mo_l_coef_good(ao_num, mo_num), mo_r_coef_good(ao_num,mo_num))
- double precision, allocatable :: mo_r_coef_new(:,:)
- double precision :: norm
- print*,'***************************************'
- print*,'***************************************'
- print*,'THRESHOLD FOR DEGENERACIES ::: ',thr_deg
- print*,'***************************************'
- print*,'***************************************'
- print*,'Starting with the following TC energy gradient :',grad_non_hermit
- mo_r_coef_good = mo_r_coef
- mo_l_coef_good = mo_l_coef
- allocate(mo_r_coef_new(ao_num, mo_num))
- mo_r_coef_new = mo_r_coef
- do i = 1, mo_num
-  norm = 1.d0/dsqrt(overlap_mo_r(i,i))
-  do j = 1, ao_num
-   mo_r_coef_new(j,i) *= norm
+! ---
+
+subroutine routine_save_rotated_mos(thr_deg, good_angles)
+
+  implicit none
+
+  double precision, intent(in)  :: thr_deg
+  logical,          intent(out) :: good_angles
+
+  integer                       :: i, j, k, n_degen_list, m, n, n_degen, ilast, ifirst
+  double precision              :: max_angle, norm
+  integer,          allocatable :: list_degen(:,:)
+  double precision, allocatable :: new_angles(:)
+  double precision, allocatable :: mo_r_coef_good(:,:), mo_l_coef_good(:,:)
+  double precision, allocatable :: mo_r_coef_new(:,:)
+  double precision, allocatable :: fock_diag(:),s_mat(:,:)
+  double precision, allocatable :: stmp(:,:), T(:,:), Snew(:,:), smat2(:,:)
+  double precision, allocatable :: mo_l_coef_tmp(:,:), mo_r_coef_tmp(:,:), mo_l_coef_new(:,:)
+
+  good_angles = .False.
+
+  allocate(mo_l_coef_good(ao_num, mo_num), mo_r_coef_good(ao_num,mo_num))
+
+  print *, ' ***************************************'
+  print *, ' ***************************************'
+  print *, ' THRESHOLD FOR DEGENERACIES ::: ', thr_deg
+  print *, ' ***************************************'
+  print *, ' ***************************************'
+  print *, ' Starting with the following TC energy gradient :', grad_non_hermit
+
+  mo_r_coef_good = mo_r_coef
+  mo_l_coef_good = mo_l_coef
+
+  allocate(mo_r_coef_new(ao_num, mo_num))
+  mo_r_coef_new = mo_r_coef
+  do i = 1, mo_num
+    norm = 1.d0/dsqrt(overlap_mo_r(i,i))
+    do j = 1, ao_num
+      mo_r_coef_new(j,i) *= norm
+    enddo
   enddo
- enddo
- double precision, allocatable :: fock_diag(:),s_mat(:,:)
- integer, allocatable :: list_degen(:,:)
- allocate(list_degen(mo_num,0:mo_num),s_mat(mo_num,mo_num),fock_diag(mo_num))
- do i = 1, mo_num
-  fock_diag(i) = Fock_matrix_tc_mo_tot(i,i)
- enddo
+
+  allocate(list_degen(mo_num,0:mo_num), s_mat(mo_num,mo_num), fock_diag(mo_num))
+  do i = 1, mo_num
+    fock_diag(i) = Fock_matrix_tc_mo_tot(i,i)
+  enddo
+
  ! compute the overlap between the left and rescaled right
- call build_s_matrix(ao_num,mo_num,mo_r_coef_new,mo_r_coef_new,ao_overlap,s_mat)
+  call build_s_matrix(ao_num, mo_num, mo_r_coef_new, mo_r_coef_new, ao_overlap, s_mat)
 ! call give_degen(fock_diag,mo_num,thr_deg,list_degen,n_degen_list)
- call give_degen_full_list(fock_diag,mo_num,thr_deg,list_degen,n_degen_list)
- print*,'fock_matrix_mo'
- do i = 1, mo_num
-  print*,i,fock_diag(i),angle_left_right(i)
- enddo
+  call give_degen_full_list(fock_diag, mo_num, thr_deg, list_degen, n_degen_list)
+  print *, ' fock_matrix_mo'
+  do i = 1, mo_num
+    print *, i, fock_diag(i), angle_left_right(i)
+  enddo
    
- do i = 1, n_degen_list
+  do i = 1, n_degen_list
 !  ifirst = list_degen(1,i)
 !  ilast  = list_degen(2,i)
 !  n_degen = ilast - ifirst +1
-  n_degen = list_degen(i,0)
-  double precision, allocatable :: stmp(:,:),T(:,:),Snew(:,:),smat2(:,:)
-  double precision, allocatable :: mo_l_coef_tmp(:,:),mo_r_coef_tmp(:,:),mo_l_coef_new(:,:)
-  allocate(stmp(n_degen,n_degen),smat2(n_degen,n_degen))
-  allocate(mo_r_coef_tmp(ao_num,n_degen),mo_l_coef_tmp(ao_num,n_degen),mo_l_coef_new(ao_num,n_degen))
-  allocate(T(n_degen,n_degen),Snew(n_degen,n_degen))
-  do j = 1, n_degen
-   mo_r_coef_tmp(1:ao_num,j) = mo_r_coef_new(1:ao_num,list_degen(i,j))
-   mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,list_degen(i,j))
-  enddo
-  ! Orthogonalization of right functions
-  print*,'Orthogonalization of RIGHT functions'
-  print*,'------------------------------------'
-  call orthog_functions(ao_num,n_degen,mo_r_coef_tmp,ao_overlap)
 
-  ! Orthogonalization of left functions
-  print*,'Orthogonalization of LEFT functions'
-  print*,'------------------------------------'
-  call orthog_functions(ao_num,n_degen,mo_l_coef_tmp,ao_overlap)
-  print*,'Overlap lef-right '
-  call build_s_matrix(ao_num,n_degen,mo_r_coef_tmp,mo_l_coef_tmp,ao_overlap,stmp)
-  do j = 1, n_degen
-   write(*,'(100(F8.4,X))')stmp(:,j)
+    n_degen = list_degen(i,0)
+
+    allocate(stmp(n_degen,n_degen), smat2(n_degen,n_degen))
+    allocate(mo_r_coef_tmp(ao_num,n_degen), mo_l_coef_tmp(ao_num,n_degen), mo_l_coef_new(ao_num,n_degen))
+    allocate(T(n_degen,n_degen), Snew(n_degen,n_degen))
+
+    do j = 1, n_degen
+      mo_r_coef_tmp(1:ao_num,j) = mo_r_coef_new(1:ao_num,list_degen(i,j))
+      mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,list_degen(i,j))
+    enddo
+    ! Orthogonalization of right functions
+    print *, ' Orthogonalization of RIGHT functions'
+    print *, ' ------------------------------------'
+    call orthog_functions(ao_num, n_degen, mo_r_coef_tmp, ao_overlap)
+  
+    ! Orthogonalization of left functions
+    print *, ' Orthogonalization of LEFT functions'
+    print *, ' ------------------------------------'
+    call orthog_functions(ao_num, n_degen, mo_l_coef_tmp, ao_overlap)
+
+    print *, ' Overlap lef-right '
+    call build_s_matrix(ao_num, n_degen, mo_r_coef_tmp, mo_l_coef_tmp, ao_overlap, stmp)
+    do j = 1, n_degen
+     write(*,'(100(F8.4,X))') stmp(:,j)
+    enddo
+    call build_s_matrix(ao_num, n_degen, mo_l_coef_tmp, mo_l_coef_tmp, ao_overlap, stmp)
+
+    !print*,'LEFT/LEFT OVERLAP '
+    !do j = 1, n_degen
+    ! write(*,'(100(F16.10,X))')stmp(:,j)
+    !enddo
+    call build_s_matrix(ao_num, n_degen, mo_r_coef_tmp, mo_r_coef_tmp, ao_overlap, stmp)
+    !print*,'RIGHT/RIGHT OVERLAP '
+    !do j = 1, n_degen
+    ! write(*,'(100(F16.10,X))')stmp(:,j)
+    !enddo
+
+    if(maxovl_tc) then
+      T    = 0.d0
+      Snew = 0.d0
+      call maxovl(n_degen, n_degen, stmp, T, Snew)
+    !print*,'overlap after'
+    !do j = 1, n_degen
+    ! write(*,'(100(F16.10,X))')Snew(:,j)
+    !enddo
+      call dgemm( 'N', 'N', ao_num, n_degen, n_degen, 1.d0                  &
+                , mo_l_coef_tmp, size(mo_l_coef_tmp, 1), T(1,1), size(T, 1) &
+                , 0.d0, mo_l_coef_new, size(mo_l_coef_new, 1) )
+     call build_s_matrix(ao_num, n_degen, mo_l_coef_new, mo_r_coef_tmp, ao_overlap, stmp)
+    !print*,'Overlap test'
+    !do j = 1, n_degen
+    ! write(*,'(100(F16.10,X))')stmp(:,j)
+    !enddo
+    else 
+      mo_l_coef_new = mo_l_coef_tmp
+    endif
+
+    call impose_weighted_biorthog_svd(ao_num, n_degen, ao_overlap, mo_l_coef_new, mo_r_coef_tmp)
+
+    !call build_s_matrix(ao_num, n_degen, mo_l_coef_new, mo_r_coef_tmp, ao_overlap, stmp)
+    !print*,'LAST OVERLAP '
+    !do j = 1, n_degen
+    ! write(*,'(100(F16.10,X))')stmp(:,j)
+    !enddo
+    !call build_s_matrix(ao_num, n_degen, mo_l_coef_new, mo_l_coef_new, ao_overlap, stmp)
+    !print*,'LEFT OVERLAP '
+    !do j = 1, n_degen
+    ! write(*,'(100(F16.10,X))')stmp(:,j)
+    !enddo
+    !call build_s_matrix(ao_num, n_degen, mo_r_coef_tmp, mo_r_coef_tmp, ao_overlap, stmp)
+    !print*,'RIGHT OVERLAP '
+    !do j = 1, n_degen
+    ! write(*,'(100(F16.10,X))')stmp(:,j)
+    !enddo
+    do j = 1, n_degen
+!!!   mo_l_coef_good(1:ao_num,j+ifirst-1) = mo_l_coef_new(1:ao_num,j)
+!!!   mo_r_coef_good(1:ao_num,j+ifirst-1) = mo_r_coef_tmp(1:ao_num,j)
+      mo_l_coef_good(1:ao_num,list_degen(i,j)) = mo_l_coef_new(1:ao_num,j)
+      mo_r_coef_good(1:ao_num,list_degen(i,j)) = mo_r_coef_tmp(1:ao_num,j)
+    enddo
+
+    deallocate(stmp, smat2)
+    deallocate(mo_r_coef_tmp, mo_l_coef_tmp, mo_l_coef_new)
+    deallocate(T, Snew)
   enddo
-  call build_s_matrix(ao_num,n_degen,mo_l_coef_tmp,mo_l_coef_tmp,ao_overlap,stmp)
+
+  !allocate(stmp(mo_num, mo_num))
+  !call build_s_matrix(ao_num, mo_num, mo_l_coef_good, mo_r_coef_good, ao_overlap, stmp)
+  !print*,'LEFT/RIGHT OVERLAP '
+  !do j = 1, mo_num
+  ! write(*,'(100(F16.10,X))')stmp(:,j)
+  !enddo
+  !call build_s_matrix(ao_num, mo_num, mo_l_coef_good, mo_l_coef_good, ao_overlap, stmp)
   !print*,'LEFT/LEFT OVERLAP '
-  !do j = 1, n_degen
+  !do j = 1, mo_num
   ! write(*,'(100(F16.10,X))')stmp(:,j)
   !enddo
-  call build_s_matrix(ao_num,n_degen,mo_r_coef_tmp,mo_r_coef_tmp,ao_overlap,stmp)
+  !call build_s_matrix(ao_num, mo_num, mo_r_coef_good, mo_r_coef_good, ao_overlap, stmp)
   !print*,'RIGHT/RIGHT OVERLAP '
-  !do j = 1, n_degen
+  !do j = 1, mo_num
   ! write(*,'(100(F16.10,X))')stmp(:,j)
   !enddo
-  if(maxovl_tc)then
-   T    = 0.d0
-   Snew = 0.d0
-   call maxovl(n_degen, n_degen, stmp, T, Snew)
-  !print*,'overlap after'
-  !do j = 1, n_degen
-  ! write(*,'(100(F16.10,X))')Snew(:,j)
-  !enddo
-   call dgemm( 'N', 'N', ao_num, n_degen, n_degen, 1.d0               &
-             , mo_l_coef_tmp, size(mo_l_coef_tmp, 1), T(1,1), size(T, 1) &
-             , 0.d0, mo_l_coef_new, size(mo_l_coef_new, 1) )
-   call build_s_matrix(ao_num,n_degen,mo_l_coef_new,mo_r_coef_tmp,ao_overlap,stmp)
-  !print*,'Overlap test'
-  !do j = 1, n_degen
-  ! write(*,'(100(F16.10,X))')stmp(:,j)
-  !enddo
-  else 
-   mo_l_coef_new = mo_l_coef_tmp
-  endif
-  call impose_biorthog_svd_overlap(ao_num, n_degen, ao_overlap, mo_l_coef_new, mo_r_coef_tmp)
-  call build_s_matrix(ao_num,n_degen,mo_l_coef_new,mo_r_coef_tmp,ao_overlap,stmp)
- !print*,'LAST OVERLAP '
- !do j = 1, n_degen
- ! write(*,'(100(F16.10,X))')stmp(:,j)
- !enddo
-  call build_s_matrix(ao_num,n_degen,mo_l_coef_new,mo_l_coef_new,ao_overlap,stmp)
- !print*,'LEFT OVERLAP '
- !do j = 1, n_degen
- ! write(*,'(100(F16.10,X))')stmp(:,j)
- !enddo
-  call build_s_matrix(ao_num,n_degen,mo_r_coef_tmp,mo_r_coef_tmp,ao_overlap,stmp)
- !print*,'RIGHT OVERLAP '
- !do j = 1, n_degen
- ! write(*,'(100(F16.10,X))')stmp(:,j)
- !enddo
-  do j = 1, n_degen
-!   mo_l_coef_good(1:ao_num,j+ifirst-1) = mo_l_coef_new(1:ao_num,j)
-!   mo_r_coef_good(1:ao_num,j+ifirst-1) = mo_r_coef_tmp(1:ao_num,j)
-   mo_l_coef_good(1:ao_num,list_degen(i,j)) = mo_l_coef_new(1:ao_num,j)
-   mo_r_coef_good(1:ao_num,list_degen(i,j)) = mo_r_coef_tmp(1:ao_num,j)
-  enddo
-  deallocate(stmp,smat2)
-  deallocate(mo_r_coef_tmp,mo_l_coef_tmp,mo_l_coef_new)
-  deallocate(T,Snew)
- enddo
 
- allocate(stmp(mo_num, mo_num))
- call build_s_matrix(ao_num,mo_num,mo_l_coef_good,mo_r_coef_good,ao_overlap,stmp)
- !print*,'LEFT/RIGHT OVERLAP '
- !do j = 1, mo_num
- ! write(*,'(100(F16.10,X))')stmp(:,j)
- !enddo
- call build_s_matrix(ao_num,mo_num,mo_l_coef_good,mo_l_coef_good,ao_overlap,stmp)
- !print*,'LEFT/LEFT OVERLAP '
- !do j = 1, mo_num
- ! write(*,'(100(F16.10,X))')stmp(:,j)
- !enddo
- call build_s_matrix(ao_num,mo_num,mo_r_coef_good,mo_r_coef_good,ao_overlap,stmp)
- !print*,'RIGHT/RIGHT OVERLAP '
- !do j = 1, mo_num
- ! write(*,'(100(F16.10,X))')stmp(:,j)
- !enddo
   mo_r_coef = mo_r_coef_good
   mo_l_coef = mo_l_coef_good
   call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
   call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
   TOUCH mo_l_coef mo_r_coef
-  double precision, allocatable :: new_angles(:)
+
   allocate(new_angles(mo_num))
   new_angles(1:mo_num) = dabs(angle_left_right(1:mo_num))
-  double precision :: max_angle
   max_angle = maxval(new_angles)
   good_angles = max_angle.lt.45.d0
-  print*,'max_angle = ',max_angle
+  print *, ' max_angle = ', max_angle
   
 end
 
-subroutine build_s_matrix(m,n,C1,C2,overlap,smat)
- implicit none
- integer, intent(in) :: m,n
- double precision, intent(in) :: C1(m,n),C2(m,n),overlap(m,m)
- double precision, intent(out):: smat(n,n)
- integer :: i,j,k,l
- smat = 0.D0
+! ---
+
+subroutine build_s_matrix(m, n, C1, C2, overlap, smat)
+
+  implicit none
+  integer,          intent(in)  :: m, n
+  double precision, intent(in)  :: C1(m,n), C2(m,n), overlap(m,m)
+  double precision, intent(out) :: smat(n,n)
+  integer                       :: i, j, k, l
+
+  smat = 0.D0
   do i = 1, n
-   do j = 1, n
-    do k = 1, m
-     do l = 1, m
-      smat(i,j) += C1(k,i) * overlap(l,k) * C2(l,j) 
-     enddo
+    do j = 1, n
+      do k = 1, m
+        do l = 1, m
+          smat(i,j) += C1(k,i) * overlap(l,k) * C2(l,j) 
+        enddo
+      enddo
     enddo
-   enddo
   enddo
+
 end
 
+! ---
+
 subroutine orthog_functions(m,n,coef,overlap)
  implicit none
  integer, intent(in) :: m,n
@@ -217,45 +256,63 @@ subroutine orthog_functions(m,n,coef,overlap)
    coef(1,:m) *= 1.d0/dsqrt(stmp(j,j))
   enddo
   call build_s_matrix(m,n,coef,coef,overlap,stmp)
+
  !print*,'overlap after'
  !do j = 1, n
  ! write(*,'(100(F16.10,X))')stmp(:,j)
  !enddo
+
+ deallocate(stmp)
+
 end
 
-subroutine print_angles_tc
- implicit none
- integer :: i,j
- double precision :: left,right
- print*,'product of norms, angle between vectors'                                                                  
- do i = 1, mo_num
-  left = overlap_mo_l(i,i)
-  right = overlap_mo_r(i,i)
-!  print*,Fock_matrix_tc_mo_tot(i,i),left*right,angle_left_right(i)
-  print*,left*right,angle_left_right(i)
- enddo
-end
+! ---
 
-subroutine print_energy_and_mos
- implicit none
- integer :: i
-  print*,''
-  print*,'TC energy = ', TC_HF_energy
-  print*,'TC SCF energy gradient = ',grad_non_hermit
-  print*,'Max angle Left/right  = ',max_angle_left_right
-  if(max_angle_left_right.lt.45.d0)then
-   print*,'Maximum angle BELOW 45 degrees, everthing is OK !'
-  else if(max_angle_left_right.gt.45.d0.and.max_angle_left_right.lt.75.d0)then
-   print*,'Maximum angle between 45 and 75 degrees, this is not the best for TC-CI calculations ...'
-  else if(max_angle_left_right.gt.75.d0)then
-   print*,'Maximum angle between ABOVE 75 degrees, YOU WILL CERTAINLY FIND TROUBLES IN TC-CI calculations ...'
-  endif
-  print*,'Diag Fock elem, product of left/right norm, angle left/right '
+subroutine print_angles_tc()
+
+  implicit none
+  integer          :: i, j
+  double precision :: left, right
+
+  print *, ' product of norms, angle between vectors'                                                                  
   do i = 1, mo_num
-   write(*,'(I3,X,100(F16.10,X))')i,Fock_matrix_tc_mo_tot(i,i),overlap_mo_l(i,i)*overlap_mo_r(i,i),angle_left_right(i)
+    left  = overlap_mo_l(i,i)
+    right = overlap_mo_r(i,i)
+!  print*,Fock_matrix_tc_mo_tot(i,i),left*right,angle_left_right(i)
+    print *, left*right, angle_left_right(i)
   enddo
+
 end
 
+! ---
+
+subroutine print_energy_and_mos()
+
+  implicit none
+  integer :: i
+
+  print *, ' '
+  print *, ' TC energy = ', TC_HF_energy
+  print *, ' TC SCF energy gradient = ', grad_non_hermit
+  print *, ' Max angle Left/right   = ', max_angle_left_right
+
+  if(max_angle_left_right .lt. 45.d0) then
+    print *, ' Maximum angle BELOW 45 degrees, everthing is OK !'
+  else if(max_angle_left_right .gt. 45.d0 .and. max_angle_left_right .lt. 75.d0) then
+    print *, ' Maximum angle between 45 and 75 degrees, this is not the best for TC-CI calculations ...'
+  else if(max_angle_left_right .gt. 75.d0) then
+    print *, ' Maximum angle between ABOVE 75 degrees, YOU WILL CERTAINLY FIND TROUBLES IN TC-CI calculations ...'
+  endif
+
+  print *, ' Diag Fock elem, product of left/right norm, angle left/right '
+  do i = 1, mo_num
+    write(*, '(I3,X,100(F16.10,X))') i, Fock_matrix_tc_mo_tot(i,i), overlap_mo_l(i,i)*overlap_mo_r(i,i), angle_left_right(i)
+  enddo
+
+end
+
+! ---
+
 subroutine sort_by_tc_fock
  implicit none 
  integer, allocatable :: iorder(:)
@@ -276,3 +333,4 @@ subroutine sort_by_tc_fock
  touch mo_l_coef mo_r_coef
 
 end
+
diff --git a/src/tools/print_he_energy.irp.f b/src/tools/print_he_energy.irp.f
new file mode 100644
index 00000000..321b6a92
--- /dev/null
+++ b/src/tools/print_he_energy.irp.f
@@ -0,0 +1,330 @@
+
+! ---
+
+program print_he_energy
+
+  implicit none
+
+  !call print_energy1()
+  call print_energy2()
+
+end 
+
+! ---
+
+subroutine print_energy1()
+
+  implicit none
+  integer                    :: i, j, k, l
+  double precision           :: e, n, e_tmp, n_tmp
+  double precision, external :: ao_two_e_integral
+
+  e = 0.d0
+  n = 0.d0
+
+  ! --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
+
+  ! < phi_1 phi_1 | h1 | phi_1 phi_1 > 
+
+  e_tmp = 0.d0
+  n_tmp = 0.d0
+  do i = 1, ao_num
+    do j = 1, ao_num
+      e_tmp += mo_coef(i,1) * ao_one_e_integrals(i,j) * mo_coef(j,1)
+      n_tmp += mo_coef(i,1) * ao_overlap(i,j)         * mo_coef(j,1)
+    enddo
+  enddo
+
+  e += e_tmp * n_tmp
+
+  ! ---
+
+  ! < phi_1 phi_1 | h2 | phi_1 phi_1 > 
+
+  e_tmp = 0.d0
+  n_tmp = 0.d0
+  do i = 1, ao_num
+    do j = 1, ao_num
+      n_tmp += mo_coef(i,1) * ao_overlap(i,j)         * mo_coef(j,1)
+      e_tmp += mo_coef(i,1) * ao_one_e_integrals(i,j) * mo_coef(j,1)
+    enddo
+  enddo
+
+  e += e_tmp * n_tmp
+
+  ! ---
+
+  ! --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
+
+  ! ---
+
+  do i = 1, ao_num
+    do j = 1, ao_num
+      do k = 1, ao_num
+        do l = 1, ao_num
+
+          ! ( phi_1 phi_1 | phi_1 phi_1 )
+          e += mo_coef(i,1) * mo_coef(j,1) * ao_two_e_integral(i,j,k,l) * mo_coef(k,1) * mo_coef(l,1)
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  ! ---
+
+  ! --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
+
+  ! ---
+
+  ! < phi_1 phi_1 | phi_1 phi_1 >
+  e_tmp = 0.d0
+  n_tmp = 0.d0
+  do i = 1, ao_num
+    do j = 1, ao_num
+      e_tmp += mo_coef(i,1) * ao_overlap(i,j) * mo_coef(j,1)
+      n_tmp += mo_coef(i,1) * ao_overlap(i,j) * mo_coef(j,1)
+    enddo
+  enddo
+
+  n += e_tmp * n_tmp
+
+  ! ---
+
+  ! --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
+
+  e = e / n 
+  print *, ' energy = ', e
+
+end subroutine print_energy1
+
+! ---
+
+subroutine print_energy2()
+
+  implicit none
+  integer                    :: i, j, k, l
+  double precision           :: e, n, e_tmp, n_tmp
+  double precision, external :: ao_two_e_integral
+
+  e = 0.d0
+  n = 0.d0
+
+  ! --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
+
+  ! < phi_1 phi_2 | h1 | phi_1 phi_2 > 
+
+  e_tmp = 0.d0
+  n_tmp = 0.d0
+  do i = 1, ao_num
+    do j = 1, ao_num
+      e_tmp += mo_coef(i,1) * ao_one_e_integrals(i,j) * mo_coef(j,1)
+      n_tmp += mo_coef(i,2) * ao_overlap(i,j)         * mo_coef(j,2)
+    enddo
+  enddo
+
+  e += e_tmp * n_tmp
+
+  ! ---
+
+  ! < phi_1 phi_2 | h1 | phi_2 phi_1 > 
+
+  e_tmp = 0.d0
+  n_tmp = 0.d0
+  do i = 1, ao_num
+    do j = 1, ao_num
+      e_tmp += mo_coef(i,1) * ao_one_e_integrals(i,j) * mo_coef(j,2)
+      n_tmp += mo_coef(i,2) * ao_overlap(i,j)         * mo_coef(j,1)
+    enddo
+  enddo
+
+  e -= e_tmp * n_tmp
+
+  ! ---
+
+  ! < phi_2 phi_1 | h1 | phi_1 phi_2 > 
+
+  e_tmp = 0.d0
+  n_tmp = 0.d0
+  do i = 1, ao_num
+    do j = 1, ao_num
+      e_tmp += mo_coef(i,2) * ao_one_e_integrals(i,j) * mo_coef(j,1)
+      n_tmp += mo_coef(i,1) * ao_overlap(i,j)         * mo_coef(j,2)
+    enddo
+  enddo
+
+  e -= e_tmp * n_tmp
+
+  ! ---
+
+  ! < phi_2 phi_1 | h1 | phi_2 phi_1 >
+
+  e_tmp = 0.d0
+  n_tmp = 0.d0
+  do i = 1, ao_num
+    do j = 1, ao_num
+      e_tmp += mo_coef(i,2) * ao_one_e_integrals(i,j) * mo_coef(j,2)
+      n_tmp += mo_coef(i,1) * ao_overlap(i,j)         * mo_coef(j,1)
+    enddo
+  enddo
+
+  e += e_tmp * n_tmp
+
+  ! ---
+
+  ! < phi_1 phi_2 | h2 | phi_1 phi_2 > 
+
+  e_tmp = 0.d0
+  n_tmp = 0.d0
+  do i = 1, ao_num
+    do j = 1, ao_num
+      n_tmp += mo_coef(i,1) * ao_overlap(i,j)         * mo_coef(j,1)
+      e_tmp += mo_coef(i,2) * ao_one_e_integrals(i,j) * mo_coef(j,2)
+    enddo
+  enddo
+
+  e += e_tmp * n_tmp
+
+  ! ---
+
+  ! < phi_1 phi_2 | h2 | phi_2 phi_1 >
+
+  e_tmp = 0.d0
+  n_tmp = 0.d0
+  do i = 1, ao_num
+    do j = 1, ao_num
+      n_tmp += mo_coef(i,1) * ao_overlap(i,j)         * mo_coef(j,2)
+      e_tmp += mo_coef(i,2) * ao_one_e_integrals(i,j) * mo_coef(j,1)
+    enddo
+  enddo
+
+  e -= e_tmp * n_tmp
+
+  ! ---
+
+  ! < phi_2 phi_1 | h2 | phi_1 phi_2 >
+
+  e_tmp = 0.d0
+  n_tmp = 0.d0
+  do i = 1, ao_num
+    do j = 1, ao_num
+      n_tmp += mo_coef(i,2) * ao_overlap(i,j)         * mo_coef(j,1)
+      e_tmp += mo_coef(i,1) * ao_one_e_integrals(i,j) * mo_coef(j,2)
+    enddo
+  enddo
+
+  e -= e_tmp * n_tmp
+
+  ! ---
+
+  ! < phi_2 phi_1 | h2 | phi_2 phi_1 > 
+
+  e_tmp = 0.d0
+  n_tmp = 0.d0
+  do i = 1, ao_num
+    do j = 1, ao_num
+      n_tmp += mo_coef(i,2) * ao_overlap(i,j)         * mo_coef(j,2)
+      e_tmp += mo_coef(i,1) * ao_one_e_integrals(i,j) * mo_coef(j,1)
+    enddo
+  enddo
+
+  e += e_tmp * n_tmp
+
+  ! ---
+
+  ! --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
+
+  ! ---
+
+  do i = 1, ao_num
+    do j = 1, ao_num
+      do k = 1, ao_num
+        do l = 1, ao_num
+
+          ! ( phi_1 phi_1 | phi_2 phi_2 )
+          e += mo_coef(i,1) * mo_coef(j,1) * ao_two_e_integral(i,j,k,l) * mo_coef(k,2) * mo_coef(l,2)
+
+          ! ( phi_1 phi_2 | phi_2 phi_1 )
+          e -= mo_coef(i,1) * mo_coef(j,2) * ao_two_e_integral(i,j,k,l) * mo_coef(k,2) * mo_coef(l,1)
+
+          ! ( phi_2 phi_1 | phi_1 phi_2 )
+          e -= mo_coef(i,2) * mo_coef(j,1) * ao_two_e_integral(i,j,k,l) * mo_coef(k,1) * mo_coef(l,2)
+
+          ! ( phi_2 phi_2 | phi_1 phi_1 )
+          e += mo_coef(i,2) * mo_coef(j,2) * ao_two_e_integral(i,j,k,l) * mo_coef(k,1) * mo_coef(l,1)
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  ! ---
+
+  ! --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
+
+  ! ---
+
+  ! < phi_1 phi_2 | phi_1 phi_2 >
+  e_tmp = 0.d0
+  n_tmp = 0.d0
+  do i = 1, ao_num
+    do j = 1, ao_num
+      e_tmp += mo_coef(i,1) * ao_overlap(i,j) * mo_coef(j,1)
+      n_tmp += mo_coef(i,2) * ao_overlap(i,j) * mo_coef(j,2)
+    enddo
+  enddo
+
+  n += e_tmp * n_tmp
+
+  ! ---
+
+  ! < phi_1 phi_2 | phi_2 phi_1 >
+  e_tmp = 0.d0
+  n_tmp = 0.d0
+  do i = 1, ao_num
+    do j = 1, ao_num
+      e_tmp += mo_coef(i,1) * ao_overlap(i,j) * mo_coef(j,2)
+      n_tmp += mo_coef(i,2) * ao_overlap(i,j) * mo_coef(j,1)
+    enddo
+  enddo
+
+  n -= e_tmp * n_tmp
+
+  ! ---
+
+  ! < phi_2 phi_1 | phi_1 phi_2 >
+  e_tmp = 0.d0
+  n_tmp = 0.d0
+  do i = 1, ao_num
+    do j = 1, ao_num
+      e_tmp += mo_coef(i,2) * ao_overlap(i,j) * mo_coef(j,1)
+      n_tmp += mo_coef(i,1) * ao_overlap(i,j) * mo_coef(j,2)
+    enddo
+  enddo
+
+  n -= e_tmp * n_tmp
+
+  ! ---
+
+  ! < phi_2 phi_1 | phi_2 phi_1 >
+  e_tmp = 0.d0
+  n_tmp = 0.d0
+  do i = 1, ao_num
+    do j = 1, ao_num
+      e_tmp += mo_coef(i,2) * ao_overlap(i,j) * mo_coef(j,2)
+      n_tmp += mo_coef(i,1) * ao_overlap(i,j) * mo_coef(j,1)
+    enddo
+  enddo
+
+  n += e_tmp * n_tmp
+
+  ! ---
+
+  ! --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
+
+  e = e / n 
+  print *, ' energy = ', e
+
+end subroutine print_energy2
+
+! ---
diff --git a/src/utils/block_diag_degen.irp.f b/src/utils/block_diag_degen.irp.f
index f6906e23..155eef9c 100644
--- a/src/utils/block_diag_degen.irp.f
+++ b/src/utils/block_diag_degen.irp.f
@@ -1,177 +1,218 @@
 
-subroutine diag_mat_per_fock_degen(fock_diag,mat_ref,n,thr_deg,leigvec,reigvec,eigval)
- implicit none
- integer, intent(in) :: n
- double precision, intent(in) :: fock_diag(n),mat_ref(n,n),thr_deg
- double precision, intent(out):: leigvec(n,n),reigvec(n,n),eigval(n)
- BEGIN_DOC
-! subroutine that diagonalizes a matrix mat_ref BY BLOCK
-!
-! the blocks are defined by the elements having the SAME DEGENERACIES in the entries "fock_diag"
-!
-! examples : all elements having degeneracy 1 in fock_diag (i.e. not being degenerated) will be treated together
-!
-!          : all elements having degeneracy 2 in fock_diag (i.e. two elements are equal) will be treated together
-!
-!          : all elements having degeneracy 3 in fock_diag (i.e. two elements are equal) will be treated together
-!
-! etc... the advantage is to guarentee no spurious mixing because of numerical problems. 
- END_DOC
- double precision, allocatable :: leigvec_unsrtd(:,:),reigvec_unsrtd(:,:),eigval_unsrtd(:)
- integer, allocatable :: list_degen(:,:),list_same_degen(:)
- integer, allocatable :: iorder(:),list_degen_sorted(:)
- integer :: n_degen_list,n_degen,size_mat,i,j,k,icount,m,index_degen
- integer :: ii,jj,i_good,j_good,n_real
- double precision, allocatable :: mat_tmp(:,:),eigval_tmp(:),leigvec_tmp(:,:),reigvec_tmp(:,:)
+subroutine diag_mat_per_fock_degen(fock_diag, mat_ref, n, thr_d, thr_nd, thr_deg, leigvec, reigvec, eigval)
 
- allocate(leigvec_unsrtd(n,n),reigvec_unsrtd(n,n),eigval_unsrtd(n))
- leigvec_unsrtd = 0.d0
- reigvec_unsrtd = 0.d0
- eigval_unsrtd  = 0.d0
 
- allocate(list_degen(n,0:n))
+  BEGIN_DOC
+  !
+  ! subroutine that diagonalizes a matrix mat_ref BY BLOCK
+  !
+  ! the blocks are defined by the elements having the SAME DEGENERACIES in the entries "fock_diag"
+  !
+  ! examples : all elements having degeneracy 1 in fock_diag (i.e. not being degenerated) will be treated together
+  !
+  !          : all elements having degeneracy 2 in fock_diag (i.e. two elements are equal) will be treated together
+  !
+  !          : all elements having degeneracy 3 in fock_diag (i.e. two elements are equal) will be treated together
+  !
+  ! etc... the advantage is to guarentee no spurious mixing because of numerical problems. 
+  !
+  END_DOC
 
- ! obtain degeneracies 
- call give_degen_full_list(fock_diag,n,thr_deg,list_degen,n_degen_list)
- allocate(iorder(n_degen_list),list_degen_sorted(n_degen_list))
- do i =1, n_degen_list
-  n_degen = list_degen(i,0) 
-  list_degen_sorted(i) = n_degen
-  iorder(i)=i
- enddo
- ! sort by number of degeneracies 
- call isort(list_degen_sorted,iorder,n_degen_list)
- integer :: icount_eigval
- logical, allocatable :: is_ok(:)
- allocate(is_ok(n_degen_list))
- is_ok = .True.
- icount_eigval = 0
- ! loop over degeneracies 
- do i = 1, n_degen_list
-  if(.not.is_ok(i))cycle
-  is_ok(i) = .False.
-  n_degen = list_degen_sorted(i)
-  print*,'diagonalizing for n_degen = ',n_degen
-  k = 1
-  ! group all the entries having the same degeneracies 
-!  do while (list_degen_sorted(i+k)==n_degen)
-  do m = i+1,n_degen_list
-   if(list_degen_sorted(m)==n_degen)then
-    is_ok(i+k)=.False.
-    k += 1
-   endif
+  implicit none
+  integer,          intent(in)  :: n
+  double precision, intent(in)  :: fock_diag(n), mat_ref(n,n), thr_d, thr_nd, thr_deg
+  double precision, intent(out) :: leigvec(n,n), reigvec(n,n), eigval(n)
+
+  integer                       :: n_degen_list, n_degen,size_mat, i, j, k, icount, m, index_degen
+  integer                       :: ii, jj, i_good, j_good, n_real
+  integer                       :: icount_eigval
+  logical,          allocatable :: is_ok(:)
+  integer,          allocatable :: list_degen(:,:), list_same_degen(:)
+  integer,          allocatable :: iorder(:), list_degen_sorted(:)
+  double precision, allocatable :: leigvec_unsrtd(:,:), reigvec_unsrtd(:,:), eigval_unsrtd(:)
+  double precision, allocatable :: mat_tmp(:,:), eigval_tmp(:), leigvec_tmp(:,:), reigvec_tmp(:,:)
+
+  allocate(leigvec_unsrtd(n,n), reigvec_unsrtd(n,n), eigval_unsrtd(n))
+  leigvec_unsrtd = 0.d0
+  reigvec_unsrtd = 0.d0
+  eigval_unsrtd  = 0.d0
+
+  ! obtain degeneracies 
+  allocate(list_degen(n,0:n))
+  call give_degen_full_list(fock_diag, n, thr_deg, list_degen, n_degen_list)
+
+  allocate(iorder(n_degen_list), list_degen_sorted(n_degen_list))
+  do i = 1, n_degen_list
+    n_degen = list_degen(i,0) 
+    list_degen_sorted(i) = n_degen
+    iorder(i) = i
   enddo
-  print*,'number of identical degeneracies = ',k
-  size_mat = k*n_degen  
-  print*,'size_mat = ',size_mat
-  allocate(mat_tmp(size_mat,size_mat),list_same_degen(size_mat))
-  allocate(eigval_tmp(size_mat),leigvec_tmp(size_mat,size_mat),reigvec_tmp(size_mat,size_mat))
-  ! group all the elements sharing the same degeneracy
-  icount = 0
-  do j = 1, k ! jth set of degeneracy
-   index_degen = iorder(i+j-1)
-   do m = 1, n_degen
-    icount += 1
-    list_same_degen(icount) = list_degen(index_degen,m)
-   enddo
+
+  ! sort by number of degeneracies 
+  call isort(list_degen_sorted, iorder, n_degen_list)
+
+  allocate(is_ok(n_degen_list))
+  is_ok = .True.
+  icount_eigval = 0
+
+  ! loop over degeneracies 
+  do i = 1, n_degen_list
+    if(.not.is_ok(i)) cycle
+
+    is_ok(i) = .False.
+    n_degen  = list_degen_sorted(i)
+
+    print *, ' diagonalizing for n_degen = ', n_degen
+
+    k = 1
+
+   ! group all the entries having the same degeneracies 
+!!  do while (list_degen_sorted(i+k)==n_degen)
+    do m = i+1, n_degen_list
+      if(list_degen_sorted(m)==n_degen) then
+        is_ok(i+k) = .False.
+        k += 1
+      endif
+    enddo
+
+    print *, ' number of identical degeneracies = ', k
+    size_mat = k*n_degen  
+    print *, ' size_mat = ', size_mat
+    allocate(mat_tmp(size_mat,size_mat), list_same_degen(size_mat))
+    allocate(eigval_tmp(size_mat), leigvec_tmp(size_mat,size_mat), reigvec_tmp(size_mat,size_mat))
+    ! group all the elements sharing the same degeneracy
+    icount = 0
+    do j = 1, k ! jth set of degeneracy
+      index_degen = iorder(i+j-1)
+      do m = 1, n_degen
+        icount += 1
+        list_same_degen(icount) = list_degen(index_degen,m)
+      enddo
+    enddo
+
+    print *, ' list of elements '
+    do icount = 1, size_mat
+      print *, icount, list_same_degen(icount)
+    enddo
+
+    ! you copy subset of matrix elements having all the same degeneracy in mat_tmp
+    do ii = 1, size_mat
+      i_good = list_same_degen(ii)
+      do jj = 1, size_mat
+        j_good = list_same_degen(jj)
+        mat_tmp(jj,ii) = mat_ref(j_good,i_good)
+      enddo
+    enddo
+
+    call non_hrmt_bieig( size_mat, mat_tmp, thr_d, thr_nd &
+                       , leigvec_tmp, reigvec_tmp         & 
+                       , n_real, eigval_tmp )
+
+    do ii = 1, size_mat
+      icount_eigval += 1
+      eigval_unsrtd(icount_eigval) = eigval_tmp(ii) ! copy eigenvalues 
+      do jj = 1, size_mat ! copy the eigenvectors 
+        j_good = list_same_degen(jj)
+        leigvec_unsrtd(j_good,icount_eigval) = leigvec_tmp(jj,ii)
+        reigvec_unsrtd(j_good,icount_eigval) = reigvec_tmp(jj,ii)
+      enddo
+    enddo
+
+    deallocate(mat_tmp, list_same_degen)
+    deallocate(eigval_tmp, leigvec_tmp, reigvec_tmp)
   enddo
-  print*,'list of elements '
-  do icount = 1, size_mat
-   print*,icount,list_same_degen(icount)
-  enddo
-  ! you copy subset of matrix elements having all the same degeneracy in mat_tmp
-  do ii = 1, size_mat
-   i_good = list_same_degen(ii)
-   do jj = 1, size_mat
-    j_good = list_same_degen(jj)
-    mat_tmp(jj,ii) = mat_ref(j_good,i_good)
-   enddo
-  enddo
-  call non_hrmt_bieig( size_mat, mat_tmp&
-                     , leigvec_tmp, reigvec_tmp& 
-                     , n_real, eigval_tmp)
-  do ii = 1, size_mat
-   icount_eigval += 1
-   eigval_unsrtd(icount_eigval) = eigval_tmp(ii) ! copy eigenvalues 
-   do jj = 1, size_mat ! copy the eigenvectors 
-    j_good = list_same_degen(jj)
-    leigvec_unsrtd(j_good,icount_eigval) = leigvec_tmp(jj,ii)
-    reigvec_unsrtd(j_good,icount_eigval) = reigvec_tmp(jj,ii)
-   enddo
-  enddo
-  deallocate(mat_tmp,list_same_degen)
-  deallocate(eigval_tmp,leigvec_tmp,reigvec_tmp)
- enddo
- if(icount_eigval.ne.n)then
-  print*,'pb !! (icount_eigval.ne.n)'
-  print*,'icount_eigval,n',icount_eigval,n
-  stop
- endif
+
+  if(icount_eigval .ne. n) then
+    print *, ' pb !! (icount_eigval.ne.n)'
+    print *, ' icount_eigval,n', icount_eigval, n
+    stop
+  endif
  
- deallocate(iorder)
- allocate(iorder(n))
- do i = 1, n
-  iorder(i) = i
- enddo
- call dsort(eigval_unsrtd,iorder,n)
- do i = 1, n
-  print*,'sorted eigenvalues '
-  i_good = iorder(i)
-  eigval(i) = eigval_unsrtd(i)
-  print*,'i,eigval(i) = ',i,eigval(i)
-  do j = 1, n
-   leigvec(j,i) = leigvec_unsrtd(j,i_good)
-   reigvec(j,i) = reigvec_unsrtd(j,i_good)
+  deallocate(iorder)
+  allocate(iorder(n))
+  do i = 1, n
+    iorder(i) = i
   enddo
- enddo
+  call dsort(eigval_unsrtd, iorder, n)
+
+  do i = 1, n
+    print*,'sorted eigenvalues '
+    i_good = iorder(i)
+    eigval(i) = eigval_unsrtd(i)
+    print*,'i,eigval(i) = ',i,eigval(i)
+    do j = 1, n
+      leigvec(j,i) = leigvec_unsrtd(j,i_good)
+      reigvec(j,i) = reigvec_unsrtd(j,i_good)
+    enddo
+  enddo
+
+  deallocate(leigvec_unsrtd, reigvec_unsrtd, eigval_unsrtd)
+  deallocate(list_degen)
+  deallocate(iorder, list_degen_sorted)
+  deallocate(is_ok)
+
 end
 
+! ---
+
 subroutine give_degen_full_list(A,n,thr,list_degen,n_degen_list)
- implicit none
- BEGIN_DOC
- ! you enter with an array A(n) and spits out all the elements degenerated up to thr
- !
- ! the elements of A(n) DON'T HAVE TO BE SORTED IN THE ENTRANCE: TOTALLY GENERAL 
- !
- ! list_degen(i,0) = number of degenerate entries 
- !
- ! list_degen(i,1) = index of the first degenerate entry
- !
- ! list_degen(i,2:list_degen(i,0)) = list of all other dengenerate entries 
- !
- ! if list_degen(i,0) == 1 it means that there is no degeneracy for that element
- END_DOC
- double precision,intent(in) :: A(n)
- double precision,intent(in)    :: thr
- integer, intent(in) :: n
- integer, intent(out) :: list_degen(n,0:n),n_degen_list
- logical, allocatable :: is_ok(:)
- allocate(is_ok(n))
- integer :: i,j,icount,icheck
- n_degen_list = 0
- is_ok = .True.
- do i = 1, n
-  if(.not.is_ok(i))cycle
-  n_degen_list +=1
-  is_ok(i) = .False.
-  list_degen(n_degen_list,1) = i
-  icount = 1
-  do j = i+1, n
-   if(dabs(A(i)-A(j)).lt.thr.and.is_ok(j))then
-    is_ok(j) = .False.
-    icount += 1
-    list_degen(n_degen_list,icount) = j
-   endif
+
+  BEGIN_DOC
+  ! you enter with an array A(n) and spits out all the elements degenerated up to thr
+  !
+  ! the elements of A(n) DON'T HAVE TO BE SORTED IN THE ENTRANCE: TOTALLY GENERAL 
+  !
+  ! list_degen(i,0) = number of degenerate entries 
+  !
+  ! list_degen(i,1) = index of the first degenerate entry
+  !
+  ! list_degen(i,2:list_degen(i,0)) = list of all other dengenerate entries 
+  !
+  ! if list_degen(i,0) == 1 it means that there is no degeneracy for that element
+  END_DOC
+
+  implicit none
+
+  double precision, intent(in)  :: A(n)
+  double precision, intent(in)  :: thr
+  integer,          intent(in)  :: n
+  integer,          intent(out) :: list_degen(n,0:n), n_degen_list
+  integer                       :: i, j, icount, icheck
+  logical, allocatable          :: is_ok(:)
+
+
+  allocate(is_ok(n))
+  n_degen_list = 0
+  is_ok = .True.
+  do i = 1, n
+    if(.not.is_ok(i)) cycle
+    n_degen_list +=1
+    is_ok(i) = .False.
+    list_degen(n_degen_list,1) = i
+    icount = 1
+    do j = i+1, n
+      if(dabs(A(i)-A(j)).lt.thr.and.is_ok(j)) then
+        is_ok(j) = .False.
+        icount += 1
+        list_degen(n_degen_list,icount) = j
+      endif
+    enddo
+
+    list_degen(n_degen_list,0) = icount
   enddo
-  list_degen(n_degen_list,0) = icount
- enddo
- icheck = 0
- do i = 1, n_degen_list
-  icheck += list_degen(i,0)
- enddo
- if(icheck.ne.n)then
-  print*,'pb ! :: icheck.ne.n'
-  print*,icheck,n
-  stop
- endif
+
+  icheck = 0
+  do i = 1, n_degen_list
+    icheck += list_degen(i,0)
+  enddo
+
+  if(icheck.ne.n)then
+    print *, ' pb ! :: icheck.ne.n'
+    print *, icheck, n
+    stop
+  endif
+
 end
+
+! ---
+