From 450672c82592243d21192aaef2275f87af1e126b Mon Sep 17 00:00:00 2001
From: Pierre-Francois Loos <loos@lcpqmc74.univ-tlse3.fr>
Date: Thu, 16 Nov 2023 13:32:31 +0100
Subject: [PATCH] saving work

---
 src/HF/print_GHF.f90 | 49 ++++++++++++++++++++++++++++++++++----------
 1 file changed, 38 insertions(+), 11 deletions(-)

diff --git a/src/HF/print_GHF.f90 b/src/HF/print_GHF.f90
index b140b42..e83bac0 100644
--- a/src/HF/print_GHF.f90
+++ b/src/HF/print_GHF.f90
@@ -32,7 +32,7 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
   double precision                   :: Gap
   double precision                   :: Sx ,Sy ,Sz
   double precision                   :: Sx2,Sy2,Sz2
-  double precision                   :: S2
+  double precision                   :: SmSp,SpSm,S2
 
   double precision,allocatable       :: Ca(:,:)
   double precision,allocatable       :: Cb(:,:)
@@ -82,27 +82,54 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
 
 ! Compute <S^2> = <Sx^2> + <Sy^2> + <Sz^2>
 
-  Sx2 = 0.25d0*trace_matrix(nO,Paa+Pbb) + 0.25d0*trace_matrix(nO,Pab+Pba)**2 &
-      - 0.5d0*trace_matrix(nO,matmul(Paa,Pbb) + matmul(Pab,Pab))
+  SpSm = 0d0
+  do i=1,nO
+    do j=1,nO
+      SpSm = SpSm + Pab(i,i)*Pba(j,j) - Pab(i,j)*Pba(j,i)
+    end do
+  end do
+  SpSm = 0.5d0*(trace_matrix(nO,Paa) + SpSm)
 
-  Sx2 = trace_matrix(
+! Sx2 = 0.25d0*trace_matrix(nO,Paa+Pbb) + 0.25d0*trace_matrix(nO,Pab+Pba)**2 &
+!     - 0.5d0*trace_matrix(nO,matmul(Paa,Pbb) + matmul(Pab,Pab))
 
-  Sy2 = 0.25d0*trace_matrix(nO,Paa+Pbb) - 0.25d0*trace_matrix(nO,Pab-Pba)**2 &
-      - 0.5d0*trace_matrix(nO,matmul(Paa,Pbb) - matmul(Pab,Pab))
+  SmSp = 0d0
+  do i=1,nO
+    do j=1,nO
+      SmSp = SmSp + Pba(i,i)*Pab(j,j) - Pba(i,j)*Pab(j,i)
+    end do
+  end do
+  SmSp = 0.5d0*(trace_matrix(nO,Pbb) + SmSp)
 
-  Sz2 = 0.25d0*trace_matrix(nO,Paa+Pbb) + 0.25d0*trace_matrix(nO,Paa-Pbb)**2 &
-      - 0.25d0*trace_matrix(nO,matmul(Paa,Paa) + matmul(Pbb,Pbb)) &
-      + 0.25d0*trace_matrix(nO,matmul(Pab,Pba) + matmul(Pba,Pab))
+! Sy2 = 0.25d0*trace_matrix(nO,Paa+Pbb) - 0.25d0*trace_matrix(nO,Pab-Pba)**2 &
+!     - 0.5d0*trace_matrix(nO,matmul(Paa,Pbb) - matmul(Pab,Pab))
+
+  Sz2 = 0d0
+  do i=1,nO
+    do j=1,nO
+      Sz2 = Sz2 + (Paa(i,i) - Pbb(i,i))*(Paa(j,j) - Pbb(j,j)) - (Paa(i,j) - Pbb(i,j))**2
+    end do
+  end do
+  Sz2 = 0.25d0*(dble(nO) + Sz2)
+  print*,'<Sz^2> = ',Sz2
+
+  Sz2 = 0.25d0*trace_matrix(nO,Paa+Pbb)    &
+      - 0.25d0*trace_matrix(nO,Paa-Pbb)**2 &
+      - 0.25d0*trace_matrix(nO,matmul(Paa-Pbb,Paa-Pbb))
 
   S2 = Sz*(Sz+1d0) + trace_matrix(nO,Pbb) + 0.25d0*trace_matrix(nO,Paa+Pbb)
-
   do i=1,nO
     do j=1,nO
       S2 = S2 - 0.25d0*(Paa(i,j) - Pbb(i,j))**2 &
               + (Pba(i,i)*Pab(j,j) - Pba(i,j)*Pab(j,i))
     end do   
   end do   
-! print*,'<S^2> = ',S2
+  print*,'<S^2> = ',S2
+
+  Sx2 = 0.5d0*(S2 - Sz2 + 0.5d0*(SmSp + SpSm))
+  print*,'<Sx^2> = ',Sx2
+  Sy2 = 0.5d0*(S2 - Sz2 - 0.5d0*(SmSp + SpSm))
+  print*,'<Sy^2> = ',Sy2
 
 ! deallocate(Paa,Pab,Pba,Pbb)