diff --git a/input/methods b/input/methods
index f582e3d..faa85b7 100644
--- a/input/methods
+++ b/input/methods
@@ -1,5 +1,5 @@
 # RHF UHF GHF ROHF 
-  F   F   T   F   
+  F   T   F   F   
 # MP2 MP3 
   F   F   
 # CCD pCCD DCD CCSD CCSD(T) 
diff --git a/input/options b/input/options
index 54c59ac..f900f68 100644
--- a/input/options
+++ b/input/options
@@ -1,5 +1,5 @@
 # HF: maxSCF thresh  DIIS guess mix   shift stab search
-      1000  0.0000001  5     1    0.0   0.0    T    T
+      10000  0.00001  5     1    0.0   1.0    F    T
 # MP: reg
       F
 # CC: maxSCF thresh     DIIS
diff --git a/src/HF/RHF_search.f90 b/src/HF/RHF_search.f90
index d27efa4..127b89e 100644
--- a/src/HF/RHF_search.f90
+++ b/src/HF/RHF_search.f90
@@ -159,30 +159,30 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
 
       step = 1d0
 
-      do mu=1,nBas
-        ia = 0
-        do i=nC+1,nO
-          kick = 0d0
-          do a=nO+1,nBas-nR
-            ia = ia + 1
-            kick = kick + AB(ia,eig)*c(mu,a)
-          end do
-          c(mu,i) = c(mu,i) + step*kick
-        end do
-      end do
-
-!     R(:,:) = 0d0
-!     ia = 0
-!     do i=nC+1,nO
-!       do a=nO+1,nBas-nR
-!         ia = ia + 1
-!         R(a,i) = +AB(ia,eig)
-!         R(i,a) = -AB(ia,eig)
+!     do mu=1,nBas
+!       ia = 0
+!       do i=nC+1,nO
+!         kick = 0d0
+!         do a=nO+1,nBas-nR
+!           ia = ia + 1
+!           kick = kick + AB(ia,eig)*c(mu,a)
+!         end do
+!         c(mu,i) = c(mu,i) + step*kick
 !       end do
 !     end do
 
-!     call matrix_exponential(nBas,R,ExpR)
-!     c = matmul(c,ExpR)
+      R(:,:) = 0d0
+      ia = 0
+      do i=nC+1,nO
+        do a=nO+1,nBas-nR
+          ia = ia + 1
+          R(a,i) = +AB(ia,eig)
+          R(i,a) = -AB(ia,eig)
+        end do
+      end do
+
+      call matrix_exponential(nBas,R,ExpR)
+      c = matmul(c,ExpR)
 
     else
  
diff --git a/src/HF/UHF_search.f90 b/src/HF/UHF_search.f90
index ecbf1ce..511dbcb 100644
--- a/src/HF/UHF_search.f90
+++ b/src/HF/UHF_search.f90
@@ -53,6 +53,8 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
   double precision,allocatable  :: A_sc(:,:)
   double precision,allocatable  :: B_sc(:,:)
   double precision,allocatable  :: AB_sc(:,:)
+  double precision,allocatable  :: R(:,:)
+  double precision,allocatable  :: ExpR(:,:)
   
   integer                       :: eig
   double precision              :: kick,step
@@ -83,7 +85,7 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
   nS_sc = nS_aa + nS_bb
 
   allocate(ERI_aaaa(nBas,nBas,nBas,nBas),ERI_aabb(nBas,nBas,nBas,nBas),ERI_bbbb(nBas,nBas,nBas,nBas))
-  allocate(Om_sc(nS_sc),A_sc(nS_sc,nS_sc),B_sc(nS_sc,nS_sc),AB_sc(nS_sc,nS_sc))
+  allocate(Om_sc(nS_sc),A_sc(nS_sc,nS_sc),B_sc(nS_sc,nS_sc),AB_sc(nS_sc,nS_sc),R(nBas,nBas),ExpR(nBas,nBas))
 
 !------------------!
 ! Search algorithm !
@@ -179,31 +181,57 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
 
       ! Spin-up kick
 
-      do mu=1,nBas
-        ia = 0
-        do i=nC(1)+1,nO(1)
-          kick = 0d0
-          do a=nO(1)+1,nBas-nR(1)
-            ia = ia + 1
-            kick = kick + AB_sc(ia,eig)*c(mu,a,1)
-          end do
-          c(mu,i,1) = c(mu,i,1) + step*kick
+!     do mu=1,nBas
+!       ia = 0
+!       do i=nC(1)+1,nO(1)
+!         kick = 0d0
+!         do a=nO(1)+1,nBas-nR(1)
+!           ia = ia + 1
+!           kick = kick + AB_sc(ia,eig)*c(mu,a,1)
+!         end do
+!         c(mu,i,1) = c(mu,i,1) + step*kick
+!       end do
+!     end do
+
+      R(:,:) = 0d0
+      ia = 0
+      do i=nC(1)+1,nO(1)
+        do a=nO(1)+1,nBas-nR(1)
+          ia = ia + 1
+          R(a,i) = +AB_sc(ia,eig)
+          R(i,a) = -AB_sc(ia,eig)
         end do
       end do
 
+      call matrix_exponential(nBas,R,ExpR)
+      c(:,:,1) = matmul(c(:,:,1),ExpR)
+
       ! Spin-down kick
 
-      do mu=1,nBas
-        ia = nS_aa
-        do i=nC(2)+1,nO(2)
-          kick = 0d0
-          do a=nO(2)+1,nBas-nR(2)
-            ia = ia + 1
-            kick = kick + AB_sc(ia,eig)*c(mu,a,2)
-          end do
-          c(mu,i,2) = c(mu,i,2) + step*kick
+!     do mu=1,nBas
+!       ia = nS_aa
+!       do i=nC(2)+1,nO(2)
+!         kick = 0d0
+!         do a=nO(2)+1,nBas-nR(2)
+!           ia = ia + 1
+!           kick = kick + AB_sc(ia,eig)*c(mu,a,2)
+!         end do
+!         c(mu,i,2) = c(mu,i,2) + step*kick
+!       end do
+!     end do
+ 
+      R(:,:) = 0d0
+      ia = nS_aa
+      do i=nC(2)+1,nO(2)
+        do a=nO(2)+1,nBas-nR(2)
+          ia = ia + 1
+          R(a,i) = +AB_sc(ia,eig)
+          R(i,a) = -AB_sc(ia,eig)
         end do
       end do
+
+      call matrix_exponential(nBas,R,ExpR)
+      c(:,:,2) = matmul(c(:,:,2),ExpR)
  
     else
  
diff --git a/src/utils/utils.f90 b/src/utils/utils.f90
index 345836c..961b419 100644
--- a/src/utils/utils.f90
+++ b/src/utils/utils.f90
@@ -92,12 +92,10 @@ subroutine matrix_exponential(N,A,ExpA)
   W(:,:) = - matmul(A,A)
   call diagonalize_matrix(N,W,tau)
 
-  call matout(N,1,tau)
 ! do i=1,N
 !   tau(i) = max(abs(tau(i)),1d-14)
 ! end do
   tau(:) = sqrt(abs(tau(:)))
-  call matout(N,1,tau)
 
 ! Construct cos part