diff --git a/input/methods b/input/methods
index 91eee5d..b71abbf 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 1f2fdc6..e86f9b7 100644
--- a/input/options
+++ b/input/options
@@ -1,5 +1,5 @@
 # HF: maxSCF thresh    DIIS guess mix_guess level_shift stability
-      5000   0.0000001  5     1   0.0       0.0         T
+      10000   0.000001  5     2    0.0       0.0         T
 # MP: reg
       F
 # CC: maxSCF thresh     DIIS
diff --git a/src/AOtoMO/AOtoMO_integral_transform_GHF.f90 b/src/AOtoMO/AOtoMO_integral_transform_GHF.f90
new file mode 100644
index 0000000..e11e689
--- /dev/null
+++ b/src/AOtoMO/AOtoMO_integral_transform_GHF.f90
@@ -0,0 +1,90 @@
+subroutine AOtoMO_integral_transform_GHF(nBas,nBas2,c1,c2,c3,c4,ERI_AO_basis,ERI_MO_basis)
+
+! AO to MO transformation of two-electron integrals via the semi-direct O(N^5) algorithm
+! bra and ket are the spin of (bra1 bra2|ket1 ket2)
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nBas2
+  double precision,intent(in)   :: ERI_AO_basis(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: c1(nBas,nBas2)
+  double precision,intent(in)   :: c2(nBas,nBas2)
+  double precision,intent(in)   :: c3(nBas,nBas2)
+  double precision,intent(in)   :: c4(nBas,nBas2)
+
+! Local variables
+
+  double precision,allocatable  :: scr(:,:,:,:)
+  integer                       :: mu,nu,la,si,i,j,k,l
+
+! Output variables
+
+  double precision,intent(out)  :: ERI_MO_basis(nBas2,nBas2,nBas2,nBas2)
+
+! Memory allocation
+
+  allocate(scr(nBas2,nBas2,nBas2,nBas2))
+
+! Four-index transform via semi-direct O(N^5) algorithm
+
+  scr(:,:,:,:) = 0d0
+
+  do l=1,nBas2
+    do si=1,nBas
+      do la=1,nBas
+        do nu=1,nBas
+          do mu=1,nBas
+            scr(mu,nu,la,l) = scr(mu,nu,la,l) + ERI_AO_basis(mu,nu,la,si)*c4(si,l)
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+  ERI_MO_basis(:,:,:,:) = 0d0
+
+  do l=1,nBas2
+    do la=1,nBas
+      do nu=1,nBas
+        do i=1,nBas2
+          do mu=1,nBas
+            ERI_MO_basis(i,nu,la,l) = ERI_MO_basis(i,nu,la,l) + c1(mu,i)*scr(mu,nu,la,l)
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+  scr(:,:,:,:) = 0d0
+
+  do l=1,nBas2
+    do k=1,nBas2
+      do la=1,nBas
+        do nu=1,nBas
+          do i=1,nBas2
+            scr(i,nu,k,l) = scr(i,nu,k,l) + ERI_MO_basis(i,nu,la,l)*c3(la,k)
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+  ERI_MO_basis(:,:,:,:) = 0d0
+
+  do l=1,nBas2
+    do k=1,nBas2
+      do j=1,nBas2
+        do i=1,nBas2
+          do nu=1,nBas
+            ERI_MO_basis(i,j,k,l) = ERI_MO_basis(i,j,k,l) + c2(nu,j)*scr(i,nu,k,l)
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+end subroutine 
diff --git a/src/GW/GG0W0.f90 b/src/GW/GG0W0.f90
new file mode 100644
index 0000000..6aae026
--- /dev/null
+++ b/src/GW/GG0W0.f90
@@ -0,0 +1,231 @@
+subroutine GG0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & 
+                 linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
+! Perform G0W0 calculation
+
+  implicit none
+  include 'parameters.h'
+  include 'quadrature.h'
+
+! Input variables
+
+  logical,intent(in)            :: doACFDT
+  logical,intent(in)            :: exchange_kernel
+  logical,intent(in)            :: doXBS
+  logical,intent(in)            :: dophBSE
+  logical,intent(in)            :: dophBSE2
+  logical,intent(in)            :: doppBSE
+  logical,intent(in)            :: TDA_W
+  logical,intent(in)            :: TDA
+  logical,intent(in)            :: dBSE
+  logical,intent(in)            :: dTDA
+  logical,intent(in)            :: linearize
+  double precision,intent(in)   :: eta
+  logical,intent(in)            :: regularize
+
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: ENuc
+  double precision,intent(in)   :: ERHF
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
+  double precision,intent(in)   :: eHF(nBas)
+
+! Local variables
+
+  logical                       :: print_W = .true.
+  logical                       :: dRPA
+  integer                       :: ispin
+  double precision              :: EcRPA
+  double precision              :: EcBSE(nspin)
+  double precision              :: EcGM
+  double precision,allocatable  :: Aph(:,:)
+  double precision,allocatable  :: Bph(:,:)
+  double precision,allocatable  :: SigC(:)
+  double precision,allocatable  :: Z(:)
+  double precision,allocatable  :: Om(:)
+  double precision,allocatable  :: XpY(:,:)
+  double precision,allocatable  :: XmY(:,:)
+  double precision,allocatable  :: rho(:,:,:)
+
+  double precision,allocatable  :: eGWlin(:)
+  double precision,allocatable  :: eGW(:)
+
+! Output variables
+
+! Hello world
+
+  write(*,*)
+  write(*,*)'************************************************'
+  write(*,*)'|          One-shot G0W0 calculation           |'
+  write(*,*)'************************************************'
+  write(*,*)
+
+! Initialization
+
+  dRPA = .true.
+  EcRPA = 0d0
+
+! TDA for W
+
+  if(TDA_W) then 
+    write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
+    write(*,*)
+  end if
+
+! TDA 
+
+  if(TDA) then 
+    write(*,*) 'Tamm-Dancoff approximation activated!'
+    write(*,*)
+  end if
+
+! Spin manifold 
+
+  ispin = 3
+
+! Memory allocation
+
+  allocate(Aph(nS,nS),Bph(nS,nS),SigC(nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), & 
+           eGW(nBas),eGWlin(nBas))
+
+!-------------------!
+! Compute screening !
+!-------------------!
+
+                 call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
+  if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+
+  call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+
+  if(print_W) call print_excitation_energies('phRPA@GHF',ispin,nS,Om)
+
+!--------------------------!
+! Compute spectral weights !
+!--------------------------!
+
+  call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho)
+
+!------------------------!
+! Compute GW self-energy !
+!------------------------!
+
+  if(regularize) call GW_regularization(nBas,nC,nO,nV,nR,nS,eHF,Om,rho)
+
+  call GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
+
+!-----------------------------------!
+! Solve the quasi-particle equation !
+!-----------------------------------!
+
+  ! Linearized or graphical solution?
+
+  eGWlin(:) = eHF(:) + Z(:)*SigC(:)
+
+  if(linearize) then 
+ 
+    write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
+    write(*,*)
+
+    eGW(:) = eGWlin(:)
+
+  else 
+
+    write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
+    write(*,*)
+  
+    call GGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eHF,eGW,Z)
+
+  end if
+
+! call GW_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eHF,Om,rho)
+
+! Compute the RPA correlation energy
+
+  call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
+  if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+
+  call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+
+!--------------!
+! Dump results !
+!--------------!
+
+  call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM)
+
+! Deallocate memory
+
+  deallocate(SigC,Z,Om,XpY,XmY,rho)
+
+! Perform BSE calculation
+
+! if(dophBSE) then
+
+!   call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
+
+!   if(exchange_kernel) then
+!
+!     EcBSE(1) = 0.5d0*EcBSE(1)
+!     EcBSE(2) = 1.5d0*EcBSE(2)
+!
+!   end if
+
+!   write(*,*)
+!   write(*,*)'-------------------------------------------------------------------------------'
+!   write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 correlation energy (singlet) =',EcBSE(1),' au'
+!   write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 correlation energy (triplet) =',EcBSE(2),' au'
+!   write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 correlation energy           =',EcBSE(1) + EcBSE(2),' au'
+!   write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 total energy                 =',ENuc + ERHF + EcBSE(1) + EcBSE(2),' au'
+!   write(*,*)'-------------------------------------------------------------------------------'
+!   write(*,*)
+
+!   Compute the BSE correlation energy via the adiabatic connection 
+
+!   if(doACFDT) then
+
+!     write(*,*) '-------------------------------------------------------------'
+!     write(*,*) ' Adiabatic connection version of BSE@G0W0 correlation energy '
+!     write(*,*) '-------------------------------------------------------------'
+!     write(*,*) 
+
+!     if(doXBS) then 
+
+!       write(*,*) '*** scaled screening version (XBS) ***'
+!       write(*,*)
+
+!     end if
+
+!     call GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,EcBSE)
+
+!     write(*,*)
+!     write(*,*)'-------------------------------------------------------------------------------'
+!     write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy (singlet) =',EcBSE(1),' au'
+!     write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy (triplet) =',EcBSE(2),' au'
+!     write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy           =',EcBSE(1) + EcBSE(2),' au'
+!     write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 total energy                 =',ENuc + ERHF + EcBSE(1) + EcBSE(2),' au'
+!     write(*,*)'-------------------------------------------------------------------------------'
+!     write(*,*)
+
+!   end if
+
+! end if
+
+! if(doppBSE) then
+
+!   call GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
+
+!   write(*,*)
+!   write(*,*)'-------------------------------------------------------------------------------'
+!   write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (singlet) =',EcBSE(1),' au'
+!   write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (triplet) =',3d0*EcBSE(2),' au'
+!   write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy           =',EcBSE(1) + 3d0*EcBSE(2),' au'
+!   write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 total energy                 =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2),' au'
+!   write(*,*)'-------------------------------------------------------------------------------'
+!   write(*,*)
+
+! end if
+
+end subroutine 
diff --git a/src/GW/GGW.f90 b/src/GW/GGW.f90
new file mode 100644
index 0000000..fed9179
--- /dev/null
+++ b/src/GW/GGW.f90
@@ -0,0 +1,119 @@
+subroutine GGW(doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,      &
+               exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA, &
+               linearize,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)
+
+! GW module
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  logical                       :: doG0W0
+  logical                       :: doevGW
+  logical                       :: doqsGW
+  logical                       :: doufG0W0
+  logical                       :: doufGW
+  logical                       :: doSRGqsGW
+
+  integer,intent(in)            :: maxSCF
+  integer,intent(in)            :: max_diis
+  double precision,intent(in)   :: thresh
+  logical,intent(in)            :: doACFDT
+  logical,intent(in)            :: exchange_kernel
+  logical,intent(in)            :: doXBS
+  logical,intent(in)            :: dophBSE
+  logical,intent(in)            :: dophBSE2
+  logical,intent(in)            :: TDA_W
+  logical,intent(in)            :: TDA
+  logical,intent(in)            :: dBSE
+  logical,intent(in)            :: dTDA
+  logical,intent(in)            :: doppBSE
+  logical,intent(in)            :: linearize
+  double precision,intent(in)   :: eta
+  logical,intent(in)            :: regularize
+
+  integer,intent(in)            :: nNuc
+  double precision,intent(in)   :: ZNuc(nNuc)
+  double precision,intent(in)   :: rNuc(nNuc,ncart)
+  double precision,intent(in)   :: ENuc
+
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC(nspin)
+  integer,intent(in)            :: nO(nspin)
+  integer,intent(in)            :: nV(nspin)
+  integer,intent(in)            :: nR(nspin)
+  integer,intent(in)            :: nS(nspin)
+
+  double precision,intent(in)   :: EHF
+  double precision,intent(in)   :: epsHF(nBas)
+  double precision,intent(in)   :: cHF(nBas,nBas)
+  double precision,intent(in)   :: PHF(nBas,nBas)
+  double precision,intent(in)   :: S(nBas,nBas)
+  double precision,intent(in)   :: T(nBas,nBas)
+  double precision,intent(in)   :: V(nBas,nBas)
+  double precision,intent(in)   :: Hc(nBas,nBas)
+  double precision,intent(in)   :: X(nBas,nBas)
+  double precision,intent(in)   :: ERI_AO(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: dipole_int_AO(nBas,nBas,ncart)
+  double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
+
+! Local variables
+
+  double precision              :: start_GW     ,end_GW       ,t_GW
+
+!------------------------------------------------------------------------
+! Perform G0W0 calculatiom
+!------------------------------------------------------------------------
+
+  if(doG0W0) then
+    
+    call wall_time(start_GW)
+    call GG0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
+              linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
+    call wall_time(end_GW)
+  
+    t_GW = end_GW - start_GW
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0W0 = ',t_GW,' seconds'
+    write(*,*)
+
+  end if
+
+!------------------------------------------------------------------------
+! Perform evGW calculation
+!------------------------------------------------------------------------
+
+  if(doevGW) then
+
+    call wall_time(start_GW)
+!   call evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
+!             linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
+    call wall_time(end_GW)
+
+    t_GW = end_GW - start_GW
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGW = ',t_GW,' seconds'
+    write(*,*)
+
+  end if
+
+!------------------------------------------------------------------------
+! Perform qsGW calculation
+!------------------------------------------------------------------------
+
+  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 wall_time(end_GW)
+
+    t_GW = end_GW - start_GW
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGW = ',t_GW,' seconds'
+    write(*,*)
+
+  end if
+
+end subroutine
diff --git a/src/GW/GGW_QP_graph.f90 b/src/GW/GGW_QP_graph.f90
new file mode 100644
index 0000000..231368e
--- /dev/null
+++ b/src/GW/GGW_QP_graph.f90
@@ -0,0 +1,82 @@
+subroutine GGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eOld,eGW,Z)
+
+! Compute the graphical solution of the QP equation
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+
+  double precision,intent(in)   :: eta
+  double precision,intent(in)   :: eHF(nBas)
+  double precision,intent(in)   :: Om(nS)
+  double precision,intent(in)   :: rho(nBas,nBas,nS)
+
+  double precision,intent(in)   :: eGWlin(nBas)
+  double precision,intent(in)   :: eOld(nBas)
+
+! Local variables
+
+  integer                       :: p
+  integer                       :: nIt
+  integer,parameter             :: maxIt = 64
+  double precision,parameter    :: thresh = 1d-6
+  double precision,external     :: GGW_SigC,GGW_dSigC
+  double precision              :: SigC,dSigC
+  double precision              :: f,df
+  double precision              :: w
+
+! Output variables
+
+  double precision,intent(out)  :: eGW(nBas)
+  double precision,intent(out)  :: Z(nBas)
+
+! Run Newton's algorithm to find the root
+
+  write(*,*)'-----------------------------------------------------'
+  write(*,'(A5,1X,A3,1X,A15,1X,A15,1X,A10)') 'Orb.','It.','e_GWlin (eV)','e_GW (eV)','Z'
+  write(*,*)'-----------------------------------------------------'
+
+  do p=nC+1,nBas-nR
+
+    w = eGWlin(p)
+    nIt = 0
+    f = 1d0
+    
+    do while (abs(f) > thresh .and. nIt < maxIt)
+    
+      nIt = nIt + 1
+
+      SigC  = GGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho)
+      dSigC = GGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho)
+      f  = w - eHF(p) - SigC
+      df = 1d0/(1d0 - dSigC)
+    
+      w = w - df*f
+    
+    end do
+ 
+    if(nIt == maxIt) then 
+
+      eGW(p) = eGWlin(p)
+      write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6,1X,A12)') p,nIt,eGWlin(p)*HaToeV,eGW(p)*HaToeV,Z(p),'Cvg Failed!'
+
+    else
+
+      eGW(p) = w
+      Z(p)   = df
+
+      write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6)') p,nIt,eGWlin(p)*HaToeV,eGW(p)*HaToeV,Z(p)
+
+    end if
+          
+  end do
+
+end subroutine 
diff --git a/src/GW/GGW_SigC.f90 b/src/GW/GGW_SigC.f90
new file mode 100644
index 0000000..6372537
--- /dev/null
+++ b/src/GW/GGW_SigC.f90
@@ -0,0 +1,52 @@
+double precision function GGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
+
+! Compute diagonal of the correlation part of the self-energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: p
+  double precision,intent(in)   :: w
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: e(nBas)
+  double precision,intent(in)   :: Om(nS)
+  double precision,intent(in)   :: rho(nBas,nBas,nS)
+
+! Local variables
+
+  integer                       :: i,a,m
+  double precision              :: num,eps
+
+! Initialize 
+
+  GGW_SigC = 0d0
+
+! Occupied part of the correlation self-energy
+
+  do i=nC+1,nO
+    do m=1,nS
+      eps = w - e(i) + Om(m)
+      num = rho(p,i,m)**2
+      GGW_SigC = GGW_SigC + num*eps/(eps**2 + eta**2)
+    end do
+  end do
+
+! Virtual part of the correlation self-energy
+
+  do a=nO+1,nBas-nR
+    do m=1,nS
+      eps = w - e(a) - Om(m)
+      num = rho(p,a,m)**2
+      GGW_SigC = GGW_SigC + num*eps/(eps**2 + eta**2)
+    end do
+  end do
+
+end function 
diff --git a/src/GW/GGW_dSigC.f90 b/src/GW/GGW_dSigC.f90
new file mode 100644
index 0000000..e1bd2c4
--- /dev/null
+++ b/src/GW/GGW_dSigC.f90
@@ -0,0 +1,52 @@
+double precision function GGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
+
+! Compute the derivative of the correlation part of the self-energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: p
+  double precision,intent(in)   :: w
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: e(nBas)
+  double precision,intent(in)   :: Om(nS)
+  double precision,intent(in)   :: rho(nBas,nBas,nS)
+
+! Local variables
+
+  integer                       :: i,a,m
+  double precision              :: num,eps
+
+! Initialize 
+
+  GGW_dSigC = 0d0
+
+! Occupied part of the correlation self-energy
+
+  do i=nC+1,nO
+    do m=1,nS
+      eps = w - e(i) + Om(m)
+      num = rho(p,i,m)**2
+      GGW_dSigC = GGW_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+    end do
+  end do
+
+! Virtual part of the correlation self-energy
+
+  do a=nO+1,nBas-nR
+    do m=1,nS
+      eps = w - e(a) - Om(m)
+      num = rho(p,a,m)**2
+      GGW_dSigC = GGW_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+    end do
+  end do
+
+end function 
diff --git a/src/GW/GGW_self_energy_diag.f90 b/src/GW/GGW_self_energy_diag.f90
new file mode 100644
index 0000000..2a02837
--- /dev/null
+++ b/src/GW/GGW_self_energy_diag.f90
@@ -0,0 +1,90 @@
+subroutine GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
+
+! Compute diagonal of the correlation part of the self-energy and the renormalization factor
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: e(nBas)
+  double precision,intent(in)   :: Om(nS)
+  double precision,intent(in)   :: rho(nBas,nBas,nS)
+
+! Local variables
+
+  integer                       :: i,a,p,m
+  double precision              :: num,eps
+
+! Output variables
+
+  double precision,intent(out)  :: Sig(nBas)
+  double precision,intent(out)  :: Z(nBas)
+  double precision,intent(out)  :: EcGM
+
+! Initialize 
+
+  Sig(:) = 0d0
+  Z(:)   = 0d0
+
+!----------------!
+! GW self-energy !
+!----------------!
+ 
+! Occupied part of the correlation self-energy
+
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+      do m=1,nS
+
+        eps = e(p) - e(i) + Om(m)
+        num = rho(p,i,m)**2
+        Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
+        Z(p)   = Z(p)   - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+
+      end do
+    end do
+  end do
+
+! Virtual part of the correlation self-energy
+
+  do p=nC+1,nBas-nR
+    do a=nO+1,nBas-nR
+      do m=1,nS
+
+        eps = e(p) - e(a) - Om(m)
+        num = rho(p,a,m)**2
+        Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
+        Z(p)   = Z(p)   - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+
+      end do
+    end do
+  end do
+
+! Galitskii-Migdal correlation energy
+
+  EcGM = 0d0
+  do i=nC+1,nO
+    do a=nO+1,nBas-nR
+      do m=1,nS
+
+        eps = e(a) - e(i) + Om(m)
+        num = rho(a,i,m)**2
+        EcGM = EcGM - num*eps/(eps**2 + eta**2)
+
+      end do
+    end do
+  end do
+
+! Compute renormalization factor from derivative 
+
+  Z(:) = 1d0/(1d0 - Z(:))
+
+end subroutine 
diff --git a/src/HF/GHF.f90 b/src/HF/GHF.f90
index 366b884..2e73046 100644
--- a/src/HF/GHF.f90
+++ b/src/HF/GHF.f90
@@ -56,6 +56,7 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
   double precision,allocatable  :: F(:,:)
   double precision,allocatable  :: Fp(:,:)
   double precision,allocatable  :: Cp(:,:)
