From dc0d668f3829fee928b645f469a6c43f90608a8d Mon Sep 17 00:00:00 2001
From: amandadumi <amanda.e.dumi@gmail.com>
Date: Tue, 16 Jun 2020 10:57:24 -0400
Subject: [PATCH 01/17] changing molden 'Atoms' label to match coordinate units

---
 src/tools/molden.irp.f | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/tools/molden.irp.f b/src/tools/molden.irp.f
index f70ed0de..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] AU'
+  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)))),                     &

From d65c2cba57a1b794542432cef66d7a0619e17e90 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed, 2 Sep 2020 15:53:56 +0200
Subject: [PATCH 02/17] Fixed NaN bug in DIIS

---
 src/scf_utils/roothaan_hall_scf.irp.f | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f
index faf23a51..91c85f5e 100644
--- a/src/scf_utils/roothaan_hall_scf.irp.f
+++ b/src/scf_utils/roothaan_hall_scf.irp.f
@@ -208,8 +208,8 @@ END_DOC
   do j=1,dim_DIIS
     do i=1,dim_DIIS
 
-      j_DIIS = mod(iteration_SCF-j,max_dim_DIIS)+1
-      i_DIIS = mod(iteration_SCF-i,max_dim_DIIS)+1
+      j_DIIS = min(mod(iteration_SCF-j,max_dim_DIIS)+1,dim_DIIS)
+      i_DIIS = min(mod(iteration_SCF-i,max_dim_DIIS)+1,dim_DIIS)
 
 ! Compute product of two errors vectors
 

From 42b74b743f07739c31d84d5f2fe174ccfcb34d2f Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 3 Sep 2020 09:27:24 +0200
Subject: [PATCH 03/17] Stupid typo

---
 src/tools/molden.irp.f | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/tools/molden.irp.f b/src/tools/molden.irp.f
index 239baf18..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] Angs
+  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)))),                     &

From c07232aae3ccbb640d29e444b70c558c287261f6 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 3 Sep 2020 11:09:25 +0200
Subject: [PATCH 04/17] Introduce <PT2_i|PT2_j> in PT2

---
 src/cipsi/selection.irp.f | 35 ++++++++++++++++++++++-------------
 1 file changed, 22 insertions(+), 13 deletions(-)

diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f
index 071d9294..22d737a1 100644
--- a/src/cipsi/selection.irp.f
+++ b/src/cipsi/selection.irp.f
@@ -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), X(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
 
@@ -772,21 +773,26 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
         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
+        if (e_pert(istate) < 0.d0) then
+           X(istate) = -dsqrt(-e_pert(istate))
+        else
+           X(istate) = dsqrt(e_pert(istate))
+        endif
       enddo
 
-      ! Gram-Schmidt using input overlap matrix
-      do istate=1,N_states
-        do jstate=1,istate-1
-          if ( (pt2_overlap(jstate,istate) == 0.d0).or.(pt2_overlap(jstate,jstate) == 0.d0) ) cycle
-          coef(istate) = coef(istate) - pt2_overlap(jstate,istate)/pt2_overlap(jstate,jstate) * coef(jstate)
-        enddo
-      enddo
+!      ! Gram-Schmidt using input overlap matrix
+!      do istate=1,N_states
+!        do jstate=1,istate-1
+!          if ( (pt2_overlap(jstate,istate) == 0.d0).or.(pt2_overlap(jstate,jstate) == 0.d0) ) cycle
+!          coef(istate) = coef(istate) - pt2_overlap(jstate,istate)/pt2_overlap(jstate,jstate) * coef(jstate)
+!        enddo
+!      enddo
 
       do istate=1, N_states
         do jstate=1,N_states
@@ -796,10 +802,9 @@ 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
 
         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
@@ -831,7 +836,11 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
 
           case default
             ! Energy selection
-            w = w + e_pert * selection_weight(istate)
+            w = w + e_pert(istate) * selection_weight(istate)
+            do jstate=1,N_states
+              if (istate == jstate) cycle
+              w = w - dabs(X(istate))*X(jstate) * sqrt(selection_weight(istate)*selection_weight(jstate))
+            enddo
 
         end select
       end do

From fd7825ea742163b80950c8495996e4b871c95c83 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 3 Sep 2020 16:16:16 +0200
Subject: [PATCH 05/17] Changed default weight_selection to 1

---
 src/cipsi/pt2_stoch_routines.irp.f | 15 +++++++++++----
 src/cipsi/selection.irp.f          | 23 ++++++++++++++++++-----
 src/determinants/EZFIO.cfg         |  2 +-
 3 files changed, 30 insertions(+), 10 deletions(-)

diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f
index f0bd8847..faced03a 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
 
@@ -529,7 +529,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, F11.7, X, G10.3, X, G10.3, X, G10.3, X, G10.3, 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 22d737a1..e599737c 100644
--- a/src/cipsi/selection.irp.f
+++ b/src/cipsi/selection.irp.f
@@ -683,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
@@ -829,17 +834,25 @@ 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)
+            do jstate=1,N_states
+              if (istate == jstate) cycle
+              w = w + 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)
+            do jstate=1,N_states
+              if (istate == jstate) cycle
+              w = w + coef(istate)*coef(jstate) * s_weight(istate,jstate)
+            enddo
 
           case default
             ! Energy selection
