diff --git a/ocaml/Units.ml b/ocaml/Units.ml
index ab0e944c..7a7cb39e 100644
--- a/ocaml/Units.ml
+++ b/ocaml/Units.ml
@@ -4,8 +4,8 @@ type units =
 | Angstrom
 ;;
 
-let angstrom_to_bohr = 1. /. 0.52917721092
-let bohr_to_angstrom = 0.52917721092
+let angstrom_to_bohr = 1. /. 0.52917721067121
+let bohr_to_angstrom = 0.52917721067121
 ;;
 
 
diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f
index b6e959d7..8c6b3875 100644
--- a/src/ao_two_e_ints/two_e_integrals.irp.f
+++ b/src/ao_two_e_ints/two_e_integrals.irp.f
@@ -436,7 +436,7 @@ BEGIN_PROVIDER [ double precision, ao_two_e_integral_schwartz,(ao_num,ao_num)  ]
       !$OMP SCHEDULE(dynamic)
   do i=1,ao_num
     do k=1,i
-      ao_two_e_integral_schwartz(i,k) = dsqrt(ao_two_e_integral(i,k,i,k))
+      ao_two_e_integral_schwartz(i,k) = dsqrt(ao_two_e_integral(i,i,k,k))
       ao_two_e_integral_schwartz(k,i) = ao_two_e_integral_schwartz(i,k)
     enddo
   enddo
diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f
index f0bd8847..31f27e1d 100644
--- a/src/cipsi/pt2_stoch_routines.irp.f
+++ b/src/cipsi/pt2_stoch_routines.irp.f
@@ -287,9 +287,9 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
       call omp_set_nested(.false.)
 
 
-      print '(A)', '========== ================= =========== =============== =============== ================='
-      print '(A)', ' Samples        Energy        Stat. Err     Variance          Norm^2        Seconds      '
-      print '(A)', '========== ================= =========== =============== =============== ================='
+      print '(A)', '========== ======================= ===================== ===================== ==========='
+      print '(A)', ' Samples          Energy                Variance               Norm^2          Seconds'
+      print '(A)', '========== ======================= ===================== ===================== ==========='
 
       PROVIDE global_selection_buffer
 
@@ -312,26 +312,30 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
       !$OMP END PARALLEL
       call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'pt2')
 
-      print '(A)', '========== ================= =========== =============== =============== ================='
+      print '(A)', '========== ======================= ===================== ===================== ==========='
 
     do k=1,N_states
       pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate)
     enddo
-!    ! The overlap is not exactly zero because of the guiding function.
-!    ! Remove the bias
-!    do k=1,pt2_stoch_istate-1
-!      pt2_overlap(k,pt2_stoch_istate) -= pt2_data % overlap(k,pt2_stoch_istate)
-!    enddo
-print *, 'Overlap before orthogonalization'
-print *, pt2_overlap(:,pt2_stoch_istate)
-print *, 'Overlap after orthogonalization'
-print *, pt2_overlap(pt2_stoch_istate,:)
-print *, '-------'
     SOFT_TOUCH pt2_overlap
 
     enddo
     FREE pt2_stoch_istate
 
+    ! Symmetrize overlap
+    do j=2,N_states
+     do i=1,j-1
+       pt2_overlap(i,j) = 0.5d0 * (pt2_overlap(i,j) + pt2_overlap(j,i))
+       pt2_overlap(j,i) = pt2_overlap(i,j)
+     enddo
+    enddo
+
+    print *, 'Overlap of perturbed states:'
+    do k=1,N_states
+      print *, pt2_overlap(k,:)
+    enddo
+    print *, '-------'
+
     if (N_in > 0) then
       b%cur = min(N_in,b%cur)
       if (s2_eig) then
@@ -529,7 +533,14 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_
 
           if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then
             time1 = time