+  double precision,allocatable  :: H(:,:)
   double precision,allocatable  :: S(:,:)
   double precision,allocatable  :: X(:,:)
   double precision,allocatable  :: err(:,:)
@@ -85,8 +86,8 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
   nBasSq = nBas*nBas
   nBas2Sq = nBas2*nBas2
 
-  nOa = nO(1)
-  nOb = nO(2)
+  nOa  = nO(1)
+  nOb  = nO(2)
   nOcc = nOa + nOb
 
 ! Memory allocation
@@ -96,16 +97,9 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
            Faa(nBas,nBas),Fab(nBas,nBas),Fba(nBas,nBas),Fbb(nBas,nBas), &
            Paa(nBas,nBas),Pab(nBas,nBas),Pba(nBas,nBas),Pbb(nBas,nBas), &
            F(nBas2,nBas2),Fp(nBas2,nBas2),Cp(nBas2,nBas2),              &
-           S(nBas2,nBas2),X(nBas2,nBas2),err(nBas2,nBas2),              &
+           H(nBas2,nBas2),S(nBas2,nBas2),X(nBas2,nBas2),err(nBas2,nBas2), &
            err_diis(nBas2Sq,max_diis),F_diis(nBas2Sq,max_diis))
 
-! Guess coefficients and demsity matrices
-
-! do ispin=1,nspin
-!   call mo_guess(nBas,guess_type,S,Hc,X,c(:,:,ispin))
-!   P(:,:,ispin) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin)))
-! end do
-
 ! Initialization
 
   nSCF = 0
