From 243f46cbcae25e0b2b8aab18db629a0cf4740e78 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue, 13 Sep 2016 15:36:36 +0200
Subject: [PATCH] Don't mix same spin symmetry with s2eig

---
 src/Determinants/davidson.irp.f | 8 ++++----
 src/Determinants/s2.irp.f       | 6 ++++--
 2 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/src/Determinants/davidson.irp.f b/src/Determinants/davidson.irp.f
index 5bc3768e..dbf437a1 100644
--- a/src/Determinants/davidson.irp.f
+++ b/src/Determinants/davidson.irp.f
@@ -457,7 +457,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
       ! ----------------------
       
       do k=1,N_st
-          call H_u_0(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint)
+        call H_u_0(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint)
       enddo
       
       
@@ -530,7 +530,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
       !$OMP END PARALLEL
       
       write(iunit,'(X,I3,X,100(X,F16.10,X,E16.6))')  iter, to_print(:,1:N_st)
-      call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_states,converged)
+      call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_states_diag,converged)
       if (converged) then
         exit
       endif
@@ -554,14 +554,14 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
           do l=1,N_st
             c = u_dot_v(U(1,k,iter+1),U(1,l,iter2),sze)
             do i=1,sze
-              U(i,k,iter+1) -= c * U(i,l,iter2)
+              U(i,k,iter+1) = U(i,k,iter+1) - c * U(i,l,iter2)
             enddo
           enddo
         enddo
         do l=1,k-1
           c = u_dot_v(U(1,k,iter+1),U(1,l,iter+1),sze)
           do i=1,sze
-            U(i,k,iter+1) -= c * U(i,l,iter+1)
+            U(i,k,iter+1) = U(i,k,iter+1) - c * U(i,l,iter+1)
           enddo
         enddo
         call normalize( U(1,k,iter+1), sze )
diff --git a/src/Determinants/s2.irp.f b/src/Determinants/s2.irp.f
index 96e342cd..d34aad88 100644
--- a/src/Determinants/s2.irp.f
+++ b/src/Determinants/s2.irp.f
@@ -348,7 +348,8 @@ subroutine diagonalize_s2_betweenstates(keys_tmp,psi_coefs_inout,n,nmax_keys,nma
  accu_precision_of_diag = 0.d0
  do i = 1, nstates
   do j = i+1, nstates
-   if(  ( dabs(s2(i,i) - s2(j,j)) .le.1.d-10 ) .and. (dabs(s2(i,j) + dabs(s2(i,j)))) .le.1.d-10) then
+   ! Do not combine states of the same spin symmetry
+   if( dabs(s2(i,i) - s2(j,j)) .le.0.5d0) then
     s2(i,j) = 0.d0
     s2(j,i) = 0.d0
    endif
@@ -356,13 +357,14 @@ subroutine diagonalize_s2_betweenstates(keys_tmp,psi_coefs_inout,n,nmax_keys,nma
  enddo
  do i = 1, nstates
   write(*,'(10(F10.6,X))')s2(i,:)
+  s2(i,i) = s2(i,i) 
  enddo
 
  print*,'Diagonalizing the S^2 matrix'
 
  allocate(eigvalues(nstates),eigvectors(nstates,nstates))
  call lapack_diagd(eigvalues,eigvectors,s2,nstates,nstates)
- print*,'Eigenvalues of s^2'
+ print*,'Shifted Eigenvalues of s^2'
  do i = 1, nstates
   print*,'s2 = ',eigvalues(i)
   s2_eigvalues(i) = eigvalues(i)