-            w = w + e_pert(istate) * selection_weight(istate)
+            w = w + e_pert(istate) * s_weight(istate,istate)
             do jstate=1,N_states
               if (istate == jstate) cycle
-              w = w - dabs(X(istate))*X(jstate) * sqrt(selection_weight(istate)*selection_weight(jstate))
+              w = w - dabs(X(istate))*X(jstate) * s_weight(istate,jstate)
             enddo
 
         end select
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

From 8c75ad6cfab211da3519d9d7aad98676ac0a2b7f Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 3 Sep 2020 17:08:45 +0200
Subject: [PATCH 06/17] Prints

---
 src/cipsi/pt2_stoch_routines.irp.f | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f
index faced03a..72f00d70 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)', '========== ====================== ===================== ===================== ==========='
       print '(A)', ' Samples          Energy                Variance               Norm^2          Seconds'
-      print '(A)', '========== ======================= ==================== ===================== =========='
+      print '(A)', '========== ====================== ===================== ===================== ==========='
 
       PROVIDE global_selection_buffer
 
@@ -312,7 +312,7 @@ 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)

From 07785a6db19a93330d212c5e692797d7a9de2347 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 3 Sep 2020 17:13:26 +0200
Subject: [PATCH 07/17] Print overlap in deterministic PT2

---
 src/cipsi/pt2_stoch_routines.irp.f | 11 ++---------
 src/cipsi/zmq_selection.irp.f      |  6 ++++++
 2 files changed, 8 insertions(+), 9 deletions(-)

diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f
index 72f00d70..4d6ea1e8 100644
--- a/src/cipsi/pt2_stoch_routines.irp.f
+++ b/src/cipsi/pt2_stoch_routines.irp.f
@@ -317,14 +317,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
     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 *, 'Overlap of perturbed states:'
 print *, pt2_overlap(pt2_stoch_istate,:)
 print *, '-------'
     SOFT_TOUCH pt2_overlap
@@ -529,7 +522,7 @@ 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 '(I10, X, F11.7, X, G10.3, X, G10.3, X, G10.3, X, G10.3, X, G10.3, X, F10.4)', c, &
+            print '(I10, X, F11.7, X, G10.4, X, G10.3, X, G10.3, X, G10.3, 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), &
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)
 

From 5fcdbe12dfa1ba81202ab001a2b222c2b22fe055 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 3 Sep 2020 17:13:26 +0200
Subject: [PATCH 08/17] Print overlap in deterministic PT2

---
 src/cipsi/pt2_stoch_routines.irp.f | 11 ++---------
 src/cipsi/zmq_selection.irp.f      |  6 ++++++
 2 files changed, 8 insertions(+), 9 deletions(-)

diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f
index 72f00d70..7622e93f 100644
--- a/src/cipsi/pt2_stoch_routines.irp.f
+++ b/src/cipsi/pt2_stoch_routines.irp.f
@@ -317,14 +317,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
     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 *, 'Overlap of perturbed states:'
 print *, pt2_overlap(pt2_stoch_istate,:)
 print *, '-------'
     SOFT_TOUCH pt2_overlap
@@ -529,7 +522,7 @@ 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 '(I10, X, F11.7, X, G10.3, X, G10.3, X, G10.3, X, G10.3, X, G10.3, X, F10.4)', c, &
+            print '(I10, X, F11.7, X, F10.6, X, G10.3, X, F10.6, X, G10.3, 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), &
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)
 

