complex linearized evGW@cRHF

src/GW/complex_RGW.f90

@@ -1,4 +1,4 @@
-subroutine complex_RGW(dotest,docG0W0,maxSCF,thresh,max_diis,doACFDT,                              &
+subroutine complex_RGW(dotest,docG0W0,doevGW,maxSCF,thresh,max_diis,doACFDT,                    &
                exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
                linearize,eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,              &
                S,X,T,V,Hc,ERI_AO,ERI_MO,CAP_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
@@ -12,7 +12,7 @@ subroutine complex_RGW(dotest,docG0W0,maxSCF,thresh,max_diis,doACFDT,
 
   logical,intent(in)            :: dotest
 
-  logical,intent(in)            :: docG0W0
+  logical,intent(in)            :: docG0W0,doevGW
 
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
@@ -82,4 +82,22 @@ subroutine complex_RGW(dotest,docG0W0,maxSCF,thresh,max_diis,doACFDT,
     write(*,*)
 
   end if
+
+!------------------------------------------------------------------------
+! Perform evGW calculation
+!------------------------------------------------------------------------
+
+  if(doevGW) then
+
+    call wall_time(start_GW)
+    call complex_evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
+               singlet,triplet,linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
+    call wall_time(end_GW)
+
+    t_GW = end_GW - start_GW
+    write(*,'(A65,1X,F9.3,A8)') 'Total wall time for evGW = ',t_GW,' seconds'
+    write(*,*)
+
+  end if
+
 end subroutine

src/GW/complex_cRG0W0.f90

@@ -59,8 +59,6 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,
   double precision,allocatable  :: Re_Z(:)
   double precision,allocatable  :: Im_Z(:)
   complex*16,allocatable        :: Om(:)
-  double precision,allocatable  :: Re_Om(:)
-  double precision,allocatable  :: Im_Om(:)
   complex*16,allocatable        :: XpY(:,:)
   complex*16,allocatable        :: XmY(:,:)
   complex*16,allocatable        :: rho(:,:,:)
@@ -106,7 +104,7 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,
 ! Memory allocation
 
   allocate(Aph(nS,nS),Bph(nS,nS),Re_SigC(nOrb),Im_SigC(nOrb),Re_Z(nOrb),Im_Z(nOrb),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS), & 
-           Re_eGW(nOrb),Im_eGW(nOrb),Re_eGWlin(nOrb),Im_eGWlin(nOrb),Re_eHF(nOrb),Im_eHF(nOrb),Re_Om(nS),Im_Om(nS))
+           Re_eGW(nOrb),Im_eGW(nOrb),Re_eGWlin(nOrb),Im_eGWlin(nOrb),Re_eHF(nOrb),Im_eHF(nOrb),)
   Re_eHF(:) = real(eHF(:))
   Im_eHF(:) = aimag(eHF(:))
 !-------------------!
@@ -117,8 +115,6 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,
   if(.not.TDA_W) call complex_phRLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
   call complex_phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
-  Re_Om(:) = real(Om(:))
-  Im_Om(:) = aimag(Om(:))
 
   !if(print_W) call print_excitation_energies('phRPA@RHF','singlet',nS,Om)
 
@@ -152,7 +148,6 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,
 
     write(*,*) ' *** Quasiparticle energies obtained by root search *** '
     write(*,*)
-    write(*,*) "ROOT SEARCH NOT IMPLEMENTED YET" 
     call complex_RGW_QP_graph(doSRG,eta,flow,nOrb,nC,nO,nV,nR,nS,       &
             Re_eHF,Im_eHF,Om,rho,Re_eGWlin,Im_eGWlin,Re_eHF,Im_eHF,     &
             Re_eGW,Im_eGW,Re_Z,Im_Z)
@@ -160,7 +155,7 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,
 
 ! Plot self-energy, renormalization factor, and spectral function
 !
-!  if(plot_self) call RGW_plot_self_energy(nOrb,eta,nC,nO,nV,nR,nS,eHF,eHF,Om,rho)
+  if(plot_self) call complex_RGW_plot_self_energy(nOrb,eta,nC,nO,nV,nR,nS,eHF,eHF,Om,rho)
 !
 !! Cumulant expansion 
 !

src/GW/complex_evRGW.f90

@@ -0,0 +1,292 @@
+subroutine complex_evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & 
+                 singlet,triplet,linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
+
+! Perform self-consistent eigenvalue-only GW calculation
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  logical,intent(in)            :: dotest
+
+  integer,intent(in)            :: maxSCF
+  integer,intent(in)            :: max_diis
+  double precision,intent(in)   :: thresh
+  double precision,intent(in)   :: ENuc
+  double precision,intent(in)   :: ERHF
+  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)            :: singlet
+  logical,intent(in)            :: triplet
+  logical,intent(in)            :: linearize
+  double precision,intent(in)   :: eta
+  logical,intent(in)            :: doSRG
+
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nOrb
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  complex*16,intent(in)         :: eHF(nOrb)
+  complex*16,intent(in)         :: ERI(nOrb,nOrb,nOrb,nOrb)
+  complex*16,intent(in)         :: dipole_int(nOrb,nOrb,ncart)
+
+! Local variables
+
+  logical                       :: dRPA = .true.
+  integer                       :: ispin
+  integer                       :: nSCF
+  integer                       :: n_diis
+  double precision              :: flow
+  double precision              :: rcond
+  double precision              :: Conv
+  complex*16                    :: EcRPA
+  complex*16                    :: EcBSE(nspin)
+  complex*16                    :: EcGM
+  double precision              :: alpha
+  complex*16,allocatable        :: Aph(:,:)
+  complex*16,allocatable        :: Bph(:,:)
+  complex*16,allocatable        :: error_diis(:,:)
+  complex*16,allocatable        :: e_diis(:,:)
+  complex*16,allocatable        :: eGW(:)
+  complex*16,allocatable        :: eOld(:)
+  double precision,allocatable  :: Re_eGW(:)
+  double precision,allocatable  :: Im_eGW(:)
+  double precision,allocatable  :: Re_eOld(:)
+  double precision,allocatable  :: Im_eOld(:)
+  double precision,allocatable  :: Re_eHF(:)
+  double precision,allocatable  :: Im_eHF(:)
+  double precision,allocatable  :: Re_Z(:)
+  double precision,allocatable  :: Im_Z(:)
+  double precision,allocatable  :: Re_SigC(:)
+  double precision,allocatable  :: Im_SigC(:)
+  complex*16,allocatable        :: Om(:)
+  complex*16,allocatable        :: XpY(:,:)
+  complex*16,allocatable        :: XmY(:,:)
+  complex*16,allocatable        :: rho(:,:,:)
+
+! Hello world
+
+  write(*,*)
+  write(*,*)'*******************************'
+  write(*,*)'* Restricted evGW Calculation *'
+  write(*,*)'*******************************'
+  write(*,*)
+
+! TDA for W
+
+  if(TDA_W) then 
+    write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
+    write(*,*)
+  end if
+
+! SRG regularization
+
+  flow = 100d0
+
+  if(doSRG) then
+
+    write(*,*) '*** SRG regularized evGW scheme ***'
+    write(*,*)
+
+  end if
+
+! Memory allocation
+
+  allocate(Aph(nS,nS),Bph(nS,nS),Re_eGW(nOrb),Im_eGW(nOrb),eGW(nOrb),Re_eOld(nOrb),Im_eOld(nOrb),&
+           eOld(nOrb),Re_Z(nOrb),Im_Z(nOrb),Re_SigC(nOrb),Im_SigC(nOrb), & 
+           Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS),error_diis(nOrb,max_diis),e_diis(nOrb,max_diis),&
+           Re_eHF(nOrb),Im_eHF(nOrb))
+
+! Initialization
+
+  nSCF            = 0
+  ispin           = 1
+  n_diis          = 0
+  Conv            = 1d0
+  e_diis(:,:)     = 0d0
+  error_diis(:,:) = 0d0
+  Re_eHF(:)       = real(eHF(:))
+  Im_eHF(:)       = aimag(eHF(:))
+  eGW(:)          = eHF(:)
+  Re_eGW(:)       = Re_eHF(:)
+  Im_eGW(:)       = Im_eHF(:)
+  eOld(:)         = eGW(:)
+  Re_eOld(:)      = Re_eGW(:)
+  Im_eOld(:)      = Im_eGW(:)
+  Re_Z(:)         = 1d0
+  Im_Z(:)         = 0d0
+  rcond           = 0d0
+  
+!------------------------------------------------------------------------
+! Main loop
+!------------------------------------------------------------------------
+
+  do while(Conv > thresh .and. nSCF <= maxSCF)
+
+   ! Compute screening
+
+                   call complex_phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
+    if(.not.TDA_W) call complex_phRLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+
+    call complex_phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+
+    ! Compute spectral weights
+
+    call complex_RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
+
+    ! Compute correlation part of the self-energy 
+
+    call complex_RGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_eGW,Im_eGW,Om,rho,&
+            EcGM,Re_SigC,Im_SigC,Re_Z,Im_Z)
+
+    ! Solve the quasi-particle equation
+
+    if(linearize) then 
+ 
+      write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
+      write(*,*)
+
+      Re_eGW(:) = Re_eHF(:) + Re_SigC(:)
+      Im_eGW(:) = Im_eHF(:) + Im_SigC(:)
+      eGW = cmplx(Re_eGW + Im_eGW,kind=8)
+    else 
+
+      write(*,*) ' *** Quasiparticle energies obtained by root search *** '
+      write(*,*)
+  
+      call complex_RGW_QP_graph(doSRG,eta,flow,nOrb,nC,nO,nV,nR,nS,Re_eHF,Im_eHF,Om,&
+              rho,Re_eOld,Im_eOld,Re_eGW,Im_eGW,Re_Z,Im_Z)
+ 
+      eGW = cmplx(Re_eGW + Im_eGW,kind=8)
+    end if
+
+    ! Convergence criteria
+
+    Conv = maxval(abs(eGW - eOld))
+
+    ! Print results
+
+    call print_complex_evRGW(nBas,nO,nSCF,Conv,Re_eHF,Im_eHF,ENuc,ERHF,Re_SigC,Im_SigC,&
+            Re_Z,Im_Z,Re_eGW,Im_eGW,EcRPA,EcGM)
+    ! DIIS extrapolation
+
+    if(max_diis > 1) then
+
+      n_diis = min(n_diis+1,max_diis)
+      call complex_DIIS_extrapolation(rcond,nOrb,nOrb,n_diis,error_diis,e_diis,eGW-eOld,eGW)
+
+    end if
+
+    ! Save quasiparticles energy for next cycle
+    Re_eGW(:) = real(eGW(:))
+    Im_eGW(:) = aimag(eGW(:))
+    eOld(:) = eGW(:)
+    Re_eOld(:) = real(eOld(:))
+    Im_eOld(:) = aimag(eOld(:))
+
+    ! Increment
+
+    nSCF = nSCF + 1
+
+  end do
+!------------------------------------------------------------------------
+! End main loop
+!------------------------------------------------------------------------
+
+! Did it actually converge?
+
+  if(nSCF == maxSCF+1) then
+
+    write(*,*)
+    write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
+    write(*,*)'                 Convergence failed                 '
+    write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
+    write(*,*)
+
+    stop
+
+  end if
+
+!--------------------!
+! Cumulant expansion !
+!--------------------!
+!
+!! call RGWC(dotest,eta,nOrb,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,eGW,Z)
+!
+!! Deallocate memory
+!
+!  deallocate(Aph,Bph,eOld,Z,SigC,Om,XpY,XmY,rho,error_diis,e_diis)
+!
+!! Perform BSE calculation
+!
+!  if(dophBSE) then
+!
+!    call RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta, & 
+!                   nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eGW,eGW,EcBSE)
+!
+!    write(*,*)
+!    write(*,*)'-------------------------------------------------------------------------------'
+!    write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@evGW@RHF correlation energy (singlet) = ',EcBSE(1),' au'
+!    write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@evGW@RHF correlation energy (triplet) = ',EcBSE(2),' au'
+!    write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@evGW@RHF correlation energy           = ',sum(EcBSE),' au'
+!    write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@evGW@RHF total       energy           = ',ENuc + ERHF + sum(EcBSE),' au'
+!    write(*,*)'-------------------------------------------------------------------------------'
+!    write(*,*)
+!
+!!   Compute the BSE correlation energy via the adiabatic connection 
+!
+!    if(doACFDT) then
+!
+!      call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcBSE)
+!
+!      write(*,*)
+!      write(*,*)'-------------------------------------------------------------------------------'
+!      write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@evGW@RHF correlation energy (singlet) = ',EcBSE(1),' au'
+!      write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@evGW@RHF correlation energy (triplet) = ',EcBSE(2),' au'
+!      write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@evGW@RHF correlation energy           = ',sum(EcBSE),' au'
+!      write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@evGW@RHF total       energy           = ',ENuc + ERHF + sum(EcBSE),' au'
+!      write(*,*)'-------------------------------------------------------------------------------'
+!      write(*,*)
+!
+!    end if
+!
+!  end if
+!
+!  if(doppBSE) then
+!
+!    call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
+!
+!    write(*,*)
+!    write(*,*)'-------------------------------------------------------------------------------'
+!    write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@evGW@RHF correlation energy (singlet) = ',EcBSE(1),' au'
+!    write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@evGW@RHF correlation energy (triplet) = ',EcBSE(2),' au'
+!    write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@evGW@RHF correlation energy           = ',sum(EcBSE),' au'
+!    write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@evGW@RHF total       energy           = ',ENuc + ERHF + sum(EcBSE),' au'
+!    write(*,*)'-------------------------------------------------------------------------------'
+!    write(*,*)
+!
+!  end if
+!
+!! Testing zone
+!
+!  if(dotest) then
+!
+!    call dump_test_value('R','evGW correlation energy',EcRPA)
+!    call dump_test_value('R','evGW HOMO energy',eGW(nO))
+!    call dump_test_value('R','evGW LUMO energy',eGW(nO+1))
+!
+!  end if
+!
+end subroutine 

