debug qsGW and qsGF

input/methods

@@ -1,5 +1,5 @@
 # RHF UHF KS  MOM
-  T   F   F   F
+  F   T   F   F
 # MP2* MP3 MP2-F12
   F   F   F
 # CCD DCD CCSD CCSD(T) 
@@ -11,9 +11,9 @@
 # RPA* RPAx* ppRPA 
   F   F    F    
 # G0F2* evGF2* qsGF2* G0F3 evGF3
-  F     F     F      F    F
+  F     F     T      F    F
 # G0W0* evGW* qsGW*
-  T    F    F
+  F    F    F
 # G0T0 evGT qsGT 
   F    F    F
 # MCMP2

input/options

@@ -1,5 +1,5 @@
 # HF: maxSCF thresh   DIIS n_diis guess_type ortho_type mix_guess stability
-       128    0.0000001  T    5      1          1         T         F
+       1024    0.0000001  T    5      2          1         T         F
 # MP: 
       
 # CC: maxSCF thresh   DIIS n_diis
@@ -7,12 +7,12 @@
 # spin: TDA singlet triplet spin_conserved spin_flip 
          F   T       T       T              T 
 # GF: maxSCF thresh  DIIS n_diis lin eta renorm 
-       256  0.00001 T    5      T   0.001  3      
+       256  0.00001 T    5      T   0.0   3      
 # GW/GT: maxSCF thresh  DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 
-         256   0.00001   T    5     T  0.0  F      F     F     F   F  
+         256   0.0000001   T    5     T  0.0  F      F     F     F   F  
 # ACFDT: AC Kx  XBS
          F  F   T
 # BSE: BSE dBSE dTDA evDyn
-        T    T    T    F
+        F    T    T    F
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
          1000000 100000 10    0.3 10000 1234  T

mol/h2.xyz

@@ -1,4 +1,4 @@
 2
 
 H 0.0 0.0 0.0
-H 0.0 0.0 0.71
+H 0.0 0.0 1.0

src/GF/G0W0.f90

