diff --git a/src/non_hermit_dav/biorthog.irp.f b/src/non_hermit_dav/biorthog.irp.f
index d6add005..12764f44 100644
--- a/src/non_hermit_dav/biorthog.irp.f
+++ b/src/non_hermit_dav/biorthog.irp.f
@@ -329,6 +329,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)
     if(dabs(WI(i)) .lt. thr) then
       n_good += 1
     else
@@ -392,8 +393,8 @@ subroutine non_hrmt_bieig(n, A, leigvec, reigvec, n_real_eigv, eigval)
   ! -------------------------------------------------------------------------------------
   !                               check bi-orthogonality
 
-  thr_d  = 1d-8
-  thr_nd = 1d-8
+  thr_d  = 1d-10 ! -7
+  thr_nd = 1d-10 ! -7
 
   allocate( S(n_real_eigv,n_real_eigv) )
   call check_biorthog(n, n_real_eigv, leigvec, reigvec, accu_d, accu_nd, S, .false.)
@@ -422,8 +423,8 @@ subroutine non_hrmt_bieig(n, A, leigvec, reigvec, n_real_eigv, eigval)
 
     ! ---
 
-    !call impose_orthog_degen_eigvec(n, eigval, reigvec)
-    !call impose_orthog_degen_eigvec(n, eigval, leigvec)
+!   call impose_orthog_degen_eigvec(n, eigval, reigvec)
+!   call impose_orthog_degen_eigvec(n, eigval, leigvec)
 
     call impose_biorthog_degen_eigvec(n, eigval, leigvec, reigvec)
 
@@ -1040,6 +1041,72 @@ subroutine split_matrix_degen(aw,n,shift)
 
 end
 
+subroutine give_degen(a,n,shift,list_degen,n_degen_list)
+ implicit none
+ BEGIN_DOC
+ ! returns n_degen_list :: the number of degenerated SET of elements (i.e. with |A(i)-A(i+1)| below shift)
+ !
+ ! for each of these sets, list_degen(1,i) = first degenerate element of the set i, 
+ !
+ !                         list_degen(2,i) = last degenerate element of the set i.
+ END_DOC
+ double precision,intent(in) :: A(n)
+ double precision,intent(in)    :: shift
+ integer, intent(in) :: n
+ integer, intent(out) :: list_degen(2,n),n_degen_list
+ integer :: i,j,n_degen,k
+ logical :: keep_on
+ double precision,allocatable :: Aw(:)
+ list_degen = -1
+ allocate(Aw(n))
+ Aw = A
+ i=1
+ k = 0
+ do while(i.lt.n)
+  if(dabs(Aw(i)-Aw(i+1)).lt.shift)then
+   k+=1
+   j=1
+   list_degen(1,k) = i
+   keep_on = .True.
+   do while(keep_on)
+    if(i+j.gt.n)then
+     keep_on = .False.
+     exit
+    endif
+    if(dabs(Aw(i)-Aw(i+j)).lt.shift)then
+     j+=1
+    else
+     keep_on=.False.
+     exit
+    endif
+   enddo
+   n_degen = j
+   list_degen(2,k) = list_degen(1,k)-1 + n_degen
+   j=0
+   keep_on = .True.
+   do while(keep_on)
+    if(i+j+1.gt.n)then
+     keep_on = .False.
+     exit
+    endif
+    if(dabs(Aw(i+j)-Aw(i+j+1)).lt.shift)then
+     Aw(i+j) += (j-n_degen/2) * shift
+     j+=1
+    else 
+     keep_on = .False.
+     exit
+    endif
+   enddo
+   Aw(i+n_degen-1) += (n_degen-1-n_degen/2) * shift
+   i+=n_degen
+  else 
+   i+=1
+  endif
+ enddo
+ n_degen_list = k
+
+end
+
 subroutine cancel_small_elmts(aw,n,shift)
  implicit none
  BEGIN_DOC
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 2482a0d6..de3bc23c 100644
--- a/src/non_hermit_dav/lapack_diag_non_hermit.irp.f
+++ b/src/non_hermit_dav/lapack_diag_non_hermit.irp.f
@@ -1254,6 +1254,105 @@ end subroutine impose_orthog_svd
 
 ! ---
 
