diff --git a/input/methods b/input/methods
index d16be27..f7bba22 100644
--- a/input/methods
+++ b/input/methods
@@ -7,13 +7,13 @@
 # drCCD rCCD lCCD pCCD
   F      F   F    F
 # CIS* CIS(D) CID CISD
-  F    F      F   F    
+  T    F      F   F    
 # RPA* RPAx* ppRPA 
   F   F    F    
 # G0F2 evGF2 G0F3 evGF3
   F    F     F    F
 # G0W0* evGW* qsGW 
-  T    F    F
+  F    F    F
 # G0T0 evGT qsGT 
   F    F    F
 # MCMP2
diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90
index 53678e2..165330a 100644
--- a/src/HF/UHF.f90
+++ b/src/HF/UHF.f90
@@ -80,7 +80,8 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
 ! Guess coefficients and eigenvalues
 
   if(guess_type == 1) then
-
+    
+    F(:,:,:) = 0d0
     do ispin=1,nspin
       F(:,:,ispin) = Hc(:,:)
     end do
@@ -154,8 +155,9 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
     do ispin=1,nspin
       call exchange_matrix_AO_basis(nBas,P(:,:,ispin),ERI(:,:,:,:),K(:,:,ispin))
     end do
-
+ 
 !   Build Fock operator
+
     do ispin=1,nspin
       F(:,:,ispin) = Hc(:,:) + J(:,:,ispin) + J(:,:,mod(ispin,2)+1) + K(:,:,ispin)
     end do
@@ -166,7 +168,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
       err(:,:,ispin) = matmul(F(:,:,ispin),matmul(P(:,:,ispin),S(:,:))) - matmul(matmul(S(:,:),P(:,:,ispin)),F(:,:,ispin))
     end do
 
-    conv = maxval(abs(err(:,:,:)))
+    if(nSCF > 1) conv = maxval(abs(err(:,:,:)))
     
 !   DIIS extrapolation
 
diff --git a/src/LR/print_unrestricted_transition_vectors.f90 b/src/LR/print_unrestricted_transition_vectors.f90
index 7816192..a6e5d2d 100644
--- a/src/LR/print_unrestricted_transition_vectors.f90
+++ b/src/LR/print_unrestricted_transition_vectors.f90
@@ -49,7 +49,7 @@ subroutine print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,n
 
 ! Compute <S**2>
 
-  call unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,XpY,XmY,S2)
+  call unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,Omega,XpY,XmY,S2)
 
 ! Print details about spin-conserved excitations
 
@@ -60,7 +60,6 @@ subroutine print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,n
       X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:))
       Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
 
-
       print*,'-------------------------------------------------------------'
       write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & 
               ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV','  f = ',os(ia),'  <S**2> = ',S2(ia)
diff --git a/src/LR/unrestricted_S2_expval.f90 b/src/LR/unrestricted_S2_expval.f90
index c62bf2a..2389a98 100644
--- a/src/LR/unrestricted_S2_expval.f90
+++ b/src/LR/unrestricted_S2_expval.f90
@@ -1,4 +1,4 @@
-subroutine unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,XpY,XmY,S2)
+subroutine unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,Omega,XpY,XmY,S2)
 
 ! Compute <S**2> for linear response excited states
 
@@ -21,6 +21,7 @@ subroutine unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S
   integer,intent(in)            :: maxS
   double precision,intent(in)   :: c(nBas,nBas,nspin)
   double precision,intent(in)   :: S(nBas,nBas)
+  double precision,intent(in)   :: Omega(nSt)
   double precision,intent(in)   :: XpY(nSt,nSt)
   double precision,intent(in)   :: XmY(nSt,nSt)
 
@@ -152,7 +153,9 @@ subroutine unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S
       Yat(:,:) = transpose(Ya(:,:))
       Ybt(:,:) = transpose(Yb(:,:))
 
-      S2(m) = S2_gs + dble(nO(2) - nO(1)) + 1d0 &
+      S2(m) = S2_gs + dble(nO(2) - nO(1)) + 1d0  
+        
+      S2(m) = S2(m) &
 
             + trace_matrix(nV(1),matmul(Xbt,matmul(OOt,matmul(OO,Xb)))) &
             - trace_matrix(nO(2),matmul(Xb,matmul(VO,matmul(VOt,Xbt)))) &