@@ -127,6 +121,23 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
   X(     1:nBas ,     1:nBas ) = Or(1:nBas,1:nBas)
   X(nBas+1:nBas2,nBas+1:nBas2) = Or(1:nBas,1:nBas)
 
+! Construct super orthogonalization matrix
+
+  H(      :     ,       :    ) = 0d0
+  H(     1:nBas ,     1:nBas ) = Hc(1:nBas,1:nBas)
+  H(nBas+1:nBas2,nBas+1:nBas2) = Hc(1:nBas,1:nBas)
+
+! Guess coefficients and density matrices
+
+  call mo_guess(nBas2,guess_type,S,H,X,C)
+
+  P(:,:) = matmul(C(:,1:nOcc),transpose(C(:,1:nOcc)))
+
+  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)
+
 !------------------------------------------------------------------------
 ! Main SCF loop
 !------------------------------------------------------------------------
@@ -146,8 +157,6 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
 !   Build individual Hartree matrices
 
     call Hartree_matrix_AO_basis(nBas,Paa,ERI,Jaa)
-!   call Hartree_matrix_AO_basis(nBas,Pab,ERI,Jab)
-!   call Hartree_matrix_AO_basis(nBas,Pba,ERI,Jba)
     call Hartree_matrix_AO_basis(nBas,Pbb,ERI,Jbb)
 
 !   Compute individual exchange matrices
