diff --git a/input/methods b/input/methods
index faa85b7..41e9098 100644
--- a/input/methods
+++ b/input/methods
@@ -1,5 +1,5 @@
 # RHF UHF GHF ROHF 
-  F   T   F   F   
+  T   F   F   F   
 # MP2 MP3 
   F   F   
 # CCD pCCD DCD CCSD CCSD(T) 
diff --git a/src/GW/GGW.f90 b/src/GW/GGW.f90
index af1b4f0..3711831 100644
--- a/src/GW/GGW.f90
+++ b/src/GW/GGW.f90
@@ -105,9 +105,9 @@ subroutine GGW(doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,      &
   if(doqsGW) then 
 
     call wall_time(start_GW)
-!   call qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & 
-!             eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,  & 
-!             dipole_int_AO,dipole_int,PHF,cHF,epsHF)
+!   call qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & 
+!              eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,  & 
+!              dipole_int_AO,dipole_int,PHF,cHF,epsHF)
     call wall_time(end_GW)
 
     t_GW = end_GW - start_GW
diff --git a/src/GW/qsUGW.f90 b/src/GW/qsUGW.f90
index 9c43f4c..b3c7460 100644
--- a/src/GW/qsUGW.f90
+++ b/src/GW/qsUGW.f90
@@ -57,8 +57,6 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
 
 ! Local variables
 
-  logical                       :: doGWPT = .false.
-
   integer                       :: nSCF
   integer                       :: nBasSq
   integer                       :: ispin
diff --git a/src/HF/RHF_search.f90 b/src/HF/RHF_search.f90
index 144d47e..2c00bd2 100644
--- a/src/HF/RHF_search.f90
+++ b/src/HF/RHF_search.f90
@@ -48,6 +48,8 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
   double precision,allocatable  :: Bph(:,:)
   double precision,allocatable  :: AB(:,:)
   double precision,allocatable  :: Om(:)
+  double precision,allocatable  :: R(:,:)
+  double precision,allocatable  :: ExpR(:,:)
   
   integer                       :: eig
   double precision              :: kick,step
@@ -72,7 +74,8 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
 !-------------------!
 
   nS = (nO - nC)*(nV - nR)
-  allocate(ERI_MO(nBas,nBas,nBas,nBas),Aph(nS,nS),Bph(nS,nS),AB(nS,nS),Om(nS))
+  allocate(ERI_MO(nBas,nBas,nBas,nBas),Aph(nS,nS),Bph(nS,nS),AB(nS,nS),Om(nS), &
+           R(nBas,nBas),ExpR(nBas,nBas))
 
 !------------------!
 ! Search algorithm !
@@ -124,7 +127,7 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
  
     call diagonalize_matrix(nS,AB,Om)
     Om(:) = 0.5d0*Om(:)
- 
+
     write(*,*)'-------------------------------------------------------------'
     write(*,*)'|       Stability analysis: Real RHF -> Real RHF            |'
     write(*,*)'-------------------------------------------------------------'
@@ -168,6 +171,19 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
         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)
+!       end do
+!     end do
+
+!     call matrix_exponential(nBas,R,ExpR)
+!     c = matmul(c,ExpR)
+
     else
  
       write(*,'(1X,A40,1X)')           'Well done, RHF solution is stable!'
diff --git a/src/utils/utils.f90 b/src/utils/utils.f90
index 6832245..f1e26e4 100644
--- a/src/utils/utils.f90
+++ b/src/utils/utils.f90
@@ -44,6 +44,69 @@ function KroneckerDelta(i,j) result(delta)
 
 end function 
 
+!------------------------------------------------------------------------
+subroutine diagonal_matrix(N,D,A)
+
+! Construct diagonal matrix A from vector D
+
+  implicit none
+
+  integer,intent(in)            :: N
+  double precision,intent(in)   :: D(N)
+  double precision,intent(out)  :: A(N,N)
+
+  integer                       :: i
+
+  A(:,:) = 0d0
+  do i=1,N
+    A(i,i) = D(i)
+  end do
+
+end subroutine
+
+!------------------------------------------------------------------------
+subroutine matrix_exponential(N,A,ExpA)
+
+! Compute Exp(A)
+
+  implicit none
+
+  integer,intent(in)            :: N
+  double precision,intent(in)   :: A(N,N)
+  double precision,allocatable  :: W(:,:)
+  double precision,allocatable  :: tau(:)
+  double precision,allocatable  :: t(:,:)
+  double precision,intent(out)  :: ExpA(N,N)
+
+! Memory allocation
+
+  allocate(W(N,N),tau(N),t(N,N))
+
+! Initialize
+
+  ExpA(:,:) = 0d0
+
+! Diagonalize
+
+  W(:,:) = - matmul(A,A)
+  call diagonalize_matrix(N,W,tau)
+  tau(:) = sqrt(tau(:))
+
+! Construct cos part
+
+  call diagonal_matrix(N,cos(tau),t)
+  t(:,:) = matmul(t,transpose(W))
+  ExpA(:,:) = ExpA(:,:) + matmul(W,t) 
+
+! Construct sin part
+
+  call diagonal_matrix(N,sin(tau)/tau,t)
+  t(:,:) = matmul(t,transpose(W))
+  t(:,:) = matmul(t,A)
+  ExpA(:,:) = ExpA(:,:) + matmul(W,t)
+
+end subroutine
+
 !------------------------------------------------------------------------
 subroutine matout(m,n,A)