-            print '(G10.3, 2X, F16.10, 2X, G10.3, 2X, G14.6, 2X, G14.6, 2X, F10.4)', c, avg+E, eqt, avg2, avg3(pt2_stoch_istate), time-time0
+            print '(I10, X, F12.6, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.4)', c, &
+              pt2_data     % pt2(pt2_stoch_istate) +E, &
+              pt2_data_err % pt2(pt2_stoch_istate), &
+              pt2_data     % variance(pt2_stoch_istate), &
+              pt2_data_err % variance(pt2_stoch_istate), &
+              pt2_data     % overlap(pt2_stoch_istate,pt2_stoch_istate), &
+              pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate), &
+              time-time0
             if (stop_now .or. (                                      &
                   (do_exit .and. (dabs(pt2_data_err % pt2(pt2_stoch_istate)) /    &
                   (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) ) <= relative_error))) ) then
diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f
index c325020f..6c1adfde 100644
--- a/src/cipsi/selection.irp.f
+++ b/src/cipsi/selection.irp.f
@@ -1,4 +1,3 @@
-
 use bitmasks
 
 BEGIN_PROVIDER [ double precision, pt2_match_weight, (N_states) ]
@@ -144,7 +143,6 @@ BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ]
 
 END_PROVIDER
 
-
 subroutine get_mask_phase(det1, pm, Nint)
   use bitmasks
   implicit none
@@ -288,6 +286,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
   PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
   PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
   PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp
+  PROVIDE banned_excitation
 
   monoAdo = .true.
   monoBdo = .true.
@@ -607,7 +606,8 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
 
           h2 = hole_list(i2,s2)
           call apply_hole(pmask, s2,h2, mask, ok, N_int)
-          banned = .false.
+          banned(:,:,1) = banned_excitation(:,:)
+          banned(:,:,2) = banned_excitation(:,:)
           do j=1,mo_num
             bannedOrb(j, 1) = .true.
             bannedOrb(j, 2) = .true.
@@ -674,7 +674,8 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
   logical :: ok
   integer :: s1, s2, p1, p2, ib, j, istate, jstate
   integer(bit_kind) :: mask(N_int, 2), det(N_int, 2)
-  double precision :: e_pert, delta_E, val, Hii, w, tmp, alpha_h_psi, coef(N_states)
+  double precision :: e_pert(N_states), coef(N_states)
+  double precision :: delta_E, val, Hii, w, tmp, alpha_h_psi
   double precision, external :: diag_H_mat_elem_fock
   double precision :: E_shift
 
@@ -682,7 +683,12 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
   double precision, allocatable :: values(:)
   integer, allocatable          :: keys(:,:)
   integer                       :: nkeys
-
+  double precision :: s_weight(N_states,N_states)
+  do jstate=1,N_states
+    do istate=1,N_states
+      s_weight(istate,jstate) = dsqrt(selection_weight(istate)*selection_weight(jstate))
+    enddo
+  enddo
 
   if(sp == 3) then
     s1 = 1
@@ -763,23 +769,59 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
 !      call occ_pattern_of_det(det,occ,N_int)
 !      call occ_pattern_to_dets_size(occ,n,elec_alpha_num,N_int)
 
+      e_pert = 0.d0
+      coef = 0.d0
+      logical :: do_diag
+      do_diag = .False.
 
       do istate=1,N_states
         delta_E = E0(istate) - Hii + E_shift
         alpha_h_psi = mat(istate, p1, p2)
+        if (alpha_h_psi == 0.d0) cycle
+
         val = alpha_h_psi + alpha_h_psi
         tmp = dsqrt(delta_E * delta_E + val * val)
         if (delta_E < 0.d0) then
             tmp = -tmp
         endif
-        e_pert = 0.5d0 * (tmp - delta_E)
+        e_pert(istate) = 0.5d0 * (tmp - delta_E)
         if (dabs(alpha_h_psi) > 1.d-4) then
-          coef(istate) = e_pert / alpha_h_psi
+          coef(istate) = e_pert(istate) / alpha_h_psi
         else
           coef(istate) = alpha_h_psi / delta_E
         endif
       enddo
 