@@ -167,8 +176,8 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
 !  Build super Fock matrix
 
     F(     1:nBas ,     1:nBas ) = Faa(1:nBas,1:nBas)
-    F(nBas+1:nBas2,     1:nBas ) = Fab(1:nBas,1:nBas)
-    F(     1:nBas ,nBas+1:nBas2) = Fba(1:nBas,1:nBas)
+    F(     1:nBas ,nBas+1:nBas2) = Fab(1:nBas,1:nBas)
+    F(nBas+1:nBas2,     1:nBas ) = Fba(1:nBas,1:nBas)
     F(nBas+1:nBas2,nBas+1:nBas2) = Fbb(1:nBas,1:nBas)
 
 !   Check convergence 
@@ -217,11 +226,12 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
 !   if(nSCF == 1) call mix_guess(nBas,nO,mix,c)
 
 !   Compute individual density matrices
+    P(:,:) = matmul(C(:,1:nOcc),transpose(C(:,1:nOcc)))
 
-    Paa(:,:) = matmul(Ca(:,1:nOcc),transpose(Ca(:,1:nOcc)))
-    Pab(:,:) = matmul(Ca(:,1:nOcc),transpose(Cb(:,1:nOcc)))
-    Pba(:,:) = matmul(Cb(:,1:nOcc),transpose(Ca(:,1:nOcc)))
-    Pbb(:,:) = matmul(Cb(:,1:nOcc),transpose(Cb(:,1:nOcc)))
+    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)
  
 !------------------------------------------------------------------------
 !   Compute UHF energy