+subroutine impose_orthog_svd_overlap(n, m, C,overlap)
+
+  implicit none
+
+  integer,          intent(in)    :: n, m
+  double precision, intent(in   ) :: overlap(n,n)
+  double precision, intent(inout) :: C(n,m)
+
+  integer                         :: i, j, num_linear_dependencies
+  double precision                :: threshold
+  double precision, allocatable   :: S(:,:), tmp(:,:),Stmp(:,:)
+  double precision, allocatable   :: U(:,:), Vt(:,:), D(:)
+
+  print *, ' apply SVD to orthogonalize vectors'
+
+  ! ---
+
+  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), C, size(C, 1) &
+            , 0.d0, Stmp, size(Stmp, 1) )
+  call dgemm( 'T', 'N', m, m, n, 1.d0      &
+            , C, size(C, 1), Stmp, size(Stmp, 1) &
+            , 0.d0, S, size(S, 1) )
+
+  print *, ' eigenvec overlap bef SVD: '
+  do i = 1, m
+    write(*, '(1000(F16.10,X))') S(i,:)
+  enddo
+
+  ! ---
+ 
+  allocate(U(m,m), Vt(m,m), D(m))
+
+  call svd(S, m, U, m, D, Vt, m, m, m)
+
+  deallocate(S)
+
+  threshold               = 1.d-6
+  num_linear_dependencies = 0
+  do i = 1, m
+    if(abs(D(i)) <= threshold) then
+      D(i) = 0.d0
+      num_linear_dependencies += 1
+    else
+      ASSERT (D(i) > 0.d0)
+      D(i) = 1.d0 / dsqrt(D(i))
+    endif
+  enddo
+  if(num_linear_dependencies > 0) then
+    write(*,*) ' linear dependencies = ', num_linear_dependencies
+    write(*,*) ' m                   = ', m
+    stop
+  endif
+
+  ! ---
+
+  allocate(tmp(n,m))
+
+  ! tmp <-- C x U
+  call dgemm( 'N', 'N', n, m, m, 1.d0      &
+            , C, size(C, 1), U, size(U, 1) &
+            , 0.d0, tmp, size(tmp, 1) )
+
+  deallocate(U, Vt)
+
+  ! C <-- tmp x sigma^-0.5
+  do j = 1, m
+    do i = 1, n
+      C(i,j) = tmp(i,j) * D(j)
+    enddo
+  enddo
+
+  deallocate(D, tmp)
+
+  ! ---
+
+  allocate(S(m,m))
+
+  ! S = C.T x C
+  call dgemm( 'T', 'N', m, m, n, 1.d0      &
+            , C, size(C, 1), C, size(C, 1) &
+            , 0.d0, S, size(S, 1) )
+
+  print *, ' eigenvec overlap aft SVD: '
+  do i = 1, m
+    write(*, '(1000(F16.10,X))') S(i,:)
+  enddo
+
+  deallocate(S)
+
+  ! ---
+
+end subroutine impose_orthog_svd
+
+! ---
+
 subroutine impose_orthog_GramSchmidt(n, m, C)
 
   implicit none
@@ -2389,3 +2488,116 @@ end subroutine impose_biorthog_svd
 
 ! ---
 
