subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)

! CC-based GW module

  implicit none
  include 'parameters.h'

! Input variables

  integer,intent(in)            :: maxSCF
  double precision,intent(in)   :: thresh
  integer,intent(in)            :: max_diis

  integer,intent(in)            :: nBas
  integer,intent(in)            :: nOrb
  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)   :: ERHF
  double precision,intent(in)   :: eHF(nOrb)
  double precision,intent(in)   :: ERI(nOrb,nOrb,nOrb,nOrb)

! Local variables

  integer                       :: p,q,r,s
  integer                       :: i,j,k,l
  integer                       :: a,b,c,d

  integer                       :: nSCF
  double precision              :: Conv

  double precision              :: x

  double precision,allocatable  :: eGW(:)
  double precision,allocatable  :: Z(:)

  double precision,allocatable  :: del(:,:,:)
  double precision,allocatable  :: vec(:,:,:)
  double precision,allocatable  :: res(:,:,:)
  double precision,allocatable  :: amp(:,:,:)

  integer                       :: n_diis
  double precision              :: rcond
  double precision,allocatable  :: r_diis(:,:)
  double precision,allocatable  :: t_diis(:,:)


! Hello world

  write(*,*)
  write(*,*)'*****************************'
  write(*,*)'* CC-based G0W0 Calculation *'
  write(*,*)'*****************************'
  write(*,*)

! Form energy denominator and guess amplitudes

  allocate(del(nO,nV,nOrb))
  allocate(vec(nO,nV,nOrb))
  allocate(res(nO,nV,nOrb))
  allocate(amp(nO,nV,nOrb))

  allocate(eGW(nOrb),Z(nOrb))

 allocate(r_diis(nO*nV*nOrb,max_diis))
 allocate(t_diis(nO*nV*nOrb,max_diis))

! Initialization

  eGW(:) = eHF(:)
  Z(:)   = 1d0

  ! Compute energy differences
 
  do i=nC+1,nO
    do j=nC+1,nO
      do a=1,nV-nR
  
        del(i,a,j) = eHF(i) + eHF(j) - eHF(nO+a)
  
      end do
    end do
  end do
  
  do i=nC+1,nO
    do a=1,nV-nR
      do b=1,nV-nR
  
        del(i,a,nO+b) = eHF(nO+a) + eHF(nO+b) - eHF(i)
  
      end do
    end do
  end do

!-------------------------!
! Main loop over orbitals !
!-------------------------!

  do p=nC+1,nO+1

    ! Initialization
 
    Conv = 1d0
    nSCF =  0

    n_diis      = 0
    t_diis(:,:) = 0d0
    r_diis(:,:) = 0d0
    rcond       = 0d0

    amp(:,:,:) = 0d0
    res(:,:,:) = 0d0
 
  ! Compute energy differences

    do i=nC+1,nO
      do a=1,nV-nR
        do j=nC+1,nO
 
          vec(i,a,j) = sqrt(2d0)*ERI(p,nO+a,i,j)
 
        end do
      end do
    end do
 
    do i=nC+1,nO
      do a=1,nV-nR
        do b=1,nV-nR
 
          vec(i,a,nO+b) = sqrt(2d0)*ERI(p,i,nO+b,nO+a)
 
        end do
      end do
    end do

   !----------------------!
   ! Loop over amplitudes !
   !----------------------!

    write(*,*)
    write(*,*)'-------------------------------------------------------------'
    write(*,*)'| CC-based G0W0 calculation                                 |'
    write(*,*)'-------------------------------------------------------------'
    write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
              '|','#','|','HF','|','G0W0','|','Conv','|'
    write(*,*)'-------------------------------------------------------------'
 
    do while(Conv > thresh .and. nSCF < maxSCF)
 
     ! Increment 
 
      nSCF = nSCF + 1
 
      ! Compute residual for 2h1p sector
 
      res(:,:,:) = vec(:,:,:) + (del(:,:,:) - eGW(p))*amp(:,:,:)
 
      do i=nC+1,nO
        do a=1,nV-nR
          do j=nC+1,nO
 
            do k=nC+1,nO
              do c=1,nV-nR
 
                res(i,a,j) = res(i,a,j) - 2d0*ERI(j,nO+c,nO+a,k)*amp(i,c,k)
 
              end do
            end do
 
          end do
        end do
      end do
 
      ! Compute residual for 2p1h sector
 
      do i=nC+1,nO
        do a=1,nV-nR
          do b=1,nV-nR
 
            do k=nC+1,nO
              do c=1,nV-nR
 
                res(i,a,nO+b) = res(i,a,nO+b) + 2d0*ERI(nO+a,k,i,nO+c)*amp(k,c,nO+b) 
 
              end do
            end do
 
          end do
        end do
      end do
  
      ! Check convergence 
 
      Conv = maxval(abs(res))
    
      ! Update amplitudes

      amp(:,:,:) = amp(:,:,:) - res(:,:,:)/(del(:,:,:) - eHF(p))
 
      ! DIIS extrapolation

      if(max_diis > 1) then
     
        n_diis = min(n_diis+1,max_diis)
        call DIIS_extrapolation(rcond,nO*nV*nOrb,nO*nV*nOrb,n_diis,r_diis,t_diis,res,amp)
     
      end if

      ! Compute quasiparticle energy
 
      eGW(p) = eHF(p)

      do i=nC+1,nO
        do a=1,nV-nR
          do q=nC+1,nOrb-nR
 
            eGW(p) = eGW(p) + vec(i,a,q)*amp(i,a,q)
  
          end do
        end do
      end do

      ! Dump results
 
      write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X)') &
        '|',nSCF,'|',eHF(p)*HaToeV,'|',eGW(p)*HaToeV,'|',Conv,'|'
 
    end do

    write(*,*)'-------------------------------------------------------------'
    !------------------------------------------------------------------------
    ! End of SCF loop
    !------------------------------------------------------------------------

    ! Did it actually converge?

    if(nSCF == maxSCF) then
 
      write(*,*)
      write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
      write(*,*)'                 Convergence failed                 '
      write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
      write(*,*)
 
      stop
 
    end if

    write(*,*)'-------------------------------------------------------------------------------'
    write(*,*)' CC-G0W0 calculation '
    write(*,*)'-------------------------------------------------------------------------------'
    write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
              '|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|'
    write(*,*)'-------------------------------------------------------------------------------'
 
    write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X)') &
    '|',p,'|',eHF(p)*HaToeV,'|',(eGW(p)-eHF(p))*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|'
    write(*,*)'-------------------------------------------------------------------------------'

  end do

end subroutine