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 de3bc23c..0739640c 100644
--- a/src/non_hermit_dav/lapack_diag_non_hermit.irp.f
+++ b/src/non_hermit_dav/lapack_diag_non_hermit.irp.f
@@ -1281,10 +1281,10 @@ subroutine impose_orthog_svd_overlap(n, m, C,overlap)
             , 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
+! print *, ' eigenvec overlap bef SVD: '
+! do i = 1, m
+!   write(*, '(1000(F16.10,X))') S(i,:)
+! enddo
 
   ! ---
  
@@ -1340,10 +1340,10 @@ subroutine impose_orthog_svd_overlap(n, m, C,overlap)
             , 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
+! print *, ' eigenvec overlap aft SVD: '
+! do i = 1, m
+!   write(*, '(1000(F16.10,X))') S(i,:)
+! enddo
 
   deallocate(S)
 
@@ -2516,7 +2516,7 @@ subroutine impose_biorthog_svd_overlap(n, m, overlap, L, R)
 
   print *, ' overlap bef SVD: '
   do i = 1, m
-    write(*, '(1000(F16.10,X))') S(i,:)
+    write(*, '(1000(F25.16,X))') S(i,:)
   enddo
 
   ! ---
@@ -2530,6 +2530,7 @@ subroutine impose_biorthog_svd_overlap(n, m, overlap, L, R)
   threshold               = 1.d-6
   num_linear_dependencies = 0
   do i = 1, m
+    print*,'D(i) = ',D(i)
     if(abs(D(i)) <= threshold) then
       D(i) = 0.d0
       num_linear_dependencies += 1
@@ -2585,11 +2586,18 @@ 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) )
+  ! 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), R, size(R, 1) &
+            , L, size(L, 1), Stmp, size(Stmp, 1) &
             , 0.d0, S, size(S, 1) )
 
-  print *, ' overlap aft SVD: '
+  print *, ' overlap aft SVD with overlap: '
   do i = 1, m
     write(*, '(1000(F16.10,X))') S(i,:)
   enddo
diff --git a/src/tc_scf/fock_tc.irp.f b/src/tc_scf/fock_tc.irp.f
index 42b429d5..4907f2c8 100644
--- a/src/tc_scf/fock_tc.irp.f
+++ b/src/tc_scf/fock_tc.irp.f
@@ -131,3 +131,30 @@ END_PROVIDER
 
 ! ---
 
+ BEGIN_PROVIDER [ double precision, grad_non_hermit_left]
+&BEGIN_PROVIDER [ double precision, grad_non_hermit_right]
+&BEGIN_PROVIDER [ double precision, grad_non_hermit]
+ implicit none
+  integer :: i, k
+  grad_non_hermit_left = 0.d0
+  grad_non_hermit_right = 0.d0
+  do i = 1, elec_beta_num ! doc --> SOMO
+    do k = elec_beta_num+1, elec_alpha_num
+      grad_non_hermit_left+= dabs(Fock_matrix_tc_mo_tot(k,i))
+      grad_non_hermit_right+= dabs(Fock_matrix_tc_mo_tot(i,k))
+    enddo
+  enddo
+  do i = 1, elec_beta_num ! doc --> virt 
+    do k = elec_alpha_num+1, mo_num
+      grad_non_hermit_left+= dabs(Fock_matrix_tc_mo_tot(k,i))
+      grad_non_hermit_right+= dabs(Fock_matrix_tc_mo_tot(i,k))
+    enddo
+  enddo
+  do i = elec_beta_num+1, elec_alpha_num ! SOMO --> virt 
+    do k = elec_alpha_num+1, mo_num
+      grad_non_hermit_left+= dabs(Fock_matrix_tc_mo_tot(k,i))
+      grad_non_hermit_right+= dabs(Fock_matrix_tc_mo_tot(i,k))
+    enddo
+  enddo
+ grad_non_hermit = grad_non_hermit_left + grad_non_hermit_right
+END_PROVIDER 
diff --git a/src/tc_scf/minimize_tc_angles.irp.f b/src/tc_scf/minimize_tc_angles.irp.f
new file mode 100644
index 00000000..258fa114
--- /dev/null
+++ b/src/tc_scf/minimize_tc_angles.irp.f
@@ -0,0 +1,10 @@
+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 minimize_tc_orb_angles
+end
+
diff --git a/src/tc_scf/print_angle_tc_orb.irp.f b/src/tc_scf/print_angle_tc_orb.irp.f
index 51c45848..09260395 100644
--- a/src/tc_scf/print_angle_tc_orb.irp.f
+++ b/src/tc_scf/print_angle_tc_orb.irp.f
@@ -5,16 +5,5 @@ program print_angles
 !  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