+
+subroutine impose_biorthog_svd_overlap(n, m, overlap, L, R)
+
+  implicit none
+
+  integer,          intent(in)    :: n, m
+  double precision, intent(in)    :: overlap(n,n)
+  double precision, intent(inout) :: L(n,m), R(n,m)
+
+  integer                         :: i, j, num_linear_dependencies
+  double precision                :: threshold
+  double precision, allocatable   :: S(:,:), tmp(:,:),Stmp(:,:)
+  double precision, allocatable   :: U(:,:), V(:,:), Vt(:,:), D(:)
+
+  ! ---
+
+  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) &
+            , 0.d0, Stmp, size(Stmp, 1) )
+  call dgemm( 'T', 'N', m, m, n, 1.d0      &
+            , L, size(L, 1), Stmp, size(Stmp, 1) &
+            , 0.d0, S, size(S, 1) )
+
+  print *, ' overlap bef SVD: '
+  do i = 1, m
+    write(*, '(1000(F16.10,X))') S(i,:)
+  enddo
+
+  ! ---
+ 
+  allocate(U(m,m), Vt(m,m), D(m))
+
+  call svd(S, m, U, m, D, Vt, m, m, m)
+
+  deallocate(S)
+
+  threshold               = 1.d-6
+  num_linear_dependencies = 0
+  do i = 1, m
+    if(abs(D(i)) <= threshold) then
+      D(i) = 0.d0
+      num_linear_dependencies += 1
+    else
+      ASSERT (D(i) > 0.d0)
+      D(i) = 1.d0 / dsqrt(D(i))
+    endif
+  enddo
+  if(num_linear_dependencies > 0) then
+    write(*,*) ' linear dependencies = ', num_linear_dependencies
+    write(*,*) ' m                   = ', m
+    stop
+  endif
+
+  allocate(V(m,m))
+  do i = 1, m
+    do j = 1, m
+      V(j,i) = Vt(i,j)
+    enddo
+  enddo
+  deallocate(Vt)
+
+  ! ---
+
+  allocate(tmp(n,m))
+
+  ! tmp <-- R x V
+  call dgemm( 'N', 'N', n, m, m, 1.d0      &
+            , R, size(R, 1), V, size(V, 1) &
+            , 0.d0, tmp, size(tmp, 1) )
+  deallocate(V)
+  ! R <-- tmp x sigma^-0.5
+  do j = 1, m
+    do i = 1, n
+      R(i,j) = tmp(i,j) * D(j)
+    enddo
+  enddo
+
+  ! tmp <-- L x U
+  call dgemm( 'N', 'N', n, m, m, 1.d0      &
+            , L, size(L, 1), U, size(U, 1) &
+            , 0.d0, tmp, size(tmp, 1) )
+  deallocate(U)
+  ! L <-- tmp x sigma^-0.5
+  do j = 1, m
+    do i = 1, n
+      L(i,j) = tmp(i,j) * D(j)
+    enddo
+  enddo
+
+  deallocate(D, tmp)
+
+  ! ---
+
+  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) )
+
+  print *, ' overlap aft SVD: '
+  do i = 1, m
+    write(*, '(1000(F16.10,X))') S(i,:)
+  enddo
+  deallocate(S)
+
+  ! ---
+
+end subroutine impose_biorthog_svd
+
+! ---
+
diff --git a/src/some_mu_of_r/print_mu_av_tc.irp.f b/src/some_mu_of_r/print_mu_av_tc.irp.f
new file mode 100644
index 00000000..48e06a25
--- /dev/null
+++ b/src/some_mu_of_r/print_mu_av_tc.irp.f
@@ -0,0 +1,67 @@
+program print_mu_av_tc
+ implicit none
+ read_wf = .True.
+ touch read_wf
+   print*,'average_mu_lda       = ',average_mu_lda
+!   print*,'average_mu_rs        = ',average_mu_rs 
+   print*,'average_mu_rs_c      = ',average_mu_rs_c
+   print*,'average_mu_rs_c_lda  = ',average_mu_rs_c_lda
+   call plot_mu_tc
+end
+
+
+
+subroutine plot_mu_tc
+ implicit none
+ integer :: i,nx,m
+ double precision :: xmin,xmax,dx
+ double precision :: r(3)
+ double precision :: sqpi
+ double precision :: rho_a_hf, rho_b_hf, g0,rho_hf
+ double precision :: rs,grad_n,grad_n_a(3), grad_n_b(3)
+ double precision :: g0_UEG_mu_inf
+ double precision :: cst_rs,alpha_rs,mu_rs, mu_rs_c, mu_lda, mu_grad_n 
+ character*(128) :: output
+ integer :: i_unit_output,getUnitAndOpen
+ if (no_core_density)then
+  output=trim(ezfio_filename)//'.tc_mu_fc'
+ else
+  output=trim(ezfio_filename)//'.tc_mu'
+ endif
+ i_unit_output = getUnitAndOpen(output,'w')
+
+ nx = 5000
+ xmin = -10.d0
+ xmax = 10.d0
+ dx   = (xmax - xmin)/dble(nx)
+ r = 0.d0
+ r(3) = xmin
+
+ sqpi = dsqrt(dacos(-1.d0))
+ cst_rs   = (4.d0 * dacos(-1.d0)/3.d0)**(-1.d0/3.d0)
+ alpha_rs = 2.d0 * dsqrt((9.d0 * dacos(-1.d0)/4.d0)**(-1.d0/3.d0)) / sqpi
+
+  write(i_unit_output,*)'#       r(3)              rho_hf          mu_rs           mu_rs_c           mu_lda         mu_grad_n'
+ do i = 1, nx
+  rho_a_hf = 0.d0
+  grad_n   = 0.d0
+  call density_and_grad_alpha_beta(r,rho_a_hf,rho_b_hf, grad_n_a, grad_n_b)
+  do m = 1, 3
+   grad_n += grad_n_a(m)*grad_n_a(m) + grad_n_b(m)*grad_n_b(m) + 2.d0 * grad_n_a(m) * grad_n_b(m)
+  enddo
+  rho_hf = rho_a_hf + rho_b_hf
+  grad_n = dsqrt(grad_n)
+  grad_n = grad_n/(4.d0 * rho_hf)
+  rs = cst_rs * rho_hf**(-1.d0/3.d0)
+  g0 = g0_UEG_mu_inf(rho_a_hf,rho_b_hf)
+
+  mu_rs     = 1.d0/rs 
+  mu_rs_c   = alpha_rs/dsqrt(rs) 
+  mu_lda    =  - 1.d0 / (dlog(2.d0 * g0) * sqpi)  
+  mu_grad_n = grad_n 
+  write(i_unit_output,'(100(F16.10,X))')r(3),rho_hf,mu_rs,mu_rs_c,mu_lda,mu_grad_n
+  r(3) += dx
+ enddo
+
+
+end
diff --git a/src/some_mu_of_r/some_mu_of_r.irp.f b/src/some_mu_of_r/some_mu_of_r.irp.f
deleted file mode 100644
index 67d364e9..00000000
--- a/src/some_mu_of_r/some_mu_of_r.irp.f
+++ /dev/null
@@ -1,7 +0,0 @@
-program some_mu_of_r
-  implicit none
-  BEGIN_DOC
-! TODO : Put the documentation of the program here
-  END_DOC
-  print *, 'Hello world'
-end
diff --git a/src/tc_scf/diago_bi_ort_tcfock.irp.f b/src/tc_scf/diago_bi_ort_tcfock.irp.f
index 69083e33..00558de6 100644
--- a/src/tc_scf/diago_bi_ort_tcfock.irp.f
+++ b/src/tc_scf/diago_bi_ort_tcfock.irp.f
@@ -18,13 +18,13 @@
 
   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          &