+      do_diag = sum(dabs(coef)) > 0.001d0
+
+      double precision :: eigvalues(N_states+1)
+      double precision :: work(1+6*(N_states+1)+2*(N_states+1)**2)
+      integer :: iwork(3+5*(N_states+1)), info, k ,n
+
+      if (do_diag) then
+        double precision :: pt2_matrix(N_states+1,N_states+1)
+        pt2_matrix(N_states+1,N_states+1) = Hii+E_shift
+        do istate=1,N_states
+          pt2_matrix(:,istate) = 0.d0
+          pt2_matrix(istate,istate) = E0(istate)
+          pt2_matrix(istate,N_states+1) = mat(istate,p1,p2)
+          pt2_matrix(N_states+1,istate) = mat(istate,p1,p2)
+        enddo
+
+        call DSYEVD( 'V', 'U', N_states+1, pt2_matrix, N_states+1, eigvalues, &
+                     work, size(work), iwork, size(iwork), info )
+        if (info /= 0) then
+          print *, 'error in '//irp_here
+          stop -1
+        endif
+        pt2_matrix = dabs(pt2_matrix)
+        iwork(1:N_states+1) = maxloc(pt2_matrix,DIM=1)
+        do k=1,N_states
+          e_pert(iwork(k)) = eigvalues(k) - E0(iwork(k))
+        enddo
+      endif
+
+
 !      ! Gram-Schmidt using input overlap matrix
 !      do istate=1,N_states
 !        do jstate=1,istate-1
@@ -791,14 +833,13 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
       do istate=1, N_states
 
         alpha_h_psi = mat(istate, p1, p2)
-        e_pert = coef(istate) * alpha_h_psi
 
         do jstate=1,N_states
           pt2_data % overlap(jstate,istate) += coef(jstate) * alpha_h_psi
         enddo
 
         pt2_data % variance(istate) += alpha_h_psi * alpha_h_psi
-        pt2_data % pt2(istate)      += e_pert
+        pt2_data % pt2(istate)      += e_pert(istate)
 
 !!!DEBUG
 !        delta_E = E0(istate) - Hii + E_shift
@@ -823,14 +864,29 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
 
           case(5)
             ! Variance selection
-            w = w - alpha_h_psi * alpha_h_psi * selection_weight(istate)
+!            w = w - alpha_h_psi * alpha_h_psi * s_weight(istate,istate)
+            w = min(w, - alpha_h_psi * alpha_h_psi * s_weight(istate,istate))
+!            do jstate=1,N_states
+!              if (istate == jstate) cycle
+!              w = w + dabs(alpha_h_psi*mat(jstate,p1,p2)) * s_weight(istate,jstate)
+!            enddo
 
           case(6)
-            w = w - coef(istate) * coef(istate) * selection_weight(istate)
+!            w = w - coef(istate) * coef(istate) * s_weight(istate,istate)
+            w = min(w,- coef(istate) * coef(istate) * s_weight(istate,istate))
+!            do jstate=1,N_states
+!              if (istate == jstate) cycle
+!              w = w + dabs(coef(istate)*coef(jstate)) * s_weight(istate,jstate)
+!            enddo
 
           case default
             ! Energy selection
-            w = w + e_pert * selection_weight(istate)
+!            w = w + e_pert(istate) * s_weight(istate,istate)
+            w = min(w, e_pert(istate) * s_weight(istate,istate))
+!            do jstate=1,N_states
+!              if (istate == jstate) cycle
+!              w = w + dabs(X(istate)*X(jstate)) * s_weight(istate,jstate)
+!            enddo
 
         end select
       end do
diff --git a/src/cipsi/zmq_selection.irp.f b/src/cipsi/zmq_selection.irp.f
index 9ce345cf..448b409e 100644
--- a/src/cipsi/zmq_selection.irp.f
+++ b/src/cipsi/zmq_selection.irp.f
@@ -141,6 +141,12 @@ subroutine ZMQ_selection(N_in, pt2_data)
   enddo
 
   pt2_overlap(:,:) = pt2_data % overlap(:,:)
+
+  print *, 'Overlap of perturbed states:'
+  do l=1,N_states
+    print *, pt2_overlap(l,:)
+  enddo
+  print *, '-------'
   SOFT_TOUCH pt2_overlap
   call update_pt2_and_variance_weights(pt2_data, N_states)
 