From 0a5f3ac330cb9d75c641fc32ecf3b2c319558040 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 3 Sep 2020 18:12:58 +0200
Subject: [PATCH 09/17] PT2 overlap

---
 src/cipsi/pt2_stoch_routines.irp.f | 17 ++++++++++++++---
 1 file changed, 14 insertions(+), 3 deletions(-)

diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f
index 7622e93f..ab7e765a 100644
--- a/src/cipsi/pt2_stoch_routines.irp.f
+++ b/src/cipsi/pt2_stoch_routines.irp.f
@@ -317,14 +317,25 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
     do k=1,N_states
       pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate)
     enddo
-print *, 'Overlap of perturbed states:'
-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

From f0454b5b344f011520e6f2c70ded66343fa9ded0 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 3 Sep 2020 23:14:08 +0200
Subject: [PATCH 10/17] Print format of PT2

---
 src/cipsi/pt2_stoch_routines.irp.f | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f
index 1d13c737..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)', '========== ======================= ===================== ===================== ==========='
       print '(A)', ' Samples          Energy                Variance               Norm^2          Seconds'
-      print '(A)', '========== ====================== ===================== ===================== ==========='
+      print '(A)', '========== ======================= ===================== ===================== ==========='
 
       PROVIDE global_selection_buffer
 
@@ -312,7 +312,7 @@ 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)
@@ -533,7 +533,7 @@ 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 '(I10, X, F11.7, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.4)', c, &
+            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), &

From 82e4299ad6e72742a7388f4e3a889c663c7a67b0 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed, 23 Sep 2020 18:58:07 +0200
Subject: [PATCH 11/17] Introduced nxn diagonalization for PT2

---
 src/cipsi/selection.irp.f      | 89 ++++++++++++++++++++++++++--------
 src/utils/linear_algebra.irp.f |  2 +-
 2 files changed, 70 insertions(+), 21 deletions(-)

diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f
index e599737c..024dd435 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,6 +143,23 @@ BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ]
 
 END_PROVIDER
 
+BEGIN_PROVIDER [ logical, banned_excitation, (mo_num,mo_num) ]
+ implicit none
+ BEGIN_DOC
+ ! If true, the excitation is banned in the selection. Useful with local MOs.
+ END_DOC
+ banned_excitation = .False.
+ integer :: i,j
+ double precision :: buffer(mo_num)
+ do j=1,mo_num
+   call get_mo_two_e_integrals_exch_ii(j,j,mo_num,buffer,mo_integrals_map)
+   buffer = dabs(buffer)
+   do i=1,mo_num
+    banned_excitation(i,j) = buffer(i) < 1.d-15
+  enddo
+ enddo
+END_PROVIDER
+
 
 subroutine get_mask_phase(det1, pm, Nint)
   use bitmasks
@@ -288,6 +304,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 +624,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 +692,7 @@ 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(N_states), coef(N_states), X(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
@@ -769,10 +787,16 @@ 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
@@ -784,13 +808,38 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
         else
           coef(istate) = alpha_h_psi / delta_E
         endif
-        if (e_pert(istate) < 0.d0) then
-           X(istate) = -dsqrt(-e_pert(istate))
-        else
-           X(istate) = dsqrt(e_pert(istate))
-        endif
       enddo
 
+      do_diag = maxval(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 = maxloc(pt2_matrix,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
@@ -835,25 +884,25 @@ 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 * s_weight(istate,istate)
-            do jstate=1,N_states
-              if (istate == jstate) cycle
-              w = w + alpha_h_psi*mat(jstate,p1,p2) * s_weight(istate,jstate)
-            enddo
+!            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) * s_weight(istate,istate)
-            do jstate=1,N_states
-              if (istate == jstate) cycle
-              w = w + coef(istate)*coef(jstate) * s_weight(istate,jstate)
-            enddo
+!            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(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
+!            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/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f
index a8dea97a..cf43edd5 100644
--- a/src/utils/linear_algebra.irp.f
+++ b/src/utils/linear_algebra.irp.f
@@ -19,7 +19,7 @@ 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))

From 161517e0ea4d520fd4ee9a5db7ba3aa9fb3d8d61 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed, 23 Sep 2020 19:21:19 +0200
Subject: [PATCH 12/17] Using max instead of sum in pt2

---
 src/cipsi/selection.irp.f | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f
index 024dd435..27e47f03 100644
--- a/src/cipsi/selection.irp.f
+++ b/src/cipsi/selection.irp.f
@@ -810,7 +810,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
         endif
       enddo
 
-      do_diag = maxval(dabs(coef)) > 0.001d0
+      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)
@@ -883,14 +883,16 @@ 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 * s_weight(istate,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) * s_weight(istate,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)
@@ -898,7 +900,8 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
 
           case default
             ! Energy selection
