diff --git a/src/HF/print_GHF.f90 b/src/HF/print_GHF.f90
index 34eb238..4bfeb1a 100644
--- a/src/HF/print_GHF.f90
+++ b/src/HF/print_GHF.f90
@@ -41,10 +41,10 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
 
   double precision,allocatable       :: Ca(:,:)
   double precision,allocatable       :: Cb(:,:)
-  double precision,allocatable       :: Paa(:,:), Saa(:,:)
-  double precision,allocatable       :: Pab(:,:), Sab(:,:)
-  double precision,allocatable       :: Pba(:,:), Sba(:,:)
-  double precision,allocatable       :: Pbb(:,:), Sbb(:,:)
+  double precision,allocatable       :: Paa(:,:)
+  double precision,allocatable       :: Pab(:,:)
+  double precision,allocatable       :: Pba(:,:)
+  double precision,allocatable       :: Pbb(:,:)
   double precision,allocatable       :: tmp(:,:)
 
   double precision,allocatable       :: Mx(:,:)
@@ -80,6 +80,7 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
   Pba = matmul(transpose(Cb),matmul(S,Ca))
   Pbb = matmul(transpose(Cb),matmul(S,Cb))
 
+
 ! Compute components of S = (Sx,Sy,Sz)
 
   Sx = 0.5d0*(trace_matrix(nO,Pab) + trace_matrix(nO,Pba))
@@ -88,6 +89,7 @@ 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>
 
+
   SpSm = 0d0
   do i=1,nO
     do j=1,nO
@@ -124,6 +126,7 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
       - 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 &
@@ -139,12 +142,9 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
 !  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))
 
-!  Sx2 = trace_matrix(
-
 !  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 = 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))
@@ -159,77 +159,33 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
 !  end do   
 ! print*,'<S^2> = ',S2
 
-  ! TODO
-  ! check C size
-  allocate(Ca(nBas,nBas), Cb(nBas,nBas))
-  do i = 1, nBas
-    do j = 1, nBas
-      Ca(j,i) = C(j,     i)
-      Cb(j,i) = C(j,nBas+i)
-    enddo
-  enddo
 
-! allocate(Saa(nBas,nBas),Sab(nBas,nBas),Sba(nBas,nBas),Sbb(nBas,nBas))
-! allocate(tmp(nBas,nBas))
 
-  ! Saa = Ca x Sao x Ca.T
-! call dgemm("N", "N", nBas, nBas, nBas, 1.d0,  Ca,  size(Ca, 1), Sao, size(Sao, 1), 0.d0, tmp, size(tmp, 1))
-! call dgemm("N", "T", nBas, nBas, nBas, 1.d0, tmp, size(tmp, 1),  Ca,  size(Ca, 1), 0.d0, Saa, size(Saa, 1))
-
-  ! Sab = Ca x Sao x Cb.T
-! call dgemm("N", "N", nBas, nBas, nBas, 1.d0,  Ca,  size(Ca, 1), Sao, size(Sao, 1), 0.d0, tmp, size(tmp, 1))
-! call dgemm("N", "T", nBas, nBas, nBas, 1.d0, tmp, size(tmp, 1),  Cb,  size(Cb, 1), 0.d0, Sab, size(Sab, 1))
-
-  ! Sba = Cb x Sao x Ca.T
-  !     = Sab.T
-! Sba = transpose(Sab)
-
-  ! Sbb = Cb x Sao x Cb.T
-! call dgemm("N", "N", nBas, nBas, nBas, 1.d0,  Cb,  size(Cb, 1), Sao, size(Sao, 1), 0.d0, tmp, size(tmp, 1))
-! call dgemm("N", "T", nBas, nBas, nBas, 1.d0, tmp, size(tmp, 1),  Cb,  size(Cb, 1), 0.d0, Sbb, size(Sbb, 1))
-
-! deallocate(tmp)
-
-  ! TODO
-  ! nO = nb of electrons ?
-! na = 0.d0
-! nb = 0.d0
-! do i = 1, nO
-!   na = na + Saa(i,i)
-!   nb = nb + Sbb(i,i)
-! enddo
-  
-! nonco_z = dble(nO)
-! do j = 1, nO
-!   do i = 1, nO
-!     nonco_z = nonco_z - (Saa(i,j) - Sbb(i,j))**2
-!   enddo
-! enddo
-! nonco_z = 0.25d0 * nonco_z
-
-! Sz  = 0.5d0 * (na - nb)
-! Sz2 = Sz*Sz + nonco_z
-
-  ! If Na > Nb
-  !contam_uhf = nb
-  !do j = 1, nO
-  !  do i = 1, nO
-  !    contam_uhf = contam_uhf - (Sab(i,j) - Sba(j,i))
-  !  enddo
-  !enddo
-  !xy_perp = 0.d0
-  !do i = 1, nO
-  !  xy_perp = xy_perp + (Sba(i,i))**2
-  !enddo
-  !S2 = Sz * (Sz + 1.d0) + nonco_z + contam_uhf + xy_perp
-
-! contam_ghf = 0.d0
-! do j = 1, nO
-!   do i = 1, nO
-!     contam_ghf = contam_ghf - (Sab(i,i)*Sba(j,j) - Sab(i,j)*Sba(j,i))
-!   enddo
-! enddo
-! S2 = Sz * (Sz + 1.d0) + nonco_z + contam_ghf
+!  na = 0.d0
+!  nb = 0.d0
+!  do i = 1, nO
+!    na = na + Paa(i,i)
+!    nb = nb + Pbb(i,i)
+!  enddo
+!  Sz  = 0.5d0 * (na - nb)
+!  
+!  nonco_z = dble(nO)
+!  do j = 1, nO
+!    do i = 1, nO
+!      nonco_z = nonco_z - (Paa(i,j) - Pbb(i,j))**2
+!    enddo
+!  enddo
+!  nonco_z = 0.25d0 * nonco_z
+!
+!  Sz2 = Sz*Sz + nonco_z
+!
+!  contam_ghf = 0.d0
+!  do j = 1, nO
+!    do i = 1, nO
+!      contam_ghf = contam_ghf - (Pab(i,i)*Pba(j,j) - Pab(i,j)*Pba(j,i))
+!    enddo
+!  enddo
+!  S2 = Sz * (Sz + 1.d0) + nonco_z + contam_ghf