@@ -1,238 +0,0 @@
-subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & 
-                dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,            & 
-                nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0W0)
-
-! 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)            :: COHSEX
-  logical,intent(in)            :: SOSEX
-  logical,intent(in)            :: BSE
-  logical,intent(in)            :: TDA_W
-  logical,intent(in)            :: TDA
-  logical,intent(in)            :: dBSE
-  logical,intent(in)            :: dTDA
-  logical,intent(in)            :: evDyn
-  logical,intent(in)            :: singlet
-  logical,intent(in)            :: triplet
-  logical,intent(in)            :: linearize
-  double precision,intent(in)   :: eta
-
-  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
-  double precision,intent(in)   :: ENuc
-  double precision,intent(in)   :: ERHF
-  double precision,intent(in)   :: ERI_AO(nBas,nBas,nBas,nBas)
-  double precision,intent(in)   :: ERI_MO(nBas,nBas,nBas,nBas)
-  double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
-  double precision,intent(in)   :: Vxc(nBas)
-  double precision,intent(in)   :: eHF(nBas)
-  double precision,intent(in)   :: cHF(nBas,nBas)
-  double precision,intent(in)   :: PHF(nBas,nBas)
-
-! Local variables
-
-  logical                       :: print_W = .true.
-  integer                       :: ispin
-  double precision              :: EcRPA
-  double precision              :: EcBSE(nspin)
-  double precision              :: EcAC(nspin)
-  double precision              :: EcGM
-  double precision,allocatable  :: SigX(:)
-  double precision,allocatable  :: SigC(:)
-  double precision,allocatable  :: Z(:)
-  double precision,allocatable  :: OmRPA(:)
-  double precision,allocatable  :: XpY_RPA(:,:)
-  double precision,allocatable  :: XmY_RPA(:,:)
-  double precision,allocatable  :: rho_RPA(:,:,:)
-
-  double precision,allocatable  :: eG0W0lin(:)
-
-! Output variables
-
-  double precision              :: eG0W0(nBas)
-
-! Hello world
-
-  write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|          One-shot G0W0 calculation           |'
-  write(*,*)'************************************************'
-  write(*,*)
-
-! Initialization
-
-  EcRPA = 0d0
-
-! SOSEX correction
-
-  if(SOSEX) then 
-    write(*,*) 'SOSEX correction activated but BUG!'
-    stop
-  end if
-
-! COHSEX approximation
-
-  if(COHSEX) then 
-    write(*,*) 'COHSEX approximation activated!'
-    write(*,*)
-  end if
-
-! 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 = 1
-
-! Memory allocation
-
-  allocate(SigC(nBas),SigX(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS),eG0W0lin(nBas))
-
-!-------------------!
-! Compute screening !
-!-------------------!
-
-  call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0, & 
-                       eHF,ERI_MO,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
-
-  if(print_W) call print_excitation('RPA@HF      ',ispin,nS,OmRPA)
-
-!--------------------------!
-! Compute spectral weights !
-!--------------------------!
-
-  call excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA)
-
-!------------------------!
-! Compute GW self-energy !
-!------------------------!
-
-  call self_energy_exchange_diag(nBas,cHF,PHF,ERI_AO,SigX)
-  call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
-
-!--------------------------------!
-! Compute renormalization factor !
-!--------------------------------!
-
-  call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
-
-!-----------------------------------!
-! Solve the quasi-particle equation !
-!-----------------------------------!
-
-  eG0W0lin(:) = eHF(:) + Z(:)*(SigX(:) + SigC(:) - Vxc(:))
-
-  ! Linearized or graphical solution?
-
-  if(linearize) then 
- 
-    write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
-    write(*,*)
-
-    eG0W0(:) = eG0W0lin(:)
-
-  else 
-
-    write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
-    write(*,*)
-  
-    call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,OmRPA,rho_RPA,eG0W0lin,eG0W0)
-
-    ! Find all the roots of the QP equation if necessary
-
-    ! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin)
- 
-  end if
-
-! Compute the RPA correlation energy
-
-  call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI_MO,OmRPA, & 
-                       rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
-
-!--------------!
-! Dump results !
-!--------------!
-
-  call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eG0W0,EcRPA,EcGM)
-
-! Deallocate memory
-
-  deallocate(SigC,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,eG0W0lin)
-
-! Plot stuff
-
-! call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,OmRPA,rho_RPA)
-
-! Perform BSE calculation
-
-  if(BSE) then
-
-    call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eG0W0,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)') 'Tr@BSE@G0W0 correlation energy (singlet) =',EcBSE(1)
-    write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 correlation energy (triplet) =',EcBSE(2)
-    write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 correlation energy           =',EcBSE(1) + EcBSE(2)
-    write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 total energy                 =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
-    write(*,*)'-------------------------------------------------------------------------------'
-    write(*,*)
-
-!   Compute the BSE correlation energy via the adiabatic connection 
-
-    if(doACFDT) then
-
-      write(*,*) '--------------------------------------------------------------'
-      write(*,*) ' Adiabatic connection version of BSE@UG0W0 correlation energy '
-      write(*,*) '--------------------------------------------------------------'
-      write(*,*) 
-
-      if(doXBS) then 
-
-        write(*,*) '*** scaled screening version (XBS) ***'
-        write(*,*)
-
-      end if
-
-      call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF,eG0W0,EcAC)
-
-      write(*,*)
-      write(*,*)'-------------------------------------------------------------------------------'
-      write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 correlation energy (singlet) =',EcAC(1)
-      write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 correlation energy (triplet) =',EcAC(2)
-      write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 correlation energy           =',EcAC(1) + EcAC(2)
-      write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 total energy                 =',ENuc + ERHF + EcAC(1) + EcAC(2)
-      write(*,*)'-------------------------------------------------------------------------------'
-      write(*,*)
-
-    end if
-
-  end if
-
-end subroutine G0W0

src/GF/qsGF2.f90

@@ -64,6 +64,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,
   double precision,allocatable  :: c(:,:)
   double precision,allocatable  :: cp(:,:)
   double precision,allocatable  :: eGF2(:)
+  double precision,allocatable  :: eOld(:)
   double precision,allocatable  :: P(:,:)
   double precision,allocatable  :: F(:,:)
   double precision,allocatable  :: Fp(:,:)
@@ -101,7 +102,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,
 
 ! Memory allocation
 