src/GW/print_complex_evRGW.f90

@@ -0,0 +1,67 @@
+subroutine print_complex_evRGW(nBas,nO,nSCF,Conv,Re_eHF,Im_eHF,ENuc,ERHF,Re_SigC,Im_SigC,Re_Z,Im_Z,Re_eGW,Im_eGW,EcRPA,EcGM)
+
+! Print one-electron energies and other stuff for G0W0
+
+  implicit none
+  include 'parameters.h'
+
+  integer,intent(in)                 :: nBas,nO,nSCF
+  double precision,intent(in)        :: Conv
+  double precision,intent(in)        :: ENuc
+  complex*16,intent(in)              :: ERHF
+  complex*16,intent(in)              :: EcRPA
+  complex*16,intent(in)              :: EcGM
+  double precision,intent(in)        :: Re_eHF(nBas)
+  double precision,intent(in)        :: Im_eHF(nBas)
+  double precision,intent(in)        :: Re_SigC(nBas)
+  double precision,intent(in)        :: Im_SigC(nBas)
+  double precision,intent(in)        :: Re_Z(nBas)
+  double precision,intent(in)        :: Im_Z(nBas)
+  double precision,intent(in)        :: Re_eGW(nBas)
+  double precision,intent(in)        :: Im_eGW(nBas)
+
+  integer                            :: p,index_homo,index_lumo
+  double precision                   :: Re_eHOMO,Re_eLUMO,Im_eHOMO,Im_eLUMO,Re_Gap,Im_Gap
+ 
+! HOMO and LUMO
+
+  index_homo = maxloc(Re_eGW(1:nO),1)
+  Re_eHOMO = Re_eGW(index_homo)
+  Im_eHOMO = Im_eGW(index_homo)
+  index_lumo = minloc(Re_eGW(nO+1:nBas),1) + nO
+  Re_eLUMO = Re_eGW(index_lumo)
+  Im_eLUMO = Im_eGW(index_lumo)
+  Re_Gap = Re_eLUMO-Re_eHOMO
+  Im_Gap = Im_eLUMO-Im_eHOMO
+
+! Dump results
+  write(*,*)'-------------------------------------------------------------------------------'
+  if(nSCF < 10) then
+    write(*,'(1X,A20,I1,A1,I1,A17)')' Self-consistent evG',nSCF,'W',nSCF,'@cRHF calculation'
+  elseif(nSCF < 100) then
+    write(*,'(1X,A20,I2,A1,I2,A17)')' Self-consistent evG',nSCF,'W',nSCF,'@cRHF calculation'
+  else
+    write(*,'(1X,A20,I3,A1,I3,A17)')' Self-consistent evG',nSCF,'W',nSCF,'@cRHF calculation'
+  end if
+  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)') &
+            '|','#','|','Re(e_HF) (eV)','|','Re(Sig_GW) (eV)','|','Re(Z)','|','Re(e_GW) (eV)','|'
+  write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
+            '|','#','|','Im(e_HF) (eV)','|','Im(Sig_GW) (eV)','|','Im(Z)','|','Im(e_GW) (eV)','|'
+  write(*,*)'-------------------------------------------------------------------------------'
+  do p=1,nBas
+    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,'|',Re_eHF(p)*HaToeV,'|',Re_SigC(p)*HaToeV,'|',Re_Z(p),'|',Re_eGW(p)*HaToeV,'|'
+    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,'|',Im_eHF(p)*HaToeV,'|',Im_SigC(p)*HaToeV,'|',Im_Z(p),'|',Im_eGW(p)*HaToeV,'|'
+    if(p==nO) then
+      write(*,*)'-------------------------------------------------------------------------------'
+    end if
+  write(*,*)'-------------------------------------------------------------------------------'
+  end do
+  write(*,*)
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,'(2X,A10,I3)')   'Iteration ',nSCF
+  write(*,'(2X,A14,F15.5)')'Convergence = ',Conv
+  write(*,*)
+end subroutine 