@@ -243,19 +253,19 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
 
 !  Hartree energy: 16 terms?
 
-    EJaaaa = +0.5d0*trace_matrix(nBas,matmul(Paa,Jaa))
-    EJaabb = +0.5d0*trace_matrix(nBas,matmul(Paa,Jbb))
-    EJbbaa = +0.5d0*trace_matrix(nBas,matmul(Pbb,Jaa))
-    EJbbbb = +0.5d0*trace_matrix(nBas,matmul(Pbb,Jbb))
+    EJaaaa = 0.5d0*trace_matrix(nBas,matmul(Paa,Jaa))
+    EJaabb = 0.5d0*trace_matrix(nBas,matmul(Paa,Jbb))
+    EJbbaa = 0.5d0*trace_matrix(nBas,matmul(Pbb,Jaa))
+    EJbbbb = 0.5d0*trace_matrix(nBas,matmul(Pbb,Jbb))
 
     EJ = EJaaaa + EJaabb + EJbbaa + EJbbbb
   
 !   Exchange energy
 
-    Exaaaa = -0.5d0*trace_matrix(nBas,matmul(Paa,Kaa))
-    Exabba = -0.5d0*trace_matrix(nBas,matmul(Pab,Kba))
-    Exbaab = -0.5d0*trace_matrix(nBas,matmul(Pba,Kab))
-    Exbbbb = -0.5d0*trace_matrix(nBas,matmul(Pbb,Kbb))
+    Exaaaa = 0.5d0*trace_matrix(nBas,matmul(Paa,Kaa))
+    Exabba = 0.5d0*trace_matrix(nBas,matmul(Pab,Kba))
+    Exbaab = 0.5d0*trace_matrix(nBas,matmul(Pba,Kab))
+    Exbbbb = 0.5d0*trace_matrix(nBas,matmul(Pbb,Kbb))
 
     Ex = Exaaaa + Exabba + Exbaab + Exbbbb
 
diff --git a/src/HF/GHF_stability.f90 b/src/HF/GHF_stability.f90
new file mode 100644
index 0000000..718e9b5
--- /dev/null
+++ b/src/HF/GHF_stability.f90
@@ -0,0 +1,109 @@
+subroutine GHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI)
+
+! Perform a stability analysis of the GHF solution
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+
+  double precision,intent(in)   :: eHF(nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+
+! Local variables
+
+  integer,parameter             :: maxS = 20
+  integer                       :: ia
+  integer                       :: ispin
+
+  double precision,allocatable  :: A(:,:)
+  double precision,allocatable  :: B(:,:)
+  double precision,allocatable  :: AB(:,:)
+  double precision,allocatable  :: Om(:)
+
+! Memory allocation
+
+  allocate(A(nS,nS),B(nS,nS),AB(nS,nS),Om(nS))
+ 
+!-------------------------------------------------------------!
+! Stability analysis: Real GHF -> Real GHF  
+!-------------------------------------------------------------!
+ 
+  ispin = 3
+
+  call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A)
+  call phLR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B)
+
+  AB(:,:) = A(:,:) + B(:,:)
+
+  call diagonalize_matrix(nS,AB,Om)
+
+  write(*,*)'-------------------------------------------------------------'
+  write(*,*)'|       Stability analysis: Real GHF -> Real GHF            |'
+  write(*,*)'-------------------------------------------------------------'
+  write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') &
+            '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|'
+  write(*,*)'-------------------------------------------------------------'
+  do ia=1,min(nS,maxS)
+    write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') &
+      '|',ia,'|',Om(ia),'|',Om(ia)*HaToeV,'|'
+  enddo
+  write(*,*)'-------------------------------------------------------------'
+
+  if(minval(Om(:)) < 0d0) then
+
+    write(*,'(1X,A40,1X)')        'Too bad, GHF solution is unstable!'
+    write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue: ',Om(1),' au'
+
+  else
+
+    write(*,'(1X,A40,1X)')        'Well done, GHF solution is stable!'
+    write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om(1),' au'
+
+  end if
+  write(*,*)'-------------------------------------------------------------'
+  write(*,*)
+
+!-------------------------------------------------------------!
+! Stability analysis: Real GHF -> Complex GHF  
+!-------------------------------------------------------------!
+
+  AB(:,:) = A(:,:) - B(:,:)
+
+  call diagonalize_matrix(nS,AB,Om)
+
+  write(*,*)'-------------------------------------------------------------'
+  write(*,*)'|       Stability analysis: Real GHF -> Complex GHF         |'
+  write(*,*)'-------------------------------------------------------------'
+  write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') &
+            '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|'
+  write(*,*)'-------------------------------------------------------------'
+  do ia=1,min(nS,maxS)
+    write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') &
+      '|',ia,'|',Om(ia),'|',Om(ia)*HaToeV,'|'
+  enddo
+  write(*,*)'-------------------------------------------------------------'
+
+  if(minval(Om(:)) < 0d0) then
+
+    write(*,'(1X,A40,1X)')        'Too bad, GHF solution is unstable!'
+    write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue: ',Om(1),' au'
+
+  else
+
+    write(*,'(1X,A40,1X)')        'Well done, GHF solution is stable!'
+    write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om(1),' au'
+
+  end if
+  write(*,*)'-------------------------------------------------------------'
+  write(*,*)
+
+    
+end subroutine 
diff --git a/src/HF/mo_guess.f90 b/src/HF/mo_guess.f90
index 9822334..91e6e88 100644
--- a/src/HF/mo_guess.f90
+++ b/src/HF/mo_guess.f90
@@ -28,6 +28,8 @@ subroutine mo_guess(nBas,guess_type,S,Hc,X,c)
 
     call random_number(c)
 
+    c(:,:) = 2d0*c(:,:) - 1d0
+
   else
 
     print*,'Wrong guess option'