-  allocate(eGF2(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas),        &
+  allocate(eGF2(nBas),eOld(nbas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
            J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas),      & 
            error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
 
@@ -112,6 +113,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,
   ispin           = 1
   Conv            = 1d0
   P(:,:)          = PHF(:,:)
+  eOld(:)         = eHF(:)
   eGF2(:)         = eHF(:)
   c(:,:)          = cHF(:,:)
   F_diis(:,:)     = 0d0
@@ -157,16 +159,15 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,
     ! Compute commutator and convergence criteria
 
     error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
-    Conv = maxval(abs(error))
 
     ! DIIS extrapolation 
 
     n_diis = min(n_diis+1,max_diis)
+    if(abs(rcond) > 1d-7) then
       call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
-
+    else
-!   Reset DIIS if required
+      n_diis = 0
-
+    end if
-    if(abs(rcond) < 1d-15) n_diis = 0
 
     ! Diagonalize Hamiltonian in AO basis
 
@@ -174,6 +175,12 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,
     cp(:,:) = Fp(:,:)
     call diagonalize_matrix(nBas,cp,eGF2)
     c = matmul(X,cp)
+    SigCp = matmul(transpose(c),matmul(SigCp,c))
+
+    ! Save quasiparticles energy for next cycle
+
+    Conv = maxval(abs(eGF2 - eOld))
+    eOld(:) = eGF2(:)
 
     ! Compute new density matrix in the AO basis
 

src/GF/qsUGF2.f90

@@ -75,6 +75,7 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,
   double precision,allocatable  :: F_diis(:,:,:)
   double precision,allocatable  :: c(:,:,:)
   double precision,allocatable  :: cp(:,:,:)
+  double precision,allocatable  :: eOld(:,:)
   double precision,allocatable  :: eGF2(:,:)
   double precision,allocatable  :: P(:,:,:)
   double precision,allocatable  :: F(:,:,:)
@@ -117,9 +118,10 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,
   nS_bb = nS(2)
   nS_sc = nS_aa + nS_bb
 
-  allocate(eGF2(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas,nspin), &
+  allocate(eGF2(nBas,nspin),eOld(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin),F(nBas,nBas,nspin), &
-           J(nBas,nBas,nspin),K(nBas,nBas,nspin),SigC(nBas,nBas,nspin),SigCp(nBas,nBas,nspin),SigCm(nBas,nBas,nspin),         &
+           Fp(nBas,nBas,nspin),J(nBas,nBas,nspin),K(nBas,nBas,nspin),SigC(nBas,nBas,nspin),SigCp(nBas,nBas,nspin),         &
-           Z(nBas,nspin),error(nBas,nBas,nspin),error_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin))
+           SigCm(nBas,nBas,nspin),Z(nBas,nspin),error(nBas,nBas,nspin),error_diis(nBasSq,max_diis,nspin),                  &
+           F_diis(nBasSq,max_diis,nspin))
 
 ! Initialization
   
@@ -205,14 +207,14 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,
     ! DIIS extrapolation
 
     n_diis = min(n_diis+1,max_diis)
+    if(minval(rcond(:)) > 1d-7) then
       do is=1,nspin
         if(nO(is) > 1) call DIIS_extrapolation(rcond(is),nBasSq,nBasSq,n_diis,error_diis(:,1:n_diis,is), &
                                                F_diis(:,1:n_diis,is),error(:,:,is),F(:,:,is))
       end do
-
+    else 
-    ! Reset DIIS if required
+      n_diis = 0
-
+    end if
-    if(minval(rcond(:)) < 1d-15) n_diis = 0
 
     ! Transform Fock matrix in orthogonal basis
 
@@ -233,12 +235,23 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,
       c(:,:,is) = matmul(X(:,:),cp(:,:,is))
     end do
 
+    ! Back-transform self-energy
+
+    do is=1,nspin
+      SigCp(:,:,is) = matmul(transpose(c(:,:,is)),matmul(SigCp(:,:,is),c(:,:,is)))
+    end do
+
     ! Compute density matrix 
 
     do is=1,nspin
       P(:,:,is) = matmul(c(:,1:nO(is),is),transpose(c(:,1:nO(is),is)))
     end do
 
+    ! Save quasiparticles energy for next cycle
+
+    Conv = maxval(abs(eGF2(:,:) - eOld(:,:)))
+    eOld(:,:) = eGF2(:,:)
+
     !------------------------------------------------------------------------
     !   Compute total energy
     !------------------------------------------------------------------------

src/HF/RHF.f90

@@ -127,11 +127,11 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
 !   DIIS extrapolation
 
     n_diis = min(n_diis+1,max_diis)
+    if(abs(rcond) > 1d-7) then
       call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
-
+    else
-!   Reset DIIS if required
+      n_diis = 0
-
+    end if
-    if(abs(rcond) < 1d-15) n_diis = 0
 
 !  Diagonalize Fock matrix
 

src/HF/UHF.f90

@@ -179,14 +179,14 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO
 !   DIIS extrapolation
 
     n_diis = min(n_diis+1,max_diis)
+    if(minval(rcond(:)) > 1d-7) then
       do ispin=1,nspin
         if(nO(ispin) > 1) call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,1:n_diis,ispin), &
                                                 F_diis(:,1:n_diis,ispin),err(:,:,ispin),F(:,:,ispin))
       end do
-
+    else
-!   Reset DIIS if required
+      n_diis = 0
-
+    end if
-    if(minval(rcond(:)) < 1d-15) n_diis = 0
 
 !------------------------------------------------------------------------
 !   Compute UHF energy

