saving worj in ccGW

src/GW/RGW.f90

@@ -162,12 +162,13 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_dii
 ! Perform CC-based G0W0 calculation
 !------------------------------------------------------------------------
 
-  doccG0W0 = .false.
+  doccG0W0 = .true.
 
   if(doccG0W0) then
     
     call wall_time(start_GW)
-    call ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
+    call ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,nS,ERI_MO,ENuc,ERHF,eHF)
+!   call ccRG0W0_TDA(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
     call wall_time(end_GW)
   
     t_GW = end_GW - start_GW
@@ -186,8 +187,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_dii
   if(doccGW) then
     
     call wall_time(start_GW)
-!   call ccRGW(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
+    call ccRGW(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
-!   call ccRGW_mat(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
     call wall_time(end_GW)
   
     t_GW = end_GW - start_GW

src/GW/ccRG0W0.f90

@@ -1,4 +1,4 @@
-subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
+subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,nS,ERI,ENuc,ERHF,eHF)
 
 ! CC-based GW module
 
@@ -17,6 +17,7 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
   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)   :: eHF(nOrb)
@@ -27,17 +28,29 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
   integer                       :: p,q,r,s
   integer                       :: i,j,k,l
   integer                       :: a,b,c,d
+  integer                       :: m
 
+  integer                       :: isp_W
+  logical                       :: TDA_W
+  logical                       :: dRPA
   integer                       :: nSCF
   double precision              :: Conv
+  double precision              :: EcRPA
 
   double precision,allocatable  :: Sig(:)
   double precision,allocatable  :: Z(:)
 
-  double precision,allocatable  :: del(:,:,:)
+  double precision,allocatable  :: del(:,:)
-  double precision,allocatable  :: vec(:,:,:)
+  double precision,allocatable  :: vec(:,:)
-  double precision,allocatable  :: res(:,:,:)
+  double precision,allocatable  :: res(:,:)
-  double precision,allocatable  :: amp(:,:,:)
+  double precision,allocatable  :: amp(:,:)
+
+  double precision,allocatable  :: Aph(:,:)
+  double precision,allocatable  :: Bph(:,:)
+  double precision,allocatable  :: Om(:)
+  double precision,allocatable  :: XpY(:,:)
+  double precision,allocatable  :: XmY(:,:)
+  double precision,allocatable  :: rho(:,:,:)
 
   integer                       :: n_diis
   double precision              :: rcond
@@ -55,16 +68,39 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
 
 ! Memory allocation
 
-  allocate(del(nO,nV,nOrb))
+  allocate(del(nS,nOrb))
-  allocate(vec(nO,nV,nOrb))
+  allocate(vec(nS,nOrb))
-  allocate(res(nO,nV,nOrb))
+  allocate(res(nS,nOrb))
-  allocate(amp(nO,nV,nOrb))
+  allocate(amp(nS,nOrb))
 
   allocate(Sig(nOrb))
   allocate(Z(nOrb))
 
-  allocate(r_diis(nO*nV*nOrb,max_diis))
+  allocate(r_diis(nS*nOrb,max_diis))
-  allocate(t_diis(nO*nV*nOrb,max_diis))
+  allocate(t_diis(nS*nOrb,max_diis))
+
+!-------------------!
+! Compute screening !
+!-------------------!
+
+  ! Spin manifold 
+
+  isp_W = 1
+  TDA_W = .false.
+  dRPA  = .true.
+
+  ! Memory allocation
+
+  allocate(Om(nS),Aph(nS,nS),Bph(nS,nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS))
+
+  call phLR_A(isp_W,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
+  call phLR_B(isp_W,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+
+  call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+
+  call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
+
+  deallocate(Aph,Bph,XpY,XmY)
 
 ! Initialization
 
@@ -87,50 +123,28 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
     r_diis(:,:) = 0d0
     rcond       = 0d0
 
-    amp(:,:,:) = 0d0
+    amp(:,:) = 0d0
-    res(:,:,:) = 0d0
+    res(:,:) = 0d0
  
-  ! Compute energy differences
+  ! Compute energy differences and coupling blocks
  
-    do i=nC+1,nO
+    do m=1,nS
-      do j=nC+1,nO
-        do a=1,nV-nR
-    
-          del(i,a,j) = eHF(i) + eHF(j) - eHF(nO+a) - eHF(p)
-    
-        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) - eHF(p)
-    
-        end do
-      end do
-    end do
-
-    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)
+        del(m,j) = Om(m) + eHF(j) - eHF(p)
+        vec(m,j) = sqrt(2d0)*rho(p,j,m)
     
       end do
     end do
-    end do
     
-    do i=nC+1,nO
+    do m=1,nS
-      do a=1,nV-nR
       do b=1,nV-nR
     
-          vec(i,a,nO+b) = sqrt(2d0)*ERI(p,i,nO+b,nO+a)
+        del(m,nO+b) = Om(m) + eHF(nO+b) - eHF(p)
+        vec(m,nO+b) = sqrt(2d0)*rho(p,nO+b,m)
     
       end do
     end do
-    end do
 
    !----------------------!
    ! Loop over amplitudes !
@@ -152,41 +166,23 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
  
       ! Compute residual for 2h1p sector
  
-      res(:,:,:) = vec(:,:,:) + (del(:,:,:) - Sig(p))*amp(:,:,:)
+      res(:,:) = vec(:,:) + (del(:,:) - Sig(p))*amp(:,:)
  
-      do i=nC+1,nO
+      do m=1,nS
-        do a=1,nV-nR
         do j=nC+1,nO
  
-            do k=nC+1,nO
+          res(m,j) = res(m,j) - Om(m)*amp(m,j) 
-              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) &
-                                        - 2d0*ERI(i,nO+c,nO+a,k)*amp(k,c,j)
- 
-              end do
-            end do
- 
-          end do
         end do
       end do
  
       ! Compute residual for 2p1h sector
  
-      do i=nC+1,nO
+      do m=nC+1,nO
-        do a=1,nV-nR
         do b=1,nV-nR
  
-            do k=nC+1,nO
+          res(m,nO+b) = res(m,nO+b) + Om(m)*amp(m,nO+b) 
-              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) &
-                                              + 2d0*ERI(nO+b,k,i,nO+c)*amp(k,a,nO+c) 
- 
-              end do
-            end do
- 
-          end do
         end do
       end do
   
