diff --git a/input/methods b/input/methods
index 4d84f78..f032a00 100644
--- a/input/methods
+++ b/input/methods
@@ -13,7 +13,7 @@
 # G0F2* evGF2* qsGF2* G0F3 evGF3
   F     F      F      F    F
 # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW
-  T     F     F     F        F      F
+  F     F     F     F        F      F
 # G0T0pp* evGTpp* qsGTpp* G0T0eh evGTeh qsGTeh
   F       F       F       F      F      F
 # * unrestricted version available
diff --git a/src/HF/HF.f90 b/src/HF/HF.f90
index 1b62bda..e549bd3 100644
--- a/src/HF/HF.f90
+++ b/src/HF/HF.f90
@@ -39,34 +39,15 @@ subroutine HF(doRHF,doUHF,doROHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_di
 
   double precision              :: start_HF     ,end_HF       ,t_HF
 
-  integer                       :: nSCF
-  double precision              :: ET
-  double precision              :: EV
-  double precision              :: EJ
-  double precision              :: EK
-  double precision              :: dipole(ncart)
-
-  double precision              :: Conv
-  double precision              :: Gap 
-  double precision              :: rcond
-  double precision,external     :: trace_matrix
-  double precision,allocatable  :: error(:,:)
-  double precision,allocatable  :: error_diis(:,:)
-  double precision,allocatable  :: F_diis(:,:)
-  double precision,allocatable  :: J(:,:)
-  double precision,allocatable  :: K(:,:)
-  double precision,allocatable  :: cp(:,:)
-  double precision,allocatable  :: Fp(:,:)
-
 ! Output variables
 
   logical,intent(out)           :: unrestricted
 
   double precision,intent(out)  :: EHF
-  double precision,intent(out)  :: epsHF(nBas)
-  double precision,intent(out)  :: cHF(nBas,nBas)
-  double precision,intent(out)  :: PHF(nBas,nBas)
-  double precision,intent(out)  :: F(nBas,nBas)
+  double precision,intent(out)  :: epsHF(nBas,nspin)
+  double precision,intent(out)  :: cHF(nBas,nBas,nspin)
+  double precision,intent(out)  :: PHF(nBas,nBas,nspin)
+  double precision,intent(out)  :: F(nBas,nBas,nspin)
 
 !------------------------------------------------------------------------
 ! Compute RHF energy
diff --git a/src/HF/ROHF.f90 b/src/HF/ROHF.f90
index 7c22a26..69a2701 100644
--- a/src/HF/ROHF.f90
+++ b/src/HF/ROHF.f90
@@ -239,6 +239,6 @@ subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc
 ! Compute final UHF energy
 
   call dipole_moment(nBas,Ptot,nNuc,ZNuc,rNuc,dipole_int,dipole)
-  call print_ROHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole)
+  call print_ROHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole)
 
 end subroutine 
diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90
index 4f2a900..a15dacb 100644
--- a/src/HF/UHF.f90
+++ b/src/HF/UHF.f90
@@ -112,35 +112,6 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
 
     nSCF = nSCF + 1
 
-!  Transform Fock matrix in orthogonal basis
-
-    do ispin=1,nspin
-      Fp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(F(:,:,ispin),X(:,:)))
-    end do
-
-!  Diagonalize Fock matrix to get eigenvectors and eigenvalues
-
-    cp(:,:,:) = Fp(:,:,:)
-    do ispin=1,nspin
-      call diagonalize_matrix(nBas,cp(:,:,ispin),e(:,ispin))
-    end do
-    
-!   Back-transform eigenvectors in non-orthogonal basis
-
-    do ispin=1,nspin
-      c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin))
-    end do
-
-!   Mix guess for UHF solution in singlet states
-
-    if(mix .and. nSCF == 1) call mix_guess(nBas,nO,c)
-
-!   Compute density matrix 
-
-    do ispin=1,nspin
-      P(:,:,ispin) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin)))
-    end do
- 
 !   Build Coulomb repulsion
 
     do ispin=1,nspin
@@ -189,6 +160,35 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
 
     end if
 
+!  Transform Fock matrix in orthogonal basis
+
+    do ispin=1,nspin
+      Fp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(F(:,:,ispin),X(:,:)))
+    end do
+
+!  Diagonalize Fock matrix to get eigenvectors and eigenvalues
+
+    cp(:,:,:) = Fp(:,:,:)
+    do ispin=1,nspin
+      call diagonalize_matrix(nBas,cp(:,:,ispin),e(:,ispin))
+    end do
+    
+!   Back-transform eigenvectors in non-orthogonal basis
+
+    do ispin=1,nspin
+      c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin))
+    end do
+
+!   Mix guess for UHF solution in singlet states
+
+    if(mix .and. nSCF == 1) call mix_guess(nBas,nO,c)
+
+!   Compute density matrix 
+
+    do ispin=1,nspin
+      P(:,:,ispin) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin)))
+    end do
+ 
 !------------------------------------------------------------------------
 !   Compute UHF energy
 !------------------------------------------------------------------------
diff --git a/src/HF/mo_guess.f90 b/src/HF/mo_guess.f90
index 387f371..9822334 100644
--- a/src/HF/mo_guess.f90
+++ b/src/HF/mo_guess.f90
@@ -12,10 +12,6 @@ subroutine mo_guess(nBas,guess_type,S,Hc,X,c)
   double precision,intent(in)   :: Hc(nBas,nBas)
   double precision,intent(in)   :: X(nBas,nBas)
 
