diff --git a/src/HF/print_GHF.f90 b/src/HF/print_GHF.f90
index 8fde002..ef60fcd 100644
--- a/src/HF/print_GHF.f90
+++ b/src/HF/print_GHF.f90
@@ -37,6 +37,15 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole)
   double precision,allocatable       :: Pba(:,:)
   double precision,allocatable       :: Pbb(:,:)
 
+  double precision,allocatable       :: Mx(:,:)
+  double precision,allocatable       :: My(:,:)
+  double precision,allocatable       :: Mz(:,:)
+  double precision,allocatable       :: PP(:,:)
+  double precision                   :: T(3,3)
+  double precision                   :: vec(3,3)
+  double precision                   :: val(3)
+  double precision                   :: lambda
+
   double precision,external          :: trace_matrix
 
   logical                            :: dump_orb = .false.
@@ -49,43 +58,82 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole)
 
 ! Density matrices
 
-  allocate(Paa(nBas2,nBas2),Pab(nBas2,nBas2),Pba(nBas2,nBas2),Pbb(nBas2,nBas2))
+  allocate(Paa(nBas,nBas),Pab(nBas,nBas),Pba(nBas,nBas),Pbb(nBas,nBas))
 
   Paa(:,:) = P(     1:nBas ,     1:nBas )
   Pab(:,:) = P(     1:nBas ,nBas+1:nBas2)
   Pba(:,:) = P(nBas+1:nBas2,     1:nBas )
   Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2)
 
-  allocate(Ca(nBas,nBas2),Cb(nBas,nBas2))
+! allocate(Ca(nBas,nBas2),Cb(nBas,nBas2))
 
-  Ca(:,:) = C(     1:nBas ,1:nBas2)
-  Cb(:,:) = C(nBas+1:nBas2,1:nBas2)
+! Ca(:,:) = C(     1:nBas ,1:nBas2)
+! Cb(:,:) = C(nBas+1:nBas2,1:nBas2)
 
 ! Compute expectation values of S^2 (WRONG!)
 
-  Sx2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) + 0.25d0*trace_matrix(nBas,Pab+Pba)**2
-  do mu=1,nBas
-    do nu=1,nBas
-        Sx2 = Sx2 - 0.5d0*(Paa(mu,nu)*Pbb(nu,mu) + Pab(mu,nu)*Pab(nu,mu))
-    end do
-  end do
+! Sx2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) + 0.25d0*trace_matrix(nBas,Pab+Pba)**2
+! do mu=1,nBas
+!   do nu=1,nBas
+!       Sx2 = Sx2 - 0.5d0*(Paa(mu,nu)*Pbb(nu,mu) + Pab(mu,nu)*Pab(nu,mu))
+!   end do
+! end do
 
-  Sy2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) - 0.25d0*trace_matrix(nBas,Pab+Pba)**2
-  do mu=1,nBas
-    do nu=1,nBas
-        Sy2 = Sy2 - 0.5d0*(Paa(mu,nu)*Pbb(nu,mu) - Pab(mu,nu)*Pab(nu,mu))
-    end do
-  end do
+! Sy2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) - 0.25d0*trace_matrix(nBas,Pab+Pba)**2
+! do mu=1,nBas
+!   do nu=1,nBas
+!       Sy2 = Sy2 - 0.5d0*(Paa(mu,nu)*Pbb(nu,mu) - Pab(mu,nu)*Pab(nu,mu))
+!   end do
+! end do
 
-  Sz2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) + 0.25d0*trace_matrix(nBas,Pab-Pba)**2
-  do mu=1,nBas
-    do nu=1,nBas
-        Sz2 = Sz2 - 0.25d0*(Paa(mu,nu)*Pbb(nu,mu) - Pab(mu,nu)*Pab(nu,mu))
-        Sz2 = Sz2 + 0.25d0*(Pab(mu,nu)*Pba(nu,mu) - Pba(mu,nu)*Pab(nu,mu))
-    end do
-  end do
-  
-  S2 = Sx2 + Sy2 + Sz2
+! Sz2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) + 0.25d0*trace_matrix(nBas,Pab-Pba)**2
+! do mu=1,nBas
+!   do nu=1,nBas
+!       Sz2 = Sz2 - 0.25d0*(Paa(mu,nu)*Pbb(nu,mu) - Pab(mu,nu)*Pab(nu,mu))
+!       Sz2 = Sz2 + 0.25d0*(Pab(mu,nu)*Pba(nu,mu) - Pba(mu,nu)*Pab(nu,mu))
+!   end do
+! end do
+! 
+! S2 = Sx2 + Sy2 + Sz2
+
+! Checl collinearity and coplanarity 
+
+  allocate(PP(nBas,nBas),Mx(nBas,nBas),My(nBas,nBas),Mz(nBas,nBas))
+
+  PP(:,:) = 0.5d0*(Paa(:,:) + Pbb(:,:))
+  Mx(:,:) = 0.5d0*(Pba(:,:) + Pab(:,:))
+  My(:,:) = 0.5d0*(Pba(:,:) - Pab(:,:))
+  Mz(:,:) = 0.5d0*(Paa(:,:) - Pbb(:,:))
+
+  T(1,1) = trace_matrix(nBas,matmul(Mx,transpose(Mx)))
+  T(1,2) = trace_matrix(nBas,matmul(Mx,transpose(My)))
+  T(1,3) = trace_matrix(nBas,matmul(Mx,transpose(Mz)))
+  T(2,1) = trace_matrix(nBas,matmul(My,transpose(Mx)))
+  T(2,2) = trace_matrix(nBas,matmul(My,transpose(My)))
+  T(2,3) = trace_matrix(nBas,matmul(My,transpose(Mz)))
+  T(3,1) = trace_matrix(nBas,matmul(Mz,transpose(Mx)))
+  T(3,2) = trace_matrix(nBas,matmul(Mz,transpose(My)))
+  T(3,3) = trace_matrix(nBas,matmul(Mz,transpose(Mz)))
+
+  print*,'Value of Tr(P - P^2)'
+  lambda = trace_matrix(nBas,PP - matmul(PP,transpose(PP)))
+  print*,lambda
+
+  print*,'Eigenvalues of T'
+  vec(:,:) = T(:,:)
+  call diagonalize_matrix(3,vec,val)
+  print*,val
+
+  T(1,1) = - T(1,1) + lambda
+  T(2,2) = - T(2,2) + lambda 
+  T(3,3) = - T(3,3) + lambda
+
+  print*,'Eigenvalues of A'
+  vec(:,:) = T(:,:)
+  call diagonalize_matrix(3,vec,val)
+  print*,val
+
+  deallocate(PP,Mx,My,Mz)
 
 ! Dump results