@@ -196,14 +192,14 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
     
       ! Update amplitudes
 
-      amp(:,:,:) = amp(:,:,:) - res(:,:,:)/del(:,:,:)
+      amp(:,:) = amp(:,:) - res(:,:)/del(:,:)
  
       ! 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)
+        call DIIS_extrapolation(rcond,nS*nOrb,nS*nOrb,n_diis,r_diis,t_diis,res,amp)
      
       end if
 
@@ -211,15 +207,13 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
  
       Sig(p) = 0d0
 
-      do i=nC+1,nO
+      do m=1,nS
-        do a=1,nV-nR
         do q=nC+1,nOrb-nR
  
-            Sig(p) = Sig(p) + vec(i,a,q)*amp(i,a,q)
+          Sig(p) = Sig(p) + vec(m,q)*amp(m,q)
   
         end do
       end do
-      end do
 
       ! Dump results
  

src/GW/ccRG0W0_TDA.f90

@@ -0,0 +1,263 @@
+subroutine ccRG0W0_TDA(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
+
+! CC-based GW module (TDA screening)
+
+  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,allocatable  :: Sig(:)
+  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(*,*)
+
+! Memory allocation
+
+  allocate(del(nO,nV,nOrb))
+  allocate(vec(nO,nV,nOrb))
+  allocate(res(nO,nV,nOrb))
+  allocate(amp(nO,nV,nOrb))
+
+  allocate(Sig(nOrb))
+  allocate(Z(nOrb))
+
+  allocate(r_diis(nO*nV*nOrb,max_diis))
+  allocate(t_diis(nO*nV*nOrb,max_diis))
+
+! Initialization
+
+  Sig(:) = 0d0
+  Z(:)   = 1d0
+
+!-------------------------!
+! Main loop over orbitals !
+!-------------------------!
+
+  do p=nO,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 j=nC+1,nO
+        do a=1,nV-nR
+    
+          del(i,a,j) = eHF(i) + eHF(j) - eHF(nO+a) - eHF(p)
+    
+        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) - eHF(p)
+    
+        end do
+      end do
+    end do
+
+    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)') &
+              '|','#','|','Sig_c (eV)','|','e_GW (eV)','|','Conv','|'
+    write(*,*)'-------------------------------------------------------------'
+ 
+    do while(Conv > thresh .and. nSCF < maxSCF)
+ 
+     ! Increment 
+ 
+      nSCF = nSCF + 1
+ 
+      ! Compute residual for 2h1p sector
+ 
+      res(:,:,:) = vec(:,:,:) + (del(:,:,:) - Sig(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)  
+!                                       - 2d0*ERI(i,nO+c,nO+a,k)*amp(k,c,j)
+ 
+              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) 
+!                                             + 2d0*ERI(nO+b,k,i,nO+c)*amp(k,a,nO+c) 
+ 
+              end do
+            end do
+ 
+          end do
+        end do
+      end do
+  
+      ! Check convergence 
+ 
+      Conv = maxval(abs(res))
+    
+      ! Update amplitudes
+
+      amp(:,:,:) = amp(:,:,:) - res(:,:,:)/del(:,:,:)
+ 
+      ! 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
+ 
+      Sig(p) = 0d0
+
+      do i=nC+1,nO
+        do a=1,nV-nR
+          do q=nC+1,nOrb-nR
+ 
+            Sig(p) = Sig(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,'|',Sig(p)*HaToeV,'|',(eHF(p)+Sig(p))*HaToeV,'|',Conv,'|'
+ 
+    end do
+
+    write(*,*)'-------------------------------------------------------------'
+    write(*,*)
+    !------------------------------------------------------------------------
+    ! End of SCF loop
+    !------------------------------------------------------------------------
+
+    ! Did it actually converge?
+
+    if(nSCF == maxSCF) then
+ 
+      write(*,*)
+      write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
+      write(*,*)'                 Convergence failed                 '
+      write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
+      write(*,*)
+ 
+    end if
+
+    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,A15,1X,A1,1X)') &
+              '|','Orb','|','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,'|',Sig(p)*HaToeV,'|',Z(p),'|',(eHF(p)+Sig(p))*HaToeV,'|'
+    write(*,*)'-------------------------------------------------------------------------------'
+    write(*,*)
+
+  end do
+
+end subroutine 