diff --git a/src/davidson/EZFIO.cfg b/src/davidson/EZFIO.cfg
index f820730c..4343cc1e 100644
--- a/src/davidson/EZFIO.cfg
+++ b/src/davidson/EZFIO.cfg
@@ -8,7 +8,7 @@ default: 1.e-10
 type: logical
 doc: Thresholds of Davidson's algorithm is set to E(rPT2)*threshold_davidson_from_pt2
 interface: ezfio,provider,ocaml
-default: true
+default: false
 
 [n_states_diag]
 type: States_number
diff --git a/src/determinants/EZFIO.cfg b/src/determinants/EZFIO.cfg
index ef00080b..662c6fbb 100644
--- a/src/determinants/EZFIO.cfg
+++ b/src/determinants/EZFIO.cfg
@@ -44,7 +44,7 @@ default: 2
 type: integer
 doc: Weight used in the selection. 0: input state-average weight, 1: 1./(c_0^2), 2: rPT2 matching, 3: variance matching, 4: variance and rPT2 matching, 5: variance minimization and matching, 6: CI coefficients 7: input state-average multiplied by variance and rPT2 matching 8: input state-average multiplied by rPT2 matching 9: input state-average multiplied by variance matching
 interface: ezfio,provider,ocaml
-default: 2
+default: 1
 
 [threshold_generators]
 type: Threshold
diff --git a/src/mo_two_e_ints/map_integrals.irp.f b/src/mo_two_e_ints/map_integrals.irp.f
index 661add2e..2221e0c9 100644
--- a/src/mo_two_e_ints/map_integrals.irp.f
+++ b/src/mo_two_e_ints/map_integrals.irp.f
@@ -99,6 +99,10 @@ double precision function get_two_e_integral(i,j,k,l,map)
   type(map_type), intent(inout)  :: map
   real(integral_kind)            :: tmp
   PROVIDE mo_two_e_integrals_in_map mo_integrals_cache
+  if (banned_excitation(i,k) .or. banned_excitation(j,l)) then
+    get_two_e_integral = 0.d0
+    return
+  endif
   ii = l-mo_integrals_cache_min
   ii = ior(ii, k-mo_integrals_cache_min)
   ii = ior(ii, j-mo_integrals_cache_min)
@@ -159,6 +163,11 @@ subroutine get_mo_two_e_integrals(j,k,l,sze,out_val,map)
 !  return
 !DEBUG
 
+  out_val(1:sze) = 0.d0
+  if (banned_excitation(j,l)) then
+    return
+  endif
+
   ii0 = l-mo_integrals_cache_min
   ii0 = ior(ii0, k-mo_integrals_cache_min)
   ii0 = ior(ii0, j-mo_integrals_cache_min)
@@ -172,6 +181,7 @@ subroutine get_mo_two_e_integrals(j,k,l,sze,out_val,map)
   q = q+shiftr(s*s-s,1)
 
   do i=1,sze
+    if (banned_excitation(i,k)) cycle
     ii = ior(ii0, i-mo_integrals_cache_min)
     if (iand(ii, -128) == 0) then
       ii_8 = ior( shiftl(ii0_8,7), int(i,8)-mo_integrals_cache_min_8)
@@ -272,6 +282,29 @@ subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map)
 
 end
 
+BEGIN_PROVIDER [ logical, banned_excitation, (mo_num,mo_num) ]
+ implicit none
+ use map_module
+ BEGIN_DOC
+ ! If true, the excitation is banned in the selection. Useful with local MOs.
+ END_DOC
+ banned_excitation = .False.
+ integer :: i,j
+ integer(key_kind)              :: idx
+ double precision :: tmp
+! double precision :: buffer(mo_num)
+ do j=1,mo_num
+   do i=1,j-1
+    call two_e_integrals_index(i,j,j,i,idx)
+    !DIR$ FORCEINLINE
+    call map_get(mo_integrals_map,idx,tmp)
+    banned_excitation(i,j) = dabs(tmp) < 1.d-15
+    banned_excitation(j,i) = banned_excitation(i,j)
+  enddo
+ enddo
+END_PROVIDER
+
+
 
 integer*8 function get_mo_map_size()
   implicit none
diff --git a/src/tools/molden.irp.f b/src/tools/molden.irp.f
index a7d5b978..417b25ad 100644
--- a/src/tools/molden.irp.f
+++ b/src/tools/molden.irp.f
@@ -17,7 +17,7 @@ program molden
 
   write(i_unit_output,'(A)') '[Molden Format]'
 
