diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90
index ff40f9c..fd1a3e7 100644
--- a/src/GW/G0W0.f90
+++ b/src/GW/G0W0.f90
@@ -159,10 +159,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
   
     call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Om,rho,eGWlin,eGW,regularize)
 
-    ! Find all the roots of the QP equation if necessary
-
-    ! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Om,rho,eGWlin) 
- 
   end if
 
 ! Compute the RPA correlation energy
diff --git a/src/GW/QP_roots.f90 b/src/GW/QP_roots.f90
deleted file mode 100644
index 99a848a..0000000
--- a/src/GW/QP_roots.f90
+++ /dev/null
@@ -1,182 +0,0 @@
-subroutine QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGW)
-
-! Compute all the roots of the QP equation for each orbital
-
-  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)   :: Omega(nS)
-  double precision,intent(in)   :: rho(nBas,nBas,nS)
-
-! Local variables
-
-  integer                       :: i,j,a,b,x,jb,g
-  integer,parameter             :: nGrid = 1000000
-  double precision,parameter    :: wmin = -50d0
-  double precision,parameter    :: wmax = +50d0
-  double precision              :: dw
-  double precision              :: eps
-  double precision              :: left_sign,right_sign
-  double precision              :: left_QP,right_QP
-  double precision,allocatable  :: w(:)
-  double precision,allocatable  :: SigC(:)
-  double precision,allocatable  :: dSigC(:)
-
-  integer                       :: nIt
-  integer,parameter             :: maxIt = 64
-  double precision,parameter    :: thresh = 1d-6
-  double precision,external     :: SigmaC,dSigmaC
-  double precision              :: sig,dsig
-  double precision              :: sat,dsat
-  double precision              :: om
-  integer                       :: nRoot
-  double precision              :: Z
-  double precision              :: Ztot
-  double precision              :: QP
-
-! Output variables
-
-  double precision,intent(out)  :: eGW(nBas)
-
-! Construct grid
-
-  allocate(w(nGrid),SigC(nGrid),dSigC(nGrid))
-
-! Minimum and maximum frequency values
-
-  dw = (wmax - wmin)/dble(nGrid)
-
-  do g=1,nGrid
-    w(g) = wmin + dble(g)*dw
-  enddo
-
-! Main loop over the orbitals
-
-  do x=nC+1,nBas-nR
-
-  SigC(:)  = 0d0
-  dSigC(:) = 0d0
-
-    ! Loop over grid points
- 
-    do g=1,nGrid
-
-      ! Occupied part of the self-energy and spectral weight
-
-      do i=nC+1,nO
-        jb = 0
-        do j=nC+1,nO
-          do b=nO+1,nBas-nR
-
-            jb = jb + 1
-            eps = w(g) - eHF(i) + Omega(jb)
-            SigC(g)  = SigC(g)  + 2d0*rho(x,i,jb)**2*eps/(eps**2 + eta**2)
-            dSigC(g) = dSigC(g) - 2d0*rho(x,i,jb)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
-
-          enddo
-        enddo
-      enddo
- 
-      ! Virtual part of the self-energy and spectral weight
- 
-      do a=nO+1,nBas-nR
-        jb = 0
-        do j=nC+1,nO
-          do b=nO+1,nBas-nR
-
-            jb = jb + 1
-            eps = w(g) - eHF(a) - Omega(jb)
-            SigC(g)  = SigC(g)  + 2d0*rho(x,a,jb)**2*eps/(eps**2 + eta**2)
-            dSigC(g) = dSigC(g) - 2d0*rho(x,a,jb)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2
-
-          enddo                                               
-        enddo                                                 
-      enddo                                                   
-
-    enddo
-
-! Find the zeros of the QP equation
-
-    nRoot = 0
-    QP    = 0d0
-    Ztot  = 0d0
-
-    write(*,*) '-----------------'
-    write(*,'(A10,I3)') 'Orbital ',x
-    write(*,*) '-----------------'
-    write(*,*) 
-    left_QP   = w(1) - eHF(x) - SigC(1)
-    left_sign = sign(1d0,left_QP)
-
-    do g=2,nGrid
-
-      right_QP = w(g) - eHF(x) - SigC(g)
-      right_sign = sign(1d0,right_QP)
-
-      if(left_sign /= right_sign .and. left_sign == sign(1d0,-1d0)) then
-
-        nRoot = nRoot + 1
-        write(*,'(A20,I6,F10.6,F10.6,I6,F10.6,F10.6,F10.6)') & 
-          'root right here!',g-1,w(g-1),left_QP,g,w(g),right_QP,1d0/(1d0-dSigC(g))
-
-        ! Run Newton's algorithm to find the root
-
-          om = w(g-1)
-          nIt = 0
-          sat = 1d0
-    
-          do while (abs(sat) > thresh .and. nIt < maxIt)
-    
-            nIt = nIt + 1
-
-            sig  =  SigmaC(x,om,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho)
-            dsig = dSigmaC(x,om,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho)
-            sat  = om - eHF(x) - sig
-            dsat = 1d0 - dsig
-            Z    = 1d0/dsat
-    
-            write(*,'(A3,I3,A1,1X,4F15.9)') 'It.',nIt,':',om,sat,sig,Z
-    
-            om = om - sat/dsat
-    
-          end do
- 
-          if(nIt == maxIt) then 
-
-            write(*,*) 'Newton root search has not converged!'
-
-          else 
-
-            Ztot = Ztot + Z
-            QP   = QP + Z*om
-
-          end if
-          
-      end if
-
-      left_sign = right_sign
-      left_QP   = right_QP
-   
-    end do
-
-    eGW(x) = QP/Ztot
-
-    write(*,*)
-    write(*,'(A32,I3,A1,I3)') 'Number of roots for orbital ',x,':',nRoot
-    write(*,'(A32,F16.10)')   'Total weight                ',Ztot
-    write(*,'(A32,F16.10)')   'Quasiparticle energy (eV)   ',eGW(x)*HaToeV
-    write(*,*)
-
-  end do
-
-end subroutine QP_roots