src/GW/ccRG0W0_mat.f90

@@ -1,278 +0,0 @@
-subroutine ccRG0W0_mat(maxSCF,thresh,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)            :: 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
-  integer                       :: i,j,k,l
-  integer                       :: a,b,c,d
-
-  integer                       :: nSCF
-  double precision              :: Conv
-
-  double precision,allocatable  :: delta_2h1p(:)
-  double precision,allocatable  :: delta_2p1h(:)
-
-  double precision,allocatable  :: C_2h1p(:,:)
-  double precision,allocatable  :: C_2p1h(:,:)
-
-  double precision,allocatable  :: V_2h1p(:)
-  double precision,allocatable  :: V_2p1h(:)
-
-  double precision,allocatable  :: r_2h1p(:)
-  double precision,allocatable  :: r_2p1h(:)
-
-  double precision,allocatable  :: t_2h1p(:)
-  double precision,allocatable  :: t_2p1h(:)
-
-  double precision,allocatable  :: eGW(:)
-  double precision,allocatable  :: Z(:)
-
-  double precision,external     :: Kronecker_delta
-  
-  integer                       :: n1h,n1p,n1h1p,n2h1p,n2p1h
-  integer                       :: ija,klc,iab,kcd
-
-
-! Hello world
-
-  write(*,*)
-  write(*,*)'*****************************'
-  write(*,*)'* CC-based G0W0 Calculation *'
-  write(*,*)'*****************************'
-  write(*,*)
- 
-! Form energy denominator and guess amplitudes
-
-  n1h = nO
-  n1p = nV
-  n1h1p = n1h*n1p
-  n2h1p = n1h*n1h1p
-  n2p1h = n1p*n1h1p
-
-  allocate(delta_2h1p(n2h1p),delta_2p1h(n2p1h))
-  allocate(C_2h1p(n2h1p,n2h1p),C_2p1h(n2p1h,n2p1h))
-  allocate(V_2h1p(n2h1p),V_2p1h(n2p1h))
-  allocate(t_2h1p(n2h1p),t_2p1h(n2p1h))
-  allocate(r_2h1p(n2h1p),r_2p1h(n2p1h))
-  allocate(eGW(nOrb),Z(nOrb))
-
-  ! Compute C2h1p and C2p1h
-
-  ija = 0
-  do i=nC+1,nO
-    do j=nC+1,nO
-      do a=nO+1,nOrb-nR
-        ija = ija + 1
-
-        klc = 0
-        do k=nC+1,nO
-          do l=nC+1,nO
-            do c=nO+1,nOrb-nR
-              klc = klc + 1
-
-              C_2h1p(ija,klc) = ((eHF(i) + eHF(j) - eHF(a))*Kronecker_delta(j,l)*Kronecker_delta(a,c) &
-                              - 2d0*ERI(j,c,a,l))*Kronecker_delta(i,k)
-
-            end do
-          end do
-        end do
-
-      end do
-    end do
-  end do
-
-  iab = 0
-  do i=nC+1,nO
-    do a=nO+1,nOrb-nR
-      do b=nO+1,nOrb-nR
-        iab = iab + 1
-
-        kcd = 0
-        do k=nC+1,nO
-          do c=nO+1,nOrb-nR
-            do d=nO+1,nOrb-nR
-              kcd = kcd + 1
-
-              C_2p1h(iab,kcd) = ((eHF(a) + eHF(b) - eHF(i))*Kronecker_delta(i,k)*Kronecker_delta(a,c) &
-                              + 2d0*ERI(a,k,i,c))*Kronecker_delta(b,d)
-
-            end do
-          end do
-        end do
-
-      end do
-    end do
-  end do
-
-!-------------------------!
-! Main loop over orbitals !
-!-------------------------!
-
-  eGW(:) = eHF(:)
-
-  do p=nO,nO
-
-    ! Initialization
- 
-    Conv = 1d0
-    nSCF =  0
- 
-    t_2h1p(:) = 0d0
-    t_2p1h(:) = 0d0
- 
-    ! Compute energy differences
- 
-    ija = 0
-    do i=nC+1,nO
-      do j=nC+1,nO
-        do a=nO+1,nOrb-nR
-          ija = ija + 1
- 
-          delta_2h1p(ija) = eHF(i) + eHF(j) - eHF(a) - eHF(p)
- 
-        end do
-      end do
-    end do
- 
-    iab = 0
-    do i=nC+1,nO
-      do a=nO+1,nOrb-nR
-        do b=nO+1,nOrb-nR
-          iab = iab + 1
- 
-          delta_2p1h(iab) = eHF(a) + eHF(b) - eHF(i) - eHF(p)
- 
-        end do
-      end do
-    end do
-
-    klc = 0
-    do k=nC+1,nO
-      do l=nC+1,nO
-        do c=nO+1,nOrb-nR
-          klc = klc + 1
- 
-          V_2h1p(klc) = sqrt(2d0)*ERI(p,c,k,l)
- 
-        end do
-      end do
-    end do
- 
-    kcd = 0
-    do k=nC+1,nO
-      do c=nO+1,nOrb-nR
-        do d=nO+1,nOrb-nR
-          kcd = kcd + 1
- 
-          V_2p1h(kcd) = sqrt(2d0)*ERI(p,k,d,c)
- 
-        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,A10,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') &
-              '|','#','|','HF','|','G0W0','|','Conv','|'
-    write(*,*)'----------------------------------------------'
- 
-    do while(Conv > thresh .and. nSCF < maxSCF)
- 
-     ! Increment 
- 
-      nSCF = nSCF + 1
-
-      ! Compute residual for 2h1p sector
- 
-      r_2h1p = V_2h1p + matmul(C_2h1p,t_2h1p) - t_2h1p*eHF(p) &
-             - dot_product(t_2h1p,V_2h1p)*t_2h1p - t_2h1p*dot_product(V_2p1h,t_2p1h)
-
-      ! Compute residual for 2p1h sector
- 
-      r_2p1h = V_2p1h + matmul(C_2p1h,t_2p1h) - t_2p1h*eHF(p) &
-             - t_2p1h*dot_product(V_2h1p,t_2h1p) - dot_product(t_2p1h,V_2p1h)*t_2p1h
-  
-      !  Check convergence 
- 
-      Conv = max(maxval(abs(r_2h1p)),maxval(abs(r_2p1h)))
-    
-      ! Update amplitudes
- 
-      t_2h1p(:) = t_2h1p(:) - r_2h1p(:)/delta_2h1p(:)
-      t_2p1h(:) = t_2p1h(:) - r_2p1h(:)/delta_2p1h(:)
- 
-      ! Compute quasiparticle energies
- 
-      eGW(p) = eHF(p) + dot_product(V_2h1p,t_2h1p) + dot_product(V_2p1h,t_2p1h)
- 
-      ! Renormalization factor
- 
-      Z(:) = 1d0
- 
-      ! Dump results
- 
-      write(*,'(1X,A1,1X,I3,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X,F10.6,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.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
-    '|',p,'|',eHF(p)*HaToeV,'|',(eGW(p)-eHF(p))*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|'
-    write(*,*)'-------------------------------------------------------------------------------'
-
-  end do
-
-end subroutine