-  write(i_unit_output,'(A)') '[Atoms] ANGSTROM'
+  write(i_unit_output,'(A)') '[Atoms] Angs'
   do i = 1, nucl_num
     write(i_unit_output,'(A2,2X,I4,2X,I4,3(2X,F15.10))')                   &
         trim(element_name(int(nucl_charge(i)))),                     &
diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f
index a8dea97a..deb4f4ca 100644
--- a/src/utils/linear_algebra.irp.f
+++ b/src/utils/linear_algebra.irp.f
@@ -11,7 +11,7 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n)
 
   integer, intent(in)             :: LDA, LDU, LDVt, m, n
   double precision, intent(in)    :: A(LDA,n)
-  double precision, intent(out)   :: U(LDU,m)
+  double precision, intent(out)   :: U(LDU,min(m,n))
   double precision,intent(out)    :: Vt(LDVt,n)
   double precision,intent(out)    :: D(min(m,n))
   double precision,allocatable    :: work(:)
@@ -19,19 +19,19 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n)
 
   double precision,allocatable    :: A_tmp(:,:)
   allocate (A_tmp(LDA,n))
-  A_tmp = A
+  A_tmp(:,:) = A(:,:)
 
   ! Find optimal size for temp arrays
   allocate(work(1))
   lwork = -1