- double precision :: left,right
- print*,'energy,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)
- enddo
+  call print_angles_tc
 end
diff --git a/src/tc_scf/routines_rotates.irp.f b/src/tc_scf/routines_rotates.irp.f
index 0603cc77..b7c2f2d7 100644
--- a/src/tc_scf/routines_rotates.irp.f
+++ b/src/tc_scf/routines_rotates.irp.f
@@ -1,10 +1,44 @@
-subroutine routine_save_rotated_mos
+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
+ 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
+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))
@@ -19,22 +53,19 @@ subroutine routine_save_rotated_mos
  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)
+  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 give_degen(fock_diag,mo_num,thr_degen_tc,list_degen,n_degen_list)
+ call give_degen(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
- 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)
+!  if(ifirst.ne.12)cycle
   ilast  = list_degen(2,i)
   n_degen = ilast - ifirst +1
   print*,'ifirst,n_degen = ',ifirst,n_degen
@@ -48,16 +79,29 @@ subroutine routine_save_rotated_mos
    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'
+  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'
-  call orthog_functions(ao_num,n_degen,mo_r_coef_tmp,ao_overlap)
+  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
@@ -77,6 +121,16 @@ subroutine routine_save_rotated_mos
   else 
    mo_l_coef_new = mo_l_coef_tmp
   endif
+  call build_s_matrix(ao_num,n_degen,mo_l_coef_new,mo_l_coef_new,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
   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 '
@@ -131,6 +185,13 @@ subroutine routine_save_rotated_mos
   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
+  
 end
 
 subroutine build_s_matrix(m,n,C1,C2,overlap,smat)
@@ -175,3 +236,28 @@ subroutine orthog_functions(m,n,coef,overlap)
    write(*,'(100(F16.10,X))')stmp(:,j)
   enddo
 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*,'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
diff --git a/src/tc_scf/tc_scf.irp.f b/src/tc_scf/tc_scf.irp.f
index bf46fc5d..bd7284ee 100644
--- a/src/tc_scf/tc_scf.irp.f
+++ b/src/tc_scf/tc_scf.irp.f
@@ -19,7 +19,7 @@ program tc_scf
   !call orthonormalize_mos
 
   call routine_scf()
-  call routine_save_rotated_mos
+  call minimize_tc_orb_angles
   call print_energy_and_mos
 
 
@@ -126,6 +126,7 @@ subroutine routine_scf()
       print*,'***'
       e_delta = dabs( TC_HF_energy - e_save )
       print*, 'it, delta E = ', it, e_delta
+      print*, 'it, gradient= ',grad_non_hermit_right
       e_save    = TC_HF_energy
       mo_l_coef = fock_tc_leigvec_ao
       mo_r_coef = fock_tc_reigvec_ao
@@ -181,13 +182,3 @@ end subroutine routine_scf
 
 ! ---
 
-subroutine print_energy_and_mos
- implicit none
- integer :: i
-  print*,'Energy converged !'
-  print*,'Final TC energy = ', TC_HF_energy
-  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
diff --git a/src/tc_scf/test_rotate_2.irp.f b/src/tc_scf/test_rotate_2.irp.f
deleted file mode 100644
index b8458925..00000000
--- a/src/tc_scf/test_rotate_2.irp.f
+++ /dev/null
@@ -1,187 +0,0 @@
-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
- 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
- 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))
- print*,'l coef'
- do i = 1, mo_num
-  write(*,'(100(F8.4,X))')mo_l_coef_good(:,i)
- enddo
- print*,'r coef'
- do i = 1, mo_num
-  write(*,'(100(F8.4,X))')mo_r_coef_good(:,i)
- enddo
- 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
-  call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef_good)
-  call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef_good)
-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