working on SRG for GW

src/GW/GG0W0.f90

@@ -40,6 +40,7 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 
   logical                       :: print_W = .true.
   logical                       :: dRPA
+  double precision              :: flow
   double precision              :: EcRPA
   double precision              :: EcBSE
   double precision              :: EcGM
@@ -68,7 +69,6 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 ! Initialization
 
   dRPA = .true.
-  EcRPA = 0d0
 
 ! TDA for W
 
@@ -79,6 +79,8 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 
 ! SRG regularization
 
+  flow = 500d0
+
   if(doSRG) then
 
     write(*,*) '*** SRG regularized qsGW scheme ***'
@@ -113,7 +115,7 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 !------------------------!
 
   if(doSRG) then 
-    call GGW_SRG_self_energy_diag(nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
+    call GGW_SRG_self_energy_diag(flow,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
   else
     call GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
   end if
@@ -142,8 +144,6 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 
   end if
 
-! call GW_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eHF,Om,rho)
-
 ! Compute the RPA correlation energy
 
                  call phGLR_A(dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
@@ -174,35 +174,6 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
-!   Compute the BSE correlation energy via the adiabatic connection 
-
-!   if(doACFDT) then
-
-!     write(*,*) '-------------------------------------------------------------'
-!     write(*,*) ' Adiabatic connection version of BSE@G0W0 correlation energy '
-!     write(*,*) '-------------------------------------------------------------'
-!     write(*,*) 
-
-!     if(doXBS) then 
-
-!       write(*,*) '*** scaled screening version (XBS) ***'
-!       write(*,*)
-
-!     end if
-
-!     call GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,EcBSE)
-
-!     write(*,*)
-!     write(*,*)'-------------------------------------------------------------------------------'
-!     write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy (singlet) =',EcBSE(1),' au'
-!     write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy (triplet) =',EcBSE(2),' au'
-!     write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy           =',EcBSE(1) + EcBSE(2),' au'
-!     write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 total energy                 =',ENuc + EGHF + EcBSE(1) + EcBSE(2),' au'
-!     write(*,*)'-------------------------------------------------------------------------------'
-!     write(*,*)
-
-!   end if
-
   end if
 
   if(doppBSE) then

src/GW/GGW_SRG_self_energy.f90

@@ -1,4 +1,4 @@
-subroutine GGW_SRG_self_energy(nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
+subroutine GGW_SRG_self_energy(flow,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
 
 ! Compute correlation part of the self-energy
 
@@ -7,6 +7,7 @@ subroutine GGW_SRG_self_energy(nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
 
 ! Input variables
 
+  double precision,intent(in)   :: flow
   integer,intent(in)            :: nOrb
   integer,intent(in)            :: nC
   integer,intent(in)            :: nO
@@ -35,7 +36,7 @@ subroutine GGW_SRG_self_energy(nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
 ! SRG flow parameter !
 !--------------------!
 
-  s = 500d0
+  s = flow
 
 !--------------------!
 ! SRG-GW self-energy !

src/GW/GGW_SRG_self_energy_diag.f90

@@ -1,4 +1,4 @@
-subroutine GGW_SRG_self_energy_diag(nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
+subroutine GGW_SRG_self_energy_diag(flow,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
 
 ! Compute correlation part of the self-energy
 
@@ -7,6 +7,7 @@ subroutine GGW_SRG_self_energy_diag(nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
 
 ! Input variables
 
+  double precision,intent(in)   :: flow
   integer,intent(in)            :: nOrb
   integer,intent(in)            :: nC
   integer,intent(in)            :: nO
@@ -35,7 +36,7 @@ subroutine GGW_SRG_self_energy_diag(nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
 ! SRG flow parameter !
 !--------------------!
 
-  s = 500d0
+  s = flow
 
 !--------------------!
 ! SRG-GW self-energy !

src/GW/RG0W0.f90

@@ -46,7 +46,7 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
   logical                       :: plot_self = .false.
   logical                       :: dRPA_W
   integer                       :: isp_W
-  double precision              :: lambda
+  double precision              :: flow
   double precision              :: EcRPA
   double precision              :: EcBSE(nspin)
   double precision              :: EcGM
@@ -73,10 +73,6 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
   write(*,*)'*******************************'
   write(*,*)
 
-! Initialization
-
-  lambda = 1d0
-
 ! Spin manifold and TDA for dynamical screening
 
   isp_W = 1
@@ -89,6 +85,8 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 
 ! SRG regularization
 
+  flow = 500d0
+
   if(doSRG) then
 
     write(*,*) '*** SRG regularized qsGW scheme ***'
@@ -105,8 +103,8 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 ! Compute screening !
 !-------------------!
 
-                 call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,lambda,eHF,ERI,Aph)
+                 call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
-  if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,lambda,ERI,Bph)
+  if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
   call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 
@@ -123,7 +121,7 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 !------------------------!
 
   if(doSRG) then 
-    call RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
+    call RGW_SRG_self_energy_diag(flow,nBas,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
   else
     call RGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
   end if
@@ -162,8 +160,8 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 
 ! Compute the RPA correlation energy
 
-                 call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,lambda,eGW,ERI,Aph)
+                 call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
-  if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,lambda,ERI,Bph)
+  if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
   call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 

src/GW/RGW_SRG_ReSigC.f90

@@ -0,0 +1,52 @@
+double precision function RGW_SRG_ReSigC(p,w,s,nBas,nC,nO,nV,nR,nS,e,Om,rho)
+
+! Compute diagonal of the correlation part of the self-energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: p
+  double precision,intent(in)   :: w
+  double precision,intent(in)   :: s
+  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)   :: e(nBas)
+  double precision,intent(in)   :: Om(nS)
+  double precision,intent(in)   :: rho(nBas,nBas,nS)
+
+! Local variables
+
+  integer                       :: i,a,m
+  double precision              :: Dpim,Dpam
+
+! Initialize 
+
+  RGW_SRG_ReSigC = 0d0
+
+! Occupied part of the correlation self-energy
+
+  do i=nC+1,nO
+    do m=1,nS
+      Dpim = w - e(i) + Om(m)
+      RGW_SRG_ReSigC = RGW_SRG_ReSigC &
+                     + 2d0*rho(p,i,m)**2*(1d0-exp(-2d0*s*Dpim*Dpim))/Dpim
+    end do
+  end do
+
+! Virtual part of the correlation self-energy
+
+  do a=nO+1,nBas-nR
+    do m=1,nS
+      Dpam = w - e(a) - Om(m)
+      RGW_SRG_ReSigC = RGW_SRG_ReSigC & 
+                     + 2d0*rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam
+    end do
+  end do
+
+end function 

src/GW/RGW_SRG_self_energy.f90

@@ -1,4 +1,4 @@
-subroutine RGW_SRG_self_energy(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
+subroutine RGW_SRG_self_energy(flow,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
 
 ! Compute correlation part of the self-energy
 
@@ -7,6 +7,7 @@ subroutine RGW_SRG_self_energy(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
 
 ! Input variables
 
+  double precision,intent(in)   :: flow
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nOrb
   integer,intent(in)            :: nC
@@ -36,7 +37,7 @@ subroutine RGW_SRG_self_energy(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
 ! SRG flow parameter !
 !--------------------!
 
-  s = 500d0
+  s = flow
 
 !--------------------!
 ! SRG-GW self-energy !

src/GW/RGW_SRG_self_energy_diag.f90

@@ -1,4 +1,4 @@
-subroutine RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
+subroutine RGW_SRG_self_energy_diag(flow,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
 
 ! Compute correlation part of the self-energy
 
@@ -7,6 +7,7 @@ subroutine RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,
 
 ! Input variables
 
+  double precision,intent(in)   :: flow
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nOrb
   integer,intent(in)            :: nC
@@ -36,7 +37,7 @@ subroutine RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,
 ! SRG flow parameter !
 !--------------------!
 
-  s = 500d0
+  s = flow
 
 !--------------------!
 ! SRG-GW self-energy !

src/GW/RGW_phACFDT.f90

@@ -62,9 +62,9 @@ subroutine RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas,
 
 ! Hello World
 
-  write(*,*) '-------------------------------------------------------------'
+  write(*,*) '-----------------------------------------------------------'
-  write(*,*) ' Adiabatic connection version of BSE@G0W0 correlation energy '
+  write(*,*) ' Adiabatic connection version of BSE@GW correlation energy '
-  write(*,*) '-------------------------------------------------------------'
+  write(*,*) '-----------------------------------------------------------'
   write(*,*)
 
 ! eXtended BSE

src/GW/RGW_phBSE.f90

@@ -190,5 +190,4 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
 
   end if
 
-
 end subroutine 

src/GW/UG0W0.f90

@@ -49,6 +49,7 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
   logical                       :: dRPA
   integer                       :: is
   integer                       :: ispin
+  double precision              :: flow
   double precision              :: EcRPA
   double precision              :: EcGM(nspin)
   double precision              :: EcBSE(nspin)
@@ -88,6 +89,8 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
 
 ! SRG regularization
 
+  flow = 500d0
+
   if(doSRG) then
 
     write(*,*) '*** SRG regularized qsGW scheme ***'
@@ -130,7 +133,7 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
 !------------------------------------------------!
 
   if(doSRG) then
-    call UGW_SRG_self_energy_diag(nBas,nC,nO,nV,nR,nSt,eHF,Om,rho,SigC,Z,EcGM)
+    call UGW_SRG_self_energy_diag(flow,nBas,nC,nO,nV,nR,nSt,eHF,Om,rho,SigC,Z,EcGM)
   else
     call UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nSt,eHF,Om,rho,SigC,Z,EcGM)
   end if
@@ -186,20 +189,9 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
 
   if(dophBSE) then
 
-    call UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,S, &
+    call UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,S, &
                    ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF,eGW,EcBSE)
 
-    if(exchange_kernel) then
- 
-      EcBSE(1) = 0.5d0*EcBSE(1)
-      EcBSE(2) = 0.5d0*EcBSE(2)
- 
-    else
- 
-      EcBSE(2) = 0.0d0
-
-    end if
-
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy (spin-conserved) = ',EcBSE(1),' au'
@@ -213,18 +205,6 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
 
     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 UGW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,dophBSE,spin_conserved,spin_flip,eta, & 
                        nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGW,EcBSE)
 

src/GW/UGW_SRG_self_energy.f90

@@ -1,4 +1,4 @@
-subroutine UGW_SRG_self_energy(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
+subroutine UGW_SRG_self_energy(flow,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
 
 ! Compute correlation part of the self-energy
 
@@ -7,6 +7,7 @@ subroutine UGW_SRG_self_energy(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
 
 ! Input variables
 
+  double precision,intent(in)   :: flow
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC(nspin)
   integer,intent(in)            :: nO(nspin)
@@ -34,7 +35,7 @@ subroutine UGW_SRG_self_energy(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
 
 ! SRG flow parameter 
 
-  s = 500d0
+  s = flow
 
 ! Initialize 
 

src/GW/UGW_SRG_self_energy_diag.f90

@@ -1,4 +1,4 @@
-subroutine UGW_SRG_self_energy_diag(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
+subroutine UGW_SRG_self_energy_diag(flow,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
 
 ! Compute correlation part of the self-energy
 
@@ -7,6 +7,7 @@ subroutine UGW_SRG_self_energy_diag(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
 
 ! Input variables
 
+  double precision,intent(in)   :: flow
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC(nspin)
   integer,intent(in)            :: nO(nspin)
@@ -34,7 +35,7 @@ subroutine UGW_SRG_self_energy_diag(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
 
 ! SRG flow parameter 
 
-  s = 500d0
+  s = flow
 
 ! Initialize 
 

src/GW/UGW_phACFDT.f90

@@ -63,6 +63,23 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
 
   allocate(Ec(nAC,nspin))
 
+
+! Hello World
+
+  write(*,*) '-----------------------------------------------------------'
+  write(*,*) ' Adiabatic connection version of BSE@GW correlation energy '
+  write(*,*) '-----------------------------------------------------------'
+  write(*,*)
+
+! eXtended BSE
+
+  if(doXBS) then
+
+    write(*,*) '*** scaled screening version (XBS) ***'
+    write(*,*)
+
+  end if
+
 ! Antisymmetrized kernel version
 
   if(exchange_kernel) then

src/GW/UGW_phBSE.f90

@@ -1,4 +1,4 @@
-subroutine UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, &
+subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, &
                      S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cW,eW,eGW,EcBSE)
 
 ! Compute the Bethe-Salpeter excitation energies
@@ -8,6 +8,7 @@ subroutine UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO
 
 ! Input variables
 
+  logical,intent(in)            :: exchange_kernel
   logical,intent(in)            :: TDA_W
   logical,intent(in)            :: TDA
   logical,intent(in)            :: dBSE
@@ -161,4 +162,18 @@ subroutine UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO
  
   end if
 
+  
+! Scale properly correlation energy if exchange is included in interaction kernel
+
+  if(exchange_kernel) then
+
+    EcBSE(1) = 0.5d0*EcBSE(1)
+    EcBSE(2) = 0.5d0*EcBSE(2)
+
+  else
+
+    EcBSE(2) = 0.0d0
+
+  end if
+
 end subroutine 

src/GW/evGGW.f90

@@ -41,10 +41,10 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 
 ! Local variables
 
-  logical                       :: linear_mixing
   logical                       :: dRPA = .true.
   integer                       :: nSCF
   integer                       :: n_diis
+  double precision              :: flow
   double precision              :: rcond
   double precision              :: Conv
   double precision              :: EcRPA
@@ -81,6 +81,8 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 
 ! SRG regularization
 
+  flow = 500d0
+
   if(doSRG) then
 
     write(*,*) '*** SRG regularized qsGW scheme ***'
@@ -88,11 +90,6 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 
   end if
 
-! Linear mixing
-
-  linear_mixing = .false.
-  alpha = 0.2d0
-
 ! Memory allocation
 
   allocate(Aph(nS,nS),Bph(nS,nS),eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas), & 
@@ -130,7 +127,7 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
     ! Compute correlation part of the self-energy 
 
     if(doSRG) then
-      call GGW_SRG_self_energy_diag(nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
+      call GGW_SRG_self_energy_diag(flow,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
     else
       call GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
     end if
@@ -163,18 +160,10 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 
     ! Linear mixing or DIIS extrapolation
 
-    if(linear_mixing) then
+    if(max_diis > 1) then
- 
-      eGW(:) = alpha*eGW(:) + (1d0 - alpha)*eOld(:)
- 
-    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)
-        call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGW-eOld,eGW)
-      else
-        n_diis = 0
-      end if
 
     end if
 
@@ -222,52 +211,8 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
-!   Compute the BSE correlation energy via the adiabatic connection 
-
-!   if(doACFDT) then
-
-!     write(*,*) '------------------------------------------------------'
-!     write(*,*) 'Adiabatic connection version of BSE correlation energy'
-!     write(*,*) '------------------------------------------------------'
-!     write(*,*)
-
-!     if(doXBS) then
-
-!       write(*,*) '*** scaled screening version (XBS) ***'
-!       write(*,*)
-
-!     end if
-
-!     call GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcBSE)
-
-!     write(*,*)
-!     write(*,*)'-------------------------------------------------------------------------------'
-!     write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (singlet) =',EcBSE(1)
-!     write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (triplet) =',EcBSE(2)
-!     write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy           =',EcBSE(1) + EcBSE(2)
-!     write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW total energy                 =',ENuc + EGHF + EcBSE(1) + EcBSE(2)
-!     write(*,*)'-------------------------------------------------------------------------------'
-!     write(*,*)
-
-!   end if
-
   end if
 
-! if(doppBSE) then
-
-!   call GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
-
-!   write(*,*)
-!   write(*,*)'-------------------------------------------------------------------------------'
-!   write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy (singlet) =',EcBSE(1)
-!   write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy (triplet) =',3d0*EcBSE(2)
-!   write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy =',EcBSE(1) + 3d0*EcBSE(2)
-!   write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW total energy =',ENuc + EGHF + EcBSE(1) + 3d0*EcBSE(2)
-!   write(*,*)'-------------------------------------------------------------------------------'
-!   write(*,*)
-
-! end if
-
 ! Testing zone
 
   if(dotest) then

src/GW/evRGW.f90

@@ -44,11 +44,11 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 
 ! Local variables
 
-  logical                       :: linear_mixing
   logical                       :: dRPA = .true.
   integer                       :: ispin
   integer                       :: nSCF
   integer                       :: n_diis
+  double precision              :: flow
   double precision              :: rcond
   double precision              :: Conv
   double precision              :: EcRPA
@@ -85,6 +85,8 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 
 ! SRG regularization
 
+  flow = 500d0
+
   if(doSRG) then
 
     write(*,*) '*** SRG regularized qsGW scheme ***'
@@ -92,11 +94,6 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 
   end if
 
-! Linear mixing
-
-  linear_mixing = .false.
-  alpha = 0.2d0
-
 ! Memory allocation
 
   allocate(Aph(nS,nS),Bph(nS,nS),eGW(nOrb),eOld(nOrb),Z(nOrb),SigC(nOrb), & 
@@ -135,7 +132,7 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
     ! Compute correlation part of the self-energy 
 
     if(doSRG) then
-      call RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
+      call RGW_SRG_self_energy_diag(flow,nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
     else
       call RGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
     end if
@@ -166,20 +163,12 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 
     call print_evRGW(nOrb,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM)
 
-    ! Linear mixing or DIIS extrapolation
+    ! DIIS extrapolation
 
-    if(linear_mixing) then
+    if(max_diis > 1) then
- 
-      eGW(:) = alpha*eGW(:) + (1d0 - alpha)*eOld(:)
- 
-    else
 
       n_diis = min(n_diis+1,max_diis)
-      if(abs(rcond) > 1d-7) then
+      call DIIS_extrapolation(rcond,nOrb,nOrb,n_diis,error_diis,e_diis,eGW-eOld,eGW)
-        call DIIS_extrapolation(rcond,nOrb,nOrb,n_diis,error_diis,e_diis,eGW-eOld,eGW)
-      else
-        n_diis = 0
-      end if
 
     end if
 

src/GW/evUGW.f90

@@ -49,11 +49,11 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
 ! Local variables
 
   logical                       :: dRPA
-  logical                       :: linear_mixing
   integer                       :: is
   integer                       :: ispin
   integer                       :: nSCF
   integer                       :: n_diis
+  double precision              :: flow
   double precision              :: rcond(nspin)
   double precision              :: Conv
   double precision              :: EcRPA(nspin)
@@ -92,6 +92,8 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
 
 ! SRG regularization
 
+  flow = 500d0
+
   if(doSRG) then
 
     write(*,*) '*** SRG regularized qsGW scheme ***'
@@ -104,11 +106,6 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
   EcRPA(:) = 0d0
   dRPA = .true.
 
-! Linear mixing
-
-  linear_mixing = .false.
-  alpha = 0.2d0
-
 ! Memory allocation
 
   nSa = nS(1)
@@ -156,7 +153,7 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
     !------------------------------------------------!
 
     if(doSRG) then
-      call UGW_SRG_self_energy_diag(nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,SigC,Z,EcGM)
+      call UGW_SRG_self_energy_diag(flow,nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,SigC,Z,EcGM)
     else
       call UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,SigC,Z,EcGM)
     end if
@@ -200,11 +197,7 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
 
     ! Linear mixing or DIIS extrapolation
 
-    if(linear_mixing) then
+    if(max_diis > 1) then
- 
-      eGW(:,:) = alpha*eGW(:,:) + (1d0 - alpha)*eOld(:,:)
- 
-    else
 
       n_diis = min(n_diis+1,max_diis)
       do is=1,nspin
@@ -212,10 +205,6 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
                                 e_diis(:,1:n_diis,is),eGW(:,is)-eOld(:,is),eGW(:,is))
       end do
 
-!    Reset DIIS if required
-
-      if(minval(rcond(:)) < 1d-15) n_diis = 0
-
     end if
 
     ! Save quasiparticles energy for next cycle
@@ -253,20 +242,9 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
 
   if(BSE) then
 
-    call UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, &
+    call UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, &
                    S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eGW,eGW,EcBSE)
 
-    if(exchange_kernel) then
-
-      EcBSE(1) = 0.5d0*EcBSE(1)
-      EcBSE(2) = 0.5d0*EcBSE(2)
-
-    else
-
-      EcBSE(2) = 0.0d0
-
-    end if
-
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@evGW@UHF correlation energy (spin-conserved) =',EcBSE(1),' au'
@@ -280,18 +258,6 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
 
     if(doACFDT) then
 
-      write(*,*) '--------------------------------------------------------------'
-      write(*,*) ' Adiabatic connection version of BSE@evUGW correlation energy '
-      write(*,*) '--------------------------------------------------------------'
-      write(*,*)
-
-      if(doXBS) then
-
-        write(*,*) '*** scaled screening version (XBS) ***'
-        write(*,*)
-
-      end if
-
       call UGW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,spin_conserved,spin_flip, &
                        eta,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eGW,eGW,EcRPA)
 

src/GW/qsGGW.f90

@@ -60,6 +60,7 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
   integer                       :: nBas2Sq
   integer                       :: ixyz
   integer                       :: n_diis
+  double precision              :: flow
   double precision              :: ET,ETaa,ETbb
   double precision              :: EV,EVaa,EVbb
   double precision              :: EJ,EJaaaa,EJaabb,EJbbaa,EJbbbb
@@ -132,6 +133,8 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 
 ! SRG regularization
 
+  flow = 500d0
+
   if(doSRG) then
 
     write(*,*) '*** SRG regularized qsGW scheme ***'
@@ -264,7 +267,7 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
     call GGW_excitation_density(nBas2,nC,nO,nR,nS,ERI_MO,XpY,rho)
 
     if(doSRG) then
-      call GGW_SRG_self_energy(nBas2,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
+      call GGW_SRG_self_energy(flow,nBas2,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
     else
       call GGW_self_energy(eta,nBas2,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
     end if
@@ -395,52 +398,8 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
-!   Compute the BSE correlation energy via the adiabatic connection 
-
-!   if(doACFDT) then
-
-!     write(*,*) '------------------------------------------------------'
-!     write(*,*) 'Adiabatic connection version of BSE correlation energy'
-!     write(*,*) '------------------------------------------------------'
-!     write(*,*)
-
-!     if(doXBS) then
-
-!       write(*,*) '*** scaled screening version (XBS) ***'
-!       write(*,*)
-
-!     end if
-
-!     call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcBSE)
-
-!     write(*,*)
-!     write(*,*)'-------------------------------------------------------------------------------'
-!     write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy (singlet) =',EcBSE(1)
-!     write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy (triplet) =',EcBSE(2)
-!     write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy           =',EcBSE(1) + EcBSE(2)
-!     write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW total energy                 =',ENuc + EqsGW + EcBSE(1) + EcBSE(2)
-!     write(*,*)'-------------------------------------------------------------------------------'
-!     write(*,*)
-
-!   end if
-
   end if
 
-! if(doppBSE) then
-!     
-!   call GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eHF,eGW,EcBSE)
-!   
-!   write(*,*)
-!   write(*,*)'-------------------------------------------------------------------------------'
-!   write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@qsGW correlation energy (singlet) =',EcBSE(1)
-!   write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@qsGW correlation energy (triplet) =',3d0*EcBSE(2)
-!   write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@qsGW correlation energy =',EcBSE(1) + 3d0*EcBSE(2)
-!   write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@qsGW total energy =',ENuc + EGHF + EcBSE(1) + 3d0*EcBSE(2)
-!   write(*,*)'-------------------------------------------------------------------------------'
-!   write(*,*)
-
-! end if
-
 ! Testing zone
 
   if(dotest) then

src/GW/qsRGW.f90

@@ -76,6 +76,8 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
   double precision,external     :: trace_matrix
   double precision              :: dipole(ncart)
 
+  double precision              :: flow
+
   logical                       :: dRPA_W  = .true.
   logical                       :: print_W = .false.
   double precision,allocatable  :: err_diis(:,:)
@@ -125,6 +127,8 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 
 ! SRG regularization
 
+  flow = 500d0
+
   if(doSRG) then
 
     write(*,*) '*** SRG regularized qsGW scheme ***'
@@ -207,7 +211,7 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
     call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI_MO,XpY,rho)
 
     if(doSRG) then 
-      call RGW_SRG_self_energy(nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
+      call RGW_SRG_self_energy(flow,nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
     else
       call RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
     end if

src/GW/qsUGW.f90

@@ -1,4 +1,4 @@
-subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
+subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
                 eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, & 
                 dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF)
 
@@ -17,7 +17,7 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
   logical,intent(in)            :: doACFDT
   logical,intent(in)            :: exchange_kernel
   logical,intent(in)            :: doXBS
-  logical,intent(in)            :: BSE
+  logical,intent(in)            :: dophBSE
   logical,intent(in)            :: TDA_W
   logical,intent(in)            :: TDA
   logical,intent(in)            :: dBSE
@@ -66,6 +66,7 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
   integer                       :: is
   integer                       :: n_diis
   integer                       :: nSa,nSb,nSt
+  double precision              :: flow
   double precision              :: dipole(ncart)
 
   double precision              :: ET(nspin)
@@ -128,6 +129,8 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
 
 ! SRG regularization
 
+  flow = 500d0
+
   if(doSRG) then
 
     write(*,*) '*** SRG regularized qsGW scheme ***'
@@ -220,7 +223,7 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
     !------------------------------------------------!
 
     if(doSRG) then
-      call UGW_SRG_self_energy(nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,SigC,Z,EcGM)
+      call UGW_SRG_self_energy(flow,nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,SigC,Z,EcGM)
     else
       call UGW_self_energy(eta,nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,SigC,Z,EcGM)
     end if
@@ -358,22 +361,11 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
 
 ! Perform BSE calculation
 
-  if(BSE) then
+  if(dophBSE) then
 
-    call UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, &
+    call UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, &
                    S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,c,eGW,eGW,EcBSE)
 
-    if(exchange_kernel) then
-
-      EcBSE(1) = 0.5d0*EcBSE(1)
-      EcBSE(2) = 0.5d0*EcBSE(2)
-
-    else
-
-      EcBSE(2) = 0.0d0
-
-    end if
-
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@qsGW@UHF correlation energy (spin-conserved) = ',EcBSE(1),' au'
@@ -387,19 +379,7 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
 
     if(doACFDT) then
 
-      write(*,*) '--------------------------------------------------------------'
+      call UGW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,dophBSE,spin_conserved,spin_flip, &
-      write(*,*) ' Adiabatic connection version of BSE@qsUGW correlation energy '
-      write(*,*) '--------------------------------------------------------------'
-      write(*,*)
-
-      if(doXBS) then
-
-        write(*,*) '*** scaled screening version (XBS) ***'
-        write(*,*)
-
-      end if
-
-      call UGW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,spin_conserved,spin_flip, &
                        eta,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eGW,eGW,EcRPA)
 
       write(*,*)