From 3a71cf0dc6aba0e914aafd61445dc64683c50b59 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed, 16 Jun 2021 09:49:53 +0200
Subject: [PATCH] CSF-based davidson as an option0

---
 src/davidson/EZFIO.cfg            | 6 ++++++
 src/davidson/diagonalize_ci.irp.f | 4 +---
 src/davidson/u0_h_u0.irp.f        | 2 +-
 src/davidson/u0_hs2_u0.irp.f      | 4 ++--
 4 files changed, 10 insertions(+), 6 deletions(-)

diff --git a/src/davidson/EZFIO.cfg b/src/davidson/EZFIO.cfg
index 4343cc1e..a3efc2be 100644
--- a/src/davidson/EZFIO.cfg
+++ b/src/davidson/EZFIO.cfg
@@ -34,6 +34,12 @@ doc: If |true|, a memory-mapped file may be used to store the W and S2 vectors i
 default: True
 interface: ezfio,provider,ocaml
 
+[csf_based]
+type: logical
+doc: If |true|, use the CSF-based algorithm
+default: False
+interface: ezfio,provider,ocaml
+
 [distributed_davidson]
 type: logical
 doc: If |true|, use the distributed algorithm
diff --git a/src/davidson/diagonalize_ci.irp.f b/src/davidson/diagonalize_ci.irp.f
index da5fb950..be4250bf 100644
--- a/src/davidson/diagonalize_ci.irp.f
+++ b/src/davidson/diagonalize_ci.irp.f
@@ -56,9 +56,7 @@ END_PROVIDER
      enddo
    enddo
 
-!   Deactivated temporarily: bug in N_csf
-!   do_csf = s2_eig .and. only_expected_s2 .and. (expected_s2 == 0.d0)
-   do_csf = .False.
+   do_csf = s2_eig .and. only_expected_s2 .and. csf_based
 
    if (diag_algorithm == "Davidson") then
 
diff --git a/src/davidson/u0_h_u0.irp.f b/src/davidson/u0_h_u0.irp.f
index 3f5113db..7ef154a3 100644
--- a/src/davidson/u0_h_u0.irp.f
+++ b/src/davidson/u0_h_u0.irp.f
@@ -46,7 +46,7 @@ subroutine u_0_H_u_0(e_0,u_0,n,keys_tmp,Nint,N_st,sze)
   do i=1,N_st
     norm = u_dot_u(u_0(1,i),n)
     if (norm /= 0.d0) then
-      e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)
+      e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n) / dsqrt(norm)
     else
       e_0(i) = 0.d0
     endif
diff --git a/src/davidson/u0_hs2_u0.irp.f b/src/davidson/u0_hs2_u0.irp.f
index 96c266c2..8f7bf06b 100644
--- a/src/davidson/u0_hs2_u0.irp.f
+++ b/src/davidson/u0_hs2_u0.irp.f
@@ -75,8 +75,8 @@ subroutine u_0_HS2_u_0(e_0,s_0,u_0,n,keys_tmp,Nint,N_st,sze)
   do i=1,N_st
     norm = u_dot_u(u_0(1,i),n)
     if (norm /= 0.d0) then
-      e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)
-      s_0(i) = u_dot_v(s_vec(1,i),u_0(1,i),n)
+      e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)/norm
+      s_0(i) = u_dot_v(s_vec(1,i),u_0(1,i),n)/norm
     else
       e_0(i) = 0.d0
       s_0(i) = 0.d0