diff --git a/src/MP/GMP.f90 b/src/MP/GMP.f90
new file mode 100644
index 0000000..2fe338d
--- /dev/null
+++ b/src/MP/GMP.f90
@@ -0,0 +1,45 @@
+subroutine GMP(doMP2,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
+
+! Moller-Plesset module
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  logical,intent(in)            :: doMP2
+
+  logical,intent(in)            :: regularize
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  double precision,intent(in)   :: ENuc
+  double precision,intent(in)   :: EHF
+  double precision,intent(in)   :: epsHF(nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+
+! Local variables
+
+  double precision              :: start_MP     ,end_MP       ,t_MP
+
+! Output variables
+
+!------------------------------------------------------------------------
+! Compute MP2 energy
+!------------------------------------------------------------------------                               
+
+  if(doMP2) then    
+       
+    call wall_time(start_MP)
+    call GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
+    call wall_time(end_MP)
+
+    t_MP = end_MP - start_MP
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2 = ',t_MP,' seconds'
+    write(*,*)
+                    
+  end if
+
+end subroutine 
diff --git a/src/MP/GMP2.f90 b/src/MP/GMP2.f90
new file mode 100644
index 0000000..0320160
--- /dev/null
+++ b/src/MP/GMP2.f90
@@ -0,0 +1,174 @@
+subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e)
+
+! Perform second-order Moller-Plesset calculation with and without regularizers
+
+  implicit none
+
+! Input variables
+
+  logical,intent(in)            :: regularize
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  double precision,intent(in)   :: ENuc
+  double precision,intent(in)   :: EHF
+  double precision,intent(in)   :: e(nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+
+! Local variables
+
+  integer                       :: i,j,a,b
+
+  double precision              :: kappa,sigm1,sigm2
+  double precision              :: Dijab
+  double precision              :: fs,fs2,fk
+
+  double precision              :: E2d,E2ds,E2ds2,E2dk
+  double precision              :: E2x,E2xs,E2xs2,E2xk
+
+  double precision              :: EcsMP2,Ecs2MP2,EckMP2
+  double precision              :: EcMP2
+
+! Output variables
+
+! Hello world
+
+  write(*,*)
+  write(*,*)'************************************************'
+  write(*,*)'|    Moller-Plesset second-order calculation   |'
+  write(*,*)'************************************************'
+  write(*,*)
+
+!---------------------------------------------!
+! Parameters for regularized MP2 calculations !
+!---------------------------------------------!
+
+  kappa = 1.1d0
+  sigm1 = 0.7d0
+  sigm2 = 0.4d0
+
+!--------------------------------------------------!
+! Compute conventinal and regularized MP2 energies !
+!--------------------------------------------------!
+
+  E2d   = 0d0
+  E2ds  = 0d0
+  E2ds2 = 0d0
+  E2dk  = 0d0
+
+  E2x   = 0d0
+  E2xs  = 0d0
+  E2xs2 = 0d0
+  E2xk  = 0d0
+
+  do i=nC+1,nO
+    do j=nC+1,nO
+      do a=nO+1,nBas-nR
+        do b=nO+1,nBas-nR
+
+          Dijab = e(a) + e(b) - e(i) - e(j) 
+
+          ! Second-order ring diagram
+
+          fs  = (1d0 - exp(-sigm1*Dijab))/Dijab
+          fs2 = (1d0 - exp(-sigm2*Dijab*Dijab))/Dijab
+          fk  = (1d0 - exp(-kappa*Dijab))**2/Dijab
+
+          E2d   = E2d   - 0.25d0*(ERI(i,j,a,b)-ERI(i,j,b,a))**2/Dijab
+          E2ds  = E2ds  - ERI(i,j,a,b)*ERI(i,j,a,b)*fs 
+          E2ds2 = E2ds2 - ERI(i,j,a,b)*ERI(i,j,a,b)*fs2
+          E2dk  = E2dk  - ERI(i,j,a,b)*ERI(i,j,a,b)*fk
+
+          ! Second-order exchange diagram
+
+          E2x   = E2x   - ERI(i,j,a,b)*ERI(i,j,b,a)/Dijab
+          E2x   = 0d0
+          E2xs  = E2xs  - ERI(i,j,a,b)*ERI(i,j,b,a)*fs 
+          E2xs2 = E2xs2 - ERI(i,j,a,b)*ERI(i,j,b,a)*fs2
+          E2xk  = E2xk  - ERI(i,j,a,b)*ERI(i,j,b,a)*fk
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  EcMP2   = E2d   - E2x
+  EcsMP2  = E2ds  - E2xs
+  Ecs2MP2 = E2ds2 - E2xs2
+  EckMP2  = E2dk  - E2xk
+
+!------------!
+! MP2 energy !
+!------------!
+
+  write(*,*)
+  write(*,'(A32)')           '---------------------------'
+  write(*,'(A32)')           ' GMP2 calculation          '
+  write(*,'(A32)')           '---------------------------'
+  write(*,'(A32,1X,F16.10)') ' GMP2 correlation energy = ',EcMP2
+  write(*,'(A32,1X,F16.10)') ' Direct part             = ',E2d
+  write(*,'(A32,1X,F16.10)') ' Exchange part           = ',-E2x
+  write(*,'(A32)')           '---------------------------'
+  write(*,'(A32,1X,F16.10)') ' GMP2 electronic  energy = ',EHF + EcMP2
+  write(*,'(A32,1X,F16.10)') ' GMP2 total       energy = ',ENuc + EHF + EcMP2
+  write(*,'(A32)')           '---------------------------'
+  write(*,*)
+
+  if(regularize) then
+
+!-------------------!
+! sigma1-MP2 energy !
+!-------------------!
+
+    write(*,*)
+    write(*,'(A32)')           '---------------------------'
+    write(*,'(A32)')           ' sigma-GMP2 calculation    '
+    write(*,'(A32)')           '---------------------------'
+    write(*,'(A32,1X,F16.10)') ' GMP2 correlation energy = ',EcsMP2
+    write(*,'(A32,1X,F16.10)') ' Direct part             = ',E2ds
+    write(*,'(A32,1X,F16.10)') ' Exchange part           = ',-E2xs
+    write(*,'(A32)')           '---------------------------'
+    write(*,'(A32,1X,F16.10)') ' GMP2 electronic  energy = ',EHF + EcsMP2
+    write(*,'(A32,1X,F16.10)') ' GMP2 total       energy = ',ENuc + EHF + EcsMP2
+    write(*,'(A32)')           '---------------------------'
+    write(*,*)
+
+!--------------------!
+! sigma^2-MP2 energy !
+!--------------------!
+
+    write(*,*)
+    write(*,'(A32)')           '---------------------------'
+    write(*,'(A32)')           ' sigma^2-GMP2 calculation  '
+    write(*,'(A32)')           '---------------------------'
+    write(*,'(A32,1X,F16.10)') ' GMP2 correlation energy = ',Ecs2MP2
+    write(*,'(A32,1X,F16.10)') ' Direct part             = ',E2ds2
+    write(*,'(A32,1X,F16.10)') ' Exchange part           = ',-E2xs2
+    write(*,'(A32)')           '---------------------------'
+    write(*,'(A32,1X,F16.10)') ' GMP2 electronic  energy = ',EHF + Ecs2MP2
+    write(*,'(A32,1X,F16.10)') ' GMP2 total       energy = ',ENuc + EHF + Ecs2MP2
+    write(*,'(A32)')           '---------------------------'
+    write(*,*)
+
+!------------------!
+! kappa-MP2 energy !
+!------------------!
+
+    write(*,*)
+    write(*,'(A32)')           '---------------------------'
+    write(*,'(A32)')           ' kappa-GMP2 calculation    '
+    write(*,'(A32)')           '---------------------------'
+    write(*,'(A32,1X,F16.10)') ' GMP2 correlation energy = ',EckMP2
+    write(*,'(A32,1X,F16.10)') ' Direct part             = ',E2dk
+    write(*,'(A32,1X,F16.10)') ' Exchange part           = ',-E2xk
+    write(*,'(A32)')           '---------------------------'
+    write(*,'(A32,1X,F16.10)') ' GMP2 electronic  energy = ',EHF + EckMP2
+    write(*,'(A32,1X,F16.10)') ' GMP2 total       energy = ',ENuc + EHF + EckMP2
+    write(*,'(A32)')           '---------------------------'
+    write(*,*)
+
+  end if
+
+end subroutine 
diff --git a/src/QuAcK/GQuAcK.f90 b/src/QuAcK/GQuAcK.f90
index 754a9fd..3c207f4 100644
--- a/src/QuAcK/GQuAcK.f90
+++ b/src/QuAcK/GQuAcK.f90
@@ -72,6 +72,8 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
   double precision              :: EHF
   double precision,allocatable  :: dipole_int_MO(:,:,:)
   double precision,allocatable  :: ERI_MO(:,:,:,:)
+  double precision,allocatable  :: ERI_tmp(:,:,:,:)
+  double precision,allocatable  :: Ca(:,:),Cb(:,:)
 
   integer                       :: ixyz
   integer                       :: nBas2
@@ -128,7 +130,17 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
   
   ! 4-index transform 
   
-  call AOtoMO_integral_transform(1,1,1,1,nBas2,cHF,ERI_AO,ERI_MO)
+  allocate(Ca(nBas,nBas2),Cb(nBas,nBas2),ERI_tmp(nBas2,nBas2,nBas2,nBas2))
+  Ca(:,:) = cHF(1:nBas,1:nBas2)
+  Cb(:,:) = cHF(nBas+1:nBas2,1:nBas2)
+  call AOtoMO_integral_transform_GHF(nBas,nBas2,Ca,Ca,Ca,Ca,ERI_AO,ERI_tmp)
+  ERI_MO(:,:,:,:) = ERI_tmp(:,:,:,:)
+  call AOtoMO_integral_transform_GHF(nBas,nBas2,Ca,Cb,Ca,Cb,ERI_AO,ERI_tmp)
+  ERI_MO(:,:,:,:) = ERI_MO(:,:,:,:) + ERI_tmp(:,:,:,:)
+  call AOtoMO_integral_transform_GHF(nBas,nBas2,Cb,Ca,Cb,Ca,ERI_AO,ERI_tmp)
+  ERI_MO(:,:,:,:) = ERI_MO(:,:,:,:) + ERI_tmp(:,:,:,:)
+  call AOtoMO_integral_transform_GHF(nBas,nBas2,Cb,Cb,Cb,Cb,ERI_AO,ERI_tmp)
+  ERI_MO(:,:,:,:) = ERI_MO(:,:,:,:) + ERI_tmp(:,:,:,:)
 
   call wall_time(end_AOtoMO)
 
@@ -143,7 +155,7 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
   if(dostab) then
 
     call wall_time(start_stab)
-!   call GHF_stability(nBas,nC,nO,nV,nR,nS,epsHF,ERI_MO)
+    call GHF_stability(nBas2,sum(nC),sum(nO),sum(nV),sum(nR),sum(nO)*sum(nV),epsHF,ERI_MO)
     call wall_time(end_stab)
 
     t_stab = end_stab - start_stab
@@ -161,7 +173,7 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
   if(doMP) then
 
     call wall_time(start_MP)
-!   call GMP(doMP2,reg_MP,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
+    call GMP(doMP2,reg_MP,nBas2,sum(nC),sum(nO),sum(nV),sum(nR),ERI_MO,ENuc,EHF,epsHF)
     call wall_time(end_MP)
 
     t_MP = end_MP - start_MP
@@ -179,8 +191,8 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
   if(doRPA) then
 
     call wall_time(start_RPA)
-!   call GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip, & 
-!             nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF,cHF,S)
+    call GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel, & 
+              nBas2,sum(nC),sum(nO),sum(nV),sum(nR),sum(nO)*sum(nV),ENuc,EHF,ERI_MO,dipole_int_MO,epsHF,cHF,S)
     call wall_time(end_RPA)
 
     t_RPA = end_RPA - start_RPA
@@ -219,10 +231,10 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
   if(doGW) then
     
     call wall_time(start_GW)
-!   call GGW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT, &
-!            exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,                     &
-!            lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,            &  
-!            ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
+    call GGW(doG0W0,doevGW,doqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT, &
+             exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,                     &
+             lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas2,sum(nC),sum(nO),sum(nV),sum(nR),sum(nO)*sum(nV), &
+             EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
     call wall_time(end_GW)
   
     t_GW = end_GW - start_GW
diff --git a/src/RPA/GRPA.f90 b/src/RPA/GRPA.f90
new file mode 100644
index 0000000..9aae7f6
--- /dev/null
+++ b/src/RPA/GRPA.f90
@@ -0,0 +1,82 @@
+subroutine GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
+
+! Random-phase approximation module
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  logical                       :: dophRPA
+  logical                       :: dophRPAx
+  logical                       :: doppRPA
+
+  logical,intent(in)            :: TDA
+  logical,intent(in)            :: doACFDT
+  logical,intent(in)            :: exchange_kernel
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC(nspin)
+  integer,intent(in)            :: nO(nspin)
+  integer,intent(in)            :: nV(nspin)
+  integer,intent(in)            :: nR(nspin)
+  integer,intent(in)            :: nS(nspin)
+  double precision,intent(in)   :: ENuc
+  double precision,intent(in)   :: EHF
+  double precision,intent(in)   :: epsHF(nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
+
+! Local variables
+
+  double precision              :: start_RPA    ,end_RPA      ,t_RPA
+
+!------------------------------------------------------------------------
+! Compute (direct) RPA excitations
+!------------------------------------------------------------------------
+
+  if(dophRPA) then
+
+    call wall_time(start_RPA)
+    call phGRPA(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
+    call wall_time(end_RPA)
+
+    t_RPA = end_RPA - start_RPA
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPA = ',t_RPA,' seconds'
+    write(*,*)
+
+  end if
+
+!------------------------------------------------------------------------
+! Compute RPAx (RPA with exchange) excitations
+!------------------------------------------------------------------------
+
+  if(dophRPAx) then
+
+    call wall_time(start_RPA)
+    call phGRPAx(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
+
+    call wall_time(end_RPA)
+
+    t_RPA = end_RPA - start_RPA
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPAx = ',t_RPA,' seconds'
+    write(*,*)
+
+  end if
+
+!------------------------------------------------------------------------
+! Compute ppRPA excitations
+!------------------------------------------------------------------------
+
+  if(doppRPA) then
+
+    call wall_time(start_RPA)
+!   call ppGRPA(TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,epsHF)
+    call wall_time(end_RPA)
+
+    t_RPA = end_RPA - start_RPA
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds'
+    write(*,*)
+
+  end if
+
+end subroutine
diff --git a/src/RPA/URPA.f90 b/src/RPA/URPA.f90
index da77aa8..a312009 100644
--- a/src/RPA/URPA.f90
+++ b/src/RPA/URPA.f90
@@ -46,8 +46,8 @@ subroutine URPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,spi
   if(dophRPA) then
 
     call wall_time(start_RPA)
-     call phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EHF, &
-                 ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
+    call phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EHF, &
+                ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
     call wall_time(end_RPA)
 
     t_RPA = end_RPA - start_RPA
@@ -63,8 +63,8 @@ subroutine URPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,spi
   if(dophRPAx) then
 
     call wall_time(start_RPA)
-     call phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EHF, &
-                  ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
+    call phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EHF, &
+                 ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
     call wall_time(end_RPA)
 
     t_RPA = end_RPA - start_RPA
diff --git a/src/RPA/phGRPA.f90 b/src/RPA/phGRPA.f90
new file mode 100644
index 0000000..3ae3e6e
--- /dev/null
+++ b/src/RPA/phGRPA.f90
@@ -0,0 +1,78 @@
+subroutine phGRPA(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
+
+! Perform a direct random phase approximation calculation
+
+  implicit none
+  include 'parameters.h'
+  include 'quadrature.h'
+
+! Input variables
+
+  logical,intent(in)            :: TDA
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: ENuc
+  double precision,intent(in)   :: EHF
+  double precision,intent(in)   :: e(nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
+
+! Local variables
+
+  integer                       :: ispin
+  logical                       :: dRPA
+  double precision,allocatable  :: Aph(:,:)
+  double precision,allocatable  :: Bph(:,:)
+  double precision,allocatable  :: Om(:)
+  double precision,allocatable  :: XpY(:,:)
+  double precision,allocatable  :: XmY(:,:)
+
+  double precision              :: EcRPA
+
+! Hello world
+
+  write(*,*)
+  write(*,*)'***********************************************'
+  write(*,*)'|  Random-phase approximation calculation     |'
+  write(*,*)'***********************************************'
+  write(*,*)
+
+! TDA 
+
+  if(TDA) then
+    write(*,*) 'Tamm-Dancoff approximation activated!'
+    write(*,*)
+  end if
+
+! Initialization
+
+  dRPA = .true.
+
+  EcRPA = 0d0
+
+! Memory allocation
+
+  allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS))
+  if(.not.TDA) allocate(Bph(nS,nS))
+
+  ispin = 3
+
+  call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph)
+  if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+
+  call phLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+  call print_excitation_energies('phRPA@GHF',ispin,nS,Om)
+  call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
+
+  write(*,*)
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy           =',EcRPA
+  write(*,'(2X,A50,F20.10)') 'Tr@phRPA total energy                 =',ENuc + EHF + EcRPA
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,*)
+
+end subroutine 
diff --git a/src/RPA/phGRPAx.f90 b/src/RPA/phGRPAx.f90
new file mode 100644
index 0000000..7fc1a7a
--- /dev/null
+++ b/src/RPA/phGRPAx.f90
@@ -0,0 +1,79 @@
+subroutine phGRPAx(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
+
+! Perform random phase approximation calculation with exchange (aka TDHF)
+
+  implicit none
+  include 'parameters.h'
+  include 'quadrature.h'
+
+! Input variables
+
+  logical,intent(in)            :: TDA
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: ENuc
+  double precision,intent(in)   :: EHF
+  double precision,intent(in)   :: e(nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
+
+! Local variables
+
+  integer                       :: ispin
+  logical                       :: dRPA
+  double precision,allocatable  :: Aph(:,:)
+  double precision,allocatable  :: Bph(:,:)
+  double precision,allocatable  :: Om(:)
+  double precision,allocatable  :: XpY(:,:)
+  double precision,allocatable  :: XmY(:,:)
+
+  double precision              :: EcRPA
+
+! Hello world
+
+  write(*,*)
+  write(*,*)'***********************************************************'
+  write(*,*)'|  Random phase approximation calculation with exchange   |'
+  write(*,*)'***********************************************************'
+  write(*,*)
+
+! TDA 
+
+  if(TDA) then
+    write(*,*) 'Tamm-Dancoff approximation activated!'
+    write(*,*) ' => RPAx + TDA = CIS '
+    write(*,*)
+  end if
+
+! Initialization
+
+  dRPA = .false.
+
+  EcRPA = 0d0
+
+! Memory allocation
+
+  allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS))
+  if(.not.TDA) allocate(Bph(nS,nS))
+
+  ispin = 3
+
+  call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph)
+  if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+
+  call phLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+  call print_excitation_energies('phRPAx@GHF',ispin,nS,Om)
+  call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
+
+  write(*,*)
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy           =',EcRPA
+  write(*,'(2X,A50,F20.10)') 'Tr@phRPAx total energy                 =',ENuc + EHF + EcRPA
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,*)
+
+end subroutine