src/MBPT/evGW.f90

@@ -195,11 +195,11 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
     else
 
       n_diis = min(n_diis+1,max_diis)
+      if(abs(rcond) > 1d-7) then
         call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGW-eOld,eGW)
-
+      else
-!    Reset DIIS if required
+        n_diis = 0
-
+      endif
-      if(abs(rcond) < 1d-15) n_diis = 0
 
     endif
 

src/MBPT/print_qsGW.f90

@@ -1,5 +1,4 @@
-subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,P,T,V,J,K,F,SigC,Z, & 
+subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,dipole)
-                      ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,dipole)
 
 ! Print one-electron energies and other stuff for qsGW
 
@@ -23,10 +22,8 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,P,T,V,J,K,F,SigC,Z, &
   double precision,intent(in)        :: eHF(nBas)
   double precision,intent(in)        :: eGW(nBas)
   double precision,intent(in)        :: c(nBas)
-  double precision,intent(in)        :: P(nBas,nBas) 
+  double precision,intent(in)        :: SigC(nBas,nBas)
-  double precision,intent(in)        :: T(nBas,nBas),V(nBas,nBas)
+  double precision,intent(in)        :: Z(nBas)
-  double precision,intent(in)        :: J(nBas,nBas),K(nBas,nBas),F(nBas,nBas)
-  double precision,intent(in)        :: Z(nBas),SigC(nBas,nBas)
   double precision,intent(in)        :: dipole(ncart)
 
 ! Local variables
@@ -67,7 +64,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,P,T,V,J,K,F,SigC,Z, &
 
   write(*,*)'-------------------------------------------------------------------------------'
   write(*,'(2X,A10,I3)')   'Iteration ',nSCF
-  write(*,'(2X,A19,F15.5)')'max(|FPS - SPF|) = ',Conv
+  write(*,'(2X,A14,F15.5)')'Convergence = ',Conv
   write(*,*)'-------------------------------------------'
   write(*,'(2X,A30,F15.6,A3)') 'qsGW HOMO      energy:',eGW(HOMO)*HaToeV,' eV'
   write(*,'(2X,A30,F15.6,A3)') 'qsGW LUMO      energy:',eGW(LUMO)*HaToeV,' eV'

src/MBPT/print_qsUGW.f90

@@ -92,7 +92,7 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, &
   write(*,*)'-------------------------------------------------------------------------------& 
               -------------------------------------------------'
   write(*,'(2X,A10,I3)')   'Iteration ',nSCF
-  write(*,'(2X,A19,F15.5)')'max(|FPS - SPF|) = ',Conv
+  write(*,'(2X,A14,F15.5)')'Convergence = ',Conv
   write(*,*)'-------------------------------------------------------------------------------& 
               -------------------------------------------------'
   write(*,'(2X,A30,F15.6,A3)') 'qsUGW HOMO      energy:',maxval(HOMO(:))*HaToeV,' eV'

src/MBPT/qsGW.f90

@@ -69,6 +69,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
   double precision,external     :: trace_matrix
   double precision              :: dipole(ncart)
 
+  logical                       :: print_W = .true.
   double precision,allocatable  :: error_diis(:,:)
   double precision,allocatable  :: F_diis(:,:)
   double precision,allocatable  :: OmRPA(:)
@@ -78,6 +79,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
   double precision,allocatable  :: c(:,:)
   double precision,allocatable  :: cp(:,:)
   double precision,allocatable  :: eGW(:)
+  double precision,allocatable  :: eOld(:)
   double precision,allocatable  :: P(:,:)
   double precision,allocatable  :: F(:,:)
   double precision,allocatable  :: Fp(:,:)
@@ -136,7 +138,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
 
 ! Memory allocation
 