src/QuAcK/QuAcK.f90

@@ -131,7 +131,8 @@ program QuAcK
                     docG0W0,docG0F2,                       &
                     doRtest,doUtest,doGtest)
   doCAP = docG0W0 .or. docG0F2 .or. docRHF ! Add different cases if they need CAP
-
+  docG0W0 = docG0W0 .or. (doG0W0 .and. docRHF)
+  docG0F2 = docG0F2 .or. (doG0F2 .and. docRHF)
 !--------------------------!
 ! Read options for methods !
 !--------------------------!

src/QuAcK/RQuAcK.f90

@@ -418,7 +418,7 @@ doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3
 
   if(doGW .and. docRHF) then
     call wall_time(start_GW)
-    call complex_RGW(dotest,docG0W0,maxSCF_GW,thresh_GW,max_diis_GW,        & 
+    call complex_RGW(dotest,docG0W0,doevGW,maxSCF_GW,thresh_GW,max_diis_GW,        & 
              doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
              lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T,               &
              V,Hc,ERI_AO,complex_ERI_MO,complex_CAP_MO,dipole_int_AO,complex_dipole_int_MO,complex_PHF,complex_cHF,complex_eHF)

src/utils/complex_DIIS_extrapolation.f90

@@ -49,7 +49,7 @@ subroutine complex_DIIS_extrapolation(rcond,n_err,n_e,n_diis,error,e,error_in,e_
   call complex_linear_solve(n_diis+1,A,b,w,rcond)
 ! Check condition number
   if (rcond < 1.0d-10) then
-    write(*,*) "!!! DIIS system ill-conditioned: rcond = ", rcond, " - Skipping DIIS update !!!"
+    !write(*,*) "!!! DIIS system ill-conditioned: rcond = ", rcond, " - Skipping DIIS update !!!"
   else
     ! Perform extrapolation only if the system is well-conditioned
     e_inout(:) = matmul(w(1:n_diis),transpose(e(:,1:n_diis)))