From 770b4f6628f920848fbd2b599c53ff11a41456b2 Mon Sep 17 00:00:00 2001
From: Kevin Gasperich <kgasperich@users.noreply.github.com>
Date: Fri, 21 Aug 2020 12:30:14 -0500
Subject: [PATCH] compute s_z2_sz in diag_S_mat_elem

---
 src/determinants/s2.irp.f | 21 ++++++++++++++++-----
 1 file changed, 16 insertions(+), 5 deletions(-)

diff --git a/src/determinants/s2.irp.f b/src/determinants/s2.irp.f
index 391d0073..c1dd06c0 100644
--- a/src/determinants/s2.irp.f
+++ b/src/determinants/s2.irp.f
@@ -8,24 +8,35 @@ double precision function diag_S_mat_elem(key_i,Nint)
   BEGIN_DOC
 ! Returns <i|S^2|i>
   END_DOC
-  integer                        :: nup, i
-  integer(bit_kind)              :: xorvec(N_int_max)
+  integer                        :: nup, ntot, i
+  integer(bit_kind)              :: xorvec(N_int_max), upvec(N_int_max)
 
   do i=1,Nint
     xorvec(i) = xor(key_i(i,1),key_i(i,2))
   enddo
 
   do i=1,Nint
-    xorvec(i) = iand(xorvec(i),key_i(i,1))
+    upvec(i) = iand(xorvec(i),key_i(i,1))
   enddo
 
+  ! nup is number of alpha unpaired
+  ! ntot is total number of unpaired
   nup = 0
+  ntot = 0
   do i=1,Nint
     if (xorvec(i) /= 0_bit_kind) then
-      nup += popcnt(xorvec(i))
+      ntot += popcnt(xorvec(i))
+      if (upvec(i) /= 0_bit_kind) then
+        nup += popcnt(upvec(i))
+      endif
     endif
   enddo
-  diag_S_mat_elem = dble(nup)
+  
+  double precision :: sz
+  sz = nup - 0.5d0*ntot
+  
+  !<S^2> = <S+ S-> + Sz(Sz-1)
+  diag_S_mat_elem = nup + sz*(sz-1)
 
 end