-  allocate(eGW(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas),        &
+  allocate(eGW(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
            J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas),     & 
            OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS),                           &
            error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
@@ -149,6 +151,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
   Conv            = 1d0
   P(:,:)          = PHF(:,:)
   eGW(:)          = eHF(:)
+  eOld(:)         = eHF(:)
   c(:,:)          = cHF(:,:)
   F_diis(:,:)     = 0d0
   error_diis(:,:) = 0d0
@@ -181,6 +184,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
 
       call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO, &
                            OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
+      if(print_W) call print_excitation('RPA@qsGW    ',ispin,nS,OmRPA)
 
     endif
 
@@ -211,21 +215,18 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
 
     F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigCp(:,:)
 
-    call matout(nBas,nBAs,SigCp)
-
     ! Compute commutator and convergence criteria
 
     error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
-    Conv = maxval(abs(error))
 
     ! DIIS extrapolation 
 
     n_diis = min(n_diis+1,max_diis)
+    if(abs(rcond) > 1d-7) then 
       call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
-
+    else
-!   Reset DIIS if required
+      n_diis = 0
-
+    end if
-    if(abs(rcond) < 1d-15) n_diis = 0
 
     ! Diagonalize Hamiltonian in AO basis
 
@@ -233,6 +234,12 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
     cp(:,:) = Fp(:,:)
     call diagonalize_matrix(nBas,cp,eGW)
     c = matmul(X,cp)
+    SigCp = matmul(transpose(c),matmul(SigCp,c))
+
+    ! Save quasiparticles energy for next cycle
+
+    Conv = maxval(abs(eGW - eOld))
+    eOld(:) = eGW(:)
 
     ! Compute new density matrix in the AO basis
 
@@ -265,7 +272,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
     ! Print results
 
     call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
-    call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,P,T,V,J,K,F,SigCp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,dipole)
+    call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigCp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,dipole)
 
   enddo
 !------------------------------------------------------------------------

src/MBPT/qsUGW.f90

@@ -92,6 +92,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
   double precision,allocatable  :: c(:,:,:)
   double precision,allocatable  :: cp(:,:,:)
   double precision,allocatable  :: eGW(:,:)
+  double precision,allocatable  :: eOld(:,:)
   double precision,allocatable  :: P(:,:,:)
   double precision,allocatable  :: F(:,:,:)
   double precision,allocatable  :: Fp(:,:,:)
@@ -154,10 +155,11 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
   nS_bb = nS(2)
   nS_sc = nS_aa + nS_bb
 
-  allocate(eGW(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas,nspin), &
+  allocate(eGW(nBas,nspin),eOld(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin),F(nBas,nBas,nspin), &
-           J(nBas,nBas,nspin),K(nBas,nBas,nspin),SigC(nBas,nBas,nspin),SigCp(nBas,nBas,nspin),SigCm(nBas,nBas,nspin),        &
+           Fp(nBas,nBas,nspin),J(nBas,nBas,nspin),K(nBas,nBas,nspin),SigC(nBas,nBas,nspin),SigCp(nBas,nBas,nspin),        &
-           Z(nBas,nspin),OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin),              &
+           SigCm(nBas,nBas,nspin),Z(nBas,nspin),OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),                   &
-           error(nBas,nBas,nspin),error_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin))
+           rho_RPA(nBas,nBas,nS_sc,nspin),error(nBas,nBas,nspin),error_diis(nBasSq,max_diis,nspin),                       & 
+           F_diis(nBasSq,max_diis,nspin))
 
 ! Initialization
   
@@ -167,6 +169,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
   Conv              = 1d0
   P(:,:,:)          = PHF(:,:,:)
   eGW(:,:)          = eHF(:,:)
+  eOld(:,:)         = eHF(:,:)
   c(:,:,:)          = cHF(:,:,:)
   F_diis(:,:,:)     = 0d0
   error_diis(:,:,:) = 0d0
@@ -268,14 +271,14 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
     ! DIIS extrapolation
 
     n_diis = min(n_diis+1,max_diis)
+    if(minval(rcond(:)) > 1d-7) then
       do is=1,nspin
         if(nO(is) > 1) call DIIS_extrapolation(rcond(is),nBasSq,nBasSq,n_diis,error_diis(:,1:n_diis,is), &
                                                   F_diis(:,1:n_diis,is),error(:,:,is),F(:,:,is))
       end do
-
+    else 
-    ! Reset DIIS if required
+      n_diis = 0
-
+    end if
-    if(minval(rcond(:)) < 1d-15) n_diis = 0
 
     ! Transform Fock matrix in orthogonal basis
 
@@ -296,12 +299,23 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
       c(:,:,is) = matmul(X(:,:),cp(:,:,is))
     end do
 
+    ! Back-transform self-energy
+
+    do is=1,nspin
+      SigCp(:,:,is) = matmul(transpose(c(:,:,is)),matmul(SigCp(:,:,is),c(:,:,is)))
+    end do
+
     ! Compute density matrix 
 
     do is=1,nspin
       P(:,:,is) = matmul(c(:,1:nO(is),is),transpose(c(:,1:nO(is),is)))
     end do
 
+    ! Save quasiparticles energy for next cycle
+
+    Conv = maxval(abs(eGW(:,:) - eOld(:,:)))
+    eOld(:,:) = eGW(:,:)
+
     !------------------------------------------------------------------------
     !   Compute total energy
     !------------------------------------------------------------------------