-! Local variables
-
-  integer                       :: nSCF
-
 ! Output variables
 
   double precision,intent(out)  :: c(nBas,nBas)
diff --git a/src/HF/print_ROHF.f90 b/src/HF/print_ROHF.f90
index eec0196..acfc694 100644
--- a/src/HF/print_ROHF.f90
+++ b/src/HF/print_ROHF.f90
@@ -1,4 +1,4 @@
-subroutine print_ROHF(nBas,nO,Ov,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole)
+subroutine print_ROHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole)
 
 ! Print one- and two-electron energies and other stuff for RoHF calculation
 
@@ -7,7 +7,6 @@ subroutine print_ROHF(nBas,nO,Ov,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole)
 
   integer,intent(in)                 :: nBas
   integer,intent(in)                 :: nO(nspin)
-  double precision,intent(in)        :: Ov(nBas,nBas)
   double precision,intent(in)        :: e(nBas)
   double precision,intent(in)        :: c(nBas,nBas)
   double precision,intent(in)        :: ENuc
@@ -23,7 +22,6 @@ subroutine print_ROHF(nBas,nO,Ov,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole)
   double precision                   :: HOMO(nspin)
   double precision                   :: LUMO(nspin)
   double precision                   :: Gap(nspin)
-  double precision                   :: S_exact,S2_exact
   double precision                   :: S,S2
 
 ! HOMO and LUMO
diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90
index 40d2e79..7219996 100644
--- a/src/QuAcK/QuAcK.f90
+++ b/src/QuAcK/QuAcK.f90
@@ -42,8 +42,8 @@ program QuAcK
   double precision,allocatable  :: dipole_int_MO(:,:,:)
   double precision,allocatable  :: dipole_int_aa(:,:,:)
   double precision,allocatable  :: dipole_int_bb(:,:,:)
-  double precision,allocatable  :: F_AO(:,:)
-  double precision,allocatable  :: F_MO(:,:)
+  double precision,allocatable  :: F_AO(:,:,:)
+  double precision,allocatable  :: F_MO(:,:,:)
   double precision,allocatable  :: ERI_AO(:,:,:,:)
   double precision,allocatable  :: ERI_MO(:,:,:,:)
   integer                       :: ixyz
@@ -185,7 +185,7 @@ program QuAcK
 
   allocate(cHF(nBas,nBas,nspin),epsHF(nBas,nspin),PHF(nBas,nBas,nspin),S(nBas,nBas),T(nBas,nBas),              &
            V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas),dipole_int_AO(nBas,nBas,ncart), & 
-           dipole_int_MO(nBas,nBas,ncart),F_AO(nBas,nBas))
+           dipole_int_MO(nBas,nBas,ncart),F_AO(nBas,nBas,nspin))
 
 ! Read integrals
 
@@ -215,8 +215,8 @@ program QuAcK
 
     call wall_time(start_HF)
     call HF(doRHF,doUHF,doROHF,doRMOM,doUMOM,unrestricted,maxSCF_HF,thresh_HF,max_diis_HF, & 
-            guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F_AO,          & 
-            ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
+            guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F_AO,ERI_AO,   & 
+            dipole_int_AO,X,EHF,epsHF,cHF,PHF)
     call wall_time(end_HF)
 
     t_HF = end_HF - start_HF
@@ -291,7 +291,7 @@ program QuAcK
 
     ! Memory allocation
    
-    allocate(ERI_MO(nBas,nBas,nBas,nBas),F_MO(nBas,nBas))
+    allocate(ERI_MO(nBas,nBas,nBas,nBas),F_MO(nBas,nBas,nspin))
 
     ! Read and transform dipole-related integrals
   
@@ -304,7 +304,7 @@ program QuAcK
    
     call AOtoMO_integral_transform(1,1,1,1,nBas,cHF,ERI_AO,ERI_MO)
 
-    F_MO(:,:) = F_AO(:,:)
+    F_MO(:,:,1) = F_AO(:,:,1)
     call AOtoMO_transform(nBas,cHF,F_MO)
 
   end if
diff --git a/src/make_ninja.py b/src/make_ninja.py
index 6c7a58d..4022c0f 100755
--- a/src/make_ninja.py
+++ b/src/make_ninja.py
@@ -57,7 +57,7 @@ else:
     compile_gfortran_mac = """
 FC = gfortran
 AR = libtool -static -o
-FFLAGS = -I$IDIR -J$IDIR -fbacktrace -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant
+FFLAGS = -I$IDIR -J$IDIR -fbacktrace -Wall -g -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant
 CC = gcc
 CXX = g++
 LAPACK=-lblas -llapack
@@ -68,7 +68,7 @@ FIX_ORDER_OF_LIBS=
     compile_gfortran_linux = """
 FC = gfortran
 AR = ar crs
-FFLAGS = -I$IDIR -J$IDIR -fbacktrace -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant
+FFLAGS = -I$IDIR -J$IDIR -fbacktrace -Wall -g -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant
 CC = gcc
 CXX = g++
 LAPACK=-lblas -llapack