-            w = w + e_pert(istate) * s_weight(istate,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)

From e8388681812b76f2a302094a3283f70ab6b5722e Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Fri, 25 Sep 2020 15:14:56 +0200
Subject: [PATCH 13/17] Fixed Schwartz and added banned-excitations

---
 ocaml/Units.ml                          |  4 +--
 src/ao_two_e_ints/two_e_integrals.irp.f |  2 +-
 src/cipsi/selection.irp.f               | 18 --------------
 src/mo_two_e_ints/map_integrals.irp.f   | 33 +++++++++++++++++++++++++
 4 files changed, 36 insertions(+), 21 deletions(-)

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/selection.irp.f b/src/cipsi/selection.irp.f
index 27e47f03..f1eb3425 100644
--- a/src/cipsi/selection.irp.f
+++ b/src/cipsi/selection.irp.f
@@ -143,24 +143,6 @@ BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ]
 
 END_PROVIDER
 
-BEGIN_PROVIDER [ logical, banned_excitation, (mo_num,mo_num) ]
- implicit none
- BEGIN_DOC
- ! If true, the excitation is banned in the selection. Useful with local MOs.
- END_DOC
- banned_excitation = .False.
- integer :: i,j
- double precision :: buffer(mo_num)
- do j=1,mo_num
-   call get_mo_two_e_integrals_exch_ii(j,j,mo_num,buffer,mo_integrals_map)
-   buffer = dabs(buffer)
-   do i=1,mo_num
-    banned_excitation(i,j) = buffer(i) < 1.d-15
-  enddo
- enddo
-END_PROVIDER
-
-
 subroutine get_mask_phase(det1, pm, Nint)
   use bitmasks
   implicit none
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

From 88b798fe22d0c3d5c718686d11772a281c31308c Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed, 30 Sep 2020 20:48:27 +0200
Subject: [PATCH 14/17] Added randomized SVD routines

---
 src/utils/linear_algebra.irp.f | 155 ++++++++++++++++++++++++++++++++-
 1 file changed, 152 insertions(+), 3 deletions(-)

diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f
index cf43edd5..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(:)
@@ -24,14 +24,14 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n)
   ! 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

From 28cb4a516702fe0da6b2a7e63b0f1fe8d5cb9b10 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Fri, 23 Oct 2020 11:35:45 +0200
Subject: [PATCH 15/17] Change default in davidson

---
 src/davidson/EZFIO.cfg | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

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

From b1484089386b4321527a4456ff0fc1fa6d737f35 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Fri, 23 Oct 2020 16:58:07 +0200
Subject: [PATCH 16/17] Fixed maxloc

---
 src/cipsi/selection.irp.f | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f
index f1eb3425..00ac915e 100644
--- a/src/cipsi/selection.irp.f
+++ b/src/cipsi/selection.irp.f
@@ -815,7 +815,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
           stop -1
         endif
         pt2_matrix = dabs(pt2_matrix)
-        iwork = maxloc(pt2_matrix,1)
+        iwork(1:N_states) = maxloc(pt2_matrix,1)
         do k=1,N_states
           e_pert(iwork(k)) = eigvalues(k) - E0(iwork(k))
         enddo

From 862a1804a4e7e3081771f8e0b9007ba750712a6f Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Sat, 24 Oct 2020 14:11:04 +0200
Subject: [PATCH 17/17] Fixed Maxloc

---
 src/cipsi/selection.irp.f | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f
index 00ac915e..08595aea 100644
--- a/src/cipsi/selection.irp.f
+++ b/src/cipsi/selection.irp.f
@@ -815,7 +815,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
           stop -1
         endif
         pt2_matrix = dabs(pt2_matrix)
-        iwork(1:N_states) = maxloc(pt2_matrix,1)
+        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