+                    , 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 &
+!                       , 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 & 
diff --git a/src/tc_scf/test_rotate_2.irp.f b/src/tc_scf/test_rotate_2.irp.f
new file mode 100644
index 00000000..5246d824
--- /dev/null
+++ b/src/tc_scf/test_rotate_2.irp.f
@@ -0,0 +1,175 @@
+program print_angles
+ implicit none
+  my_grid_becke  = .True.
+!  my_n_pt_r_grid = 30
+!  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = 10 ! small grid for quick debug
+  my_n_pt_a_grid = 14 ! small grid for quick debug
+  call routine
+end
+subroutine routine
+ implicit none
+ 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
+ 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
+ double precision, allocatable :: fock_diag(:),s_mat(:,:)
+ integer, allocatable :: list_degen(:,:)
+ allocate(list_degen(2,mo_num),s_mat(mo_num,mo_num),fock_diag(mo_num))
+ do i = 1, mo_num
+  fock_diag(i) = fock_matrix_mo(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 give_degen(fock_diag,mo_num,thr_degen_tc,list_degen,n_degen_list)
+ print*,'fock_matrix_mo'
+ do i = 1, mo_num
+  print*,i,fock_diag(i),angle_left_right(i)
+ enddo
+ print*,'Overlap '
+ do i = 1, mo_num
+  write(*,'(I2,X,100(F8.4,X))')i,s_mat(:,i)
+ enddo
+   
+ do i = 1, n_degen_list
+  ifirst = list_degen(1,i)
+  ilast  = list_degen(2,i)
+  n_degen = ilast - ifirst +1
+  print*,'ifirst,n_degen = ',ifirst,n_degen
+  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,j+ifirst-1)
+   mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,j+ifirst-1)
+  enddo
+  ! Orthogonalization of right functions
+  print*,'Orthogonalization of right functions'
+  call orthog_functions(ao_num,n_degen,mo_r_coef_tmp,ao_overlap)
+  ! Orthogonalization of left functions
+  print*,'Orthogonalization of left functions'
+  call orthog_functions(ao_num,n_degen,mo_r_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
+  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
+!  mo_l_coef_new = 0.D0
+!  do j = 1, n_degen
+!   do k = 1, n_degen
+!    do m = 1, ao_num
+!     mo_l_coef_new(m,j) += T(k,j) * mo_l_coef_tmp(m,k) 
+!    enddo 
+!   enddo
+!  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
+  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)
+  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,n_degen,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,n_degen,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,n_degen,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
+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
+  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
+    enddo
+   enddo
+  enddo
+end
+
+subroutine orthog_functions(m,n,coef,overlap)
+ implicit none
+ integer, intent(in) :: m,n
+ double precision, intent(in)    :: overlap(m,m)
+ double precision, intent(inout) :: coef(m,n)
+ double precision, allocatable :: stmp(:,:)
+ integer :: j
+ allocate(stmp(n,n))
+  call build_s_matrix(m,n,coef,coef,overlap,stmp)
+  print*,'overlap before'
+  do j = 1, n
+   write(*,'(100(F16.10,X))')stmp(:,j)
+  enddo
+  call impose_orthog_svd_overlap(m, n, coef,overlap)
+  call build_s_matrix(m,n,coef,coef,overlap,stmp)
+  do j = 1, n
+   coef(1,:m) *= 1.d0/dsqrt(stmp(j,j))
+  enddo
+  print*,'overlap after'
+  call build_s_matrix(m,n,coef,coef,overlap,stmp)
+  do j = 1, n
+   write(*,'(100(F16.10,X))')stmp(:,j)
+  enddo
+end
diff --git a/src/utils/loc.f b/src/utils/loc.f
index 32d2cd03..6e5d7345 100644
--- a/src/utils/loc.f
+++ b/src/utils/loc.f
@@ -12,8 +12,10 @@ C     W: new overlap matrix
 C
 C 
       implicit real*8(a-h,o-y),logical*1(z)
-      parameter (id1=700)
-      dimension s(id1,id1),t(id1,id1),w(id1,id1)
+!      parameter (id1=700)
+!      dimension s(id1,id1),t(id1,id1),w(id1,id1)
+      double precision, intent(inout) :: s(n,n)
+      double precision, intent(out) :: t(n,n),w(n,n)
       data small/1.d-6/
 
       zprt=.true.