-  call dgesvd('A','A', m, n, A_tmp, LDA,                             &
+  call dgesvd('S','S', m, n, A_tmp, LDA,                             &
       D, U, LDU, Vt, LDVt, work, lwork, info)
   ! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
   lwork = max(int(work(1)), 5*MIN(M,N))
   deallocate(work)
 
   allocate(work(lwork))
-  call dgesvd('A','A', m, n, A_tmp, LDA,                             &
+  call dgesvd('S','S', m, n, A_tmp, LDA,                             &
       D, U, LDU, Vt, LDVt, work, lwork, info)
   deallocate(work,A_tmp)
 
@@ -42,6 +42,128 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n)
 
 end
 
+subroutine eigSVD(A,LDA,U,LDU,D,Vt,LDVt,m,n)
+  implicit none
+  BEGIN_DOC
+! Algorithm 3 of https://arxiv.org/pdf/1810.06860.pdf
+!
+! A(m,n) = U(m,n) D(n) Vt(n,n) with m>n
+  END_DOC
+  integer, intent(in)            :: LDA, LDU, LDVt, m, n
+  double precision, intent(in)   :: A(LDA,n)
+  double precision, intent(out)  :: U(LDU,n)
+  double precision,intent(out)   :: Vt(LDVt,n)
+  double precision,intent(out)   :: D(n)
+
+  integer                        :: i,j,k
+
+  if (m<n) then
+    stop -1
+    call svd(A,LDA,U,LDU,D,Vt,LDVt,m,n)
+    return
+  endif
+
+  double precision, allocatable :: B(:,:), V(:,:)
+  allocate(B(n,n))
+  ! B = - At . A
+  call dgemm('T','N',n,n,m,-1.d0,A,size(A,1),A,size(A,1),0.d0,B,size(B,1))
+
+  ! V, D = eig(B)
+  allocate(V(n,n))
+  call lapack_diagd(D,V,B,n,n)
+  deallocate(B)
+  do j=1,n
+   do i=1,n
+     Vt(i,j) = V(j,i)
+   enddo
+  enddo
+
+  ! S = sqrt(-D)
+  ! U = A.V.S^-1
+  ! U = A.(S^-1.vt)t
+
+  do k=1,n
+    if (D(k) >= 0.d0) then
+      exit
+    endif
+    D(k) = dsqrt(-D(k))
+    call dscal(n, 1.d0/D(k), V(1,k), 1)
+  enddo
+  D(k:n) = 0.d0
+  k=k-1
+  call dgemm('N','N',m,n,k,1.d0,A,size(A,1),V,size(V,1),0.d0,U,size(U,1))
+
+end
+
+
+subroutine randomized_svd(A,LDA,U,LDU,D,Vt,LDVt,m,n,q,r)
+  implicit none
+  include 'constants.include.F'
+  BEGIN_DOC
+! Randomized SVD: rank r, q power iterations
+!
+! 1. Sample column space of A with P: Z = A.P where P is random with r+p columns.
+!
+! 2. Power iterations : Z <- X . (Xt.Z)
+!
+! 3. Z = Q.R
+!
+! 4. Compute SVD on projected Qt.X = U' . S. Vt
+!
+! 5. U = Q U'
+  END_DOC
+
+  integer, intent(in)             :: LDA, LDU, LDVt, m, n, q, r
+  double precision, intent(in)    :: A(LDA,n)
+  double precision, intent(out)   :: U(LDU,r)
+  double precision,intent(out)    :: Vt(LDVt,r)
+  double precision,intent(out)    :: D(r)
+  integer                         :: i, j, k
+
+  double precision,allocatable    :: Z(:,:), P(:,:), Y(:,:), UY(:,:)
+  double precision :: r1,r2
+  allocate(P(n,r), Z(m,r))
+
+  ! P is a normal random matrix (n,r)
+  do k=1,r
+    do i=1,n
+      call random_number(r1)
+      call random_number(r2)
+      r1 = dsqrt(-2.d0*dlog(r1))
+      r2 = dtwo_pi*r2
+      P(i,k) = r1*dcos(r2)
+    enddo
+  enddo
+
+  ! Z(m,r) = A(m,n).P(n,r)
+  call dgemm('N','N',m,r,n,1.d0,A,size(A,1),P,size(P,1),0.d0,Z,size(Z,1))
+
+  ! Power iterations
+  do k=1,q
+    ! P(n,r) = At(n,m).Z(m,r)
+    call dgemm('T','N',n,r,m,1.d0,A,size(A,1),Z,size(Z,1),0.d0,P,size(P,1))
+    ! Z(m,r) = A(m,n).P(n,r)
+    call dgemm('N','N',m,r,n,1.d0,A,size(A,1),P,size(P,1),0.d0,Z,size(Z,1))
+  enddo
+
+  deallocate(P)
+
+  ! QR factorization of Z
+  call ortho_svd(Z,size(Z,1),m,r)
+
+  allocate(Y(r,n), UY(r,r))
+  ! Y(r,n) = Zt(r,m).A(m,n)
+  call dgemm('T','N',r,n,m,1.d0,Z,size(Z,1),A,size(A,1),0.d0,Y,size(Y,1))
+
+  ! SVD of Y
+  call svd(Y,size(Y,1),UY,size(UY,1),D,Vt,size(Vt,1),r,n)
+  deallocate(Y)
+
+  ! U(m,r) = Z(m,r).UY(r,r)
+  call dgemm('N','N',m,r,r,1.d0,Z,size(Z,1),UY,size(UY,1),0.d0,U,size(U,1))
+  deallocate(UY,Z)
+
+end
 
 subroutine svd_complex(A,LDA,U,LDU,D,Vt,LDVt,m,n)
   implicit none
@@ -807,6 +929,33 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m,cutoff)
 end
 
 
+subroutine ortho_svd(A,LDA,m,n)
+  implicit none
+  BEGIN_DOC
+  ! Orthogonalization via fast SVD
+  !
+  ! A : matrix to orthogonalize
+  !
+  ! LDA : leftmost dimension of A
+  !
+  ! m : Number of rows of A
+  !
+  ! n : Number of columns of A
+  !
+  END_DOC
+  integer, intent(in)            :: m,n, LDA
+  double precision, intent(inout) :: A(LDA,n)
+  if (m < n) then
+    call ortho_qr(A,LDA,m,n)
+  endif
+  double precision, allocatable :: U(:,:), D(:), Vt(:,:)
+  allocate(U(m,n), D(n), Vt(n,n))
+  call SVD(A,LDA,U,size(U,1),D,Vt,size(Vt,1),m,n)
+  A(1:m,1:n) = U(1:m,1:n)
+  deallocate(U,D, Vt)
+
+end
+
 subroutine ortho_qr(A,LDA,m,n)
   implicit none
   BEGIN_DOC