added cG0F2

src/GF/RGF.f90

@@ -1,7 +1,8 @@
-subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,maxSCF,    &
+subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2,          &
+               renorm,maxSCF,                                                           &
                thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize, &
                eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,        &
-               S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
+               S,X,T,V,Hc,ERI_AO,ERI_MO,CAP,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
 
 ! Green's function module
 
@@ -13,6 +14,7 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max
   logical,intent(in)            :: dotest
 
   logical,intent(in)            :: doG0F2
+  logical,intent(in)            :: docG0F2
   logical,intent(in)            :: doevGF2
   logical,intent(in)            :: doqsGF2
   logical,intent(in)            :: doufG0F02
@@ -52,6 +54,7 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max
   double precision,intent(in)   :: cHF(nBas,nOrb)
   double precision,intent(in)   :: PHF(nBas,nBas)
   double precision,intent(in)   :: S(nBas,nBas)
+  double precision,intent(in)   :: CAP(nBas,nBas)
   double precision,intent(in)   :: T(nBas,nBas)
   double precision,intent(in)   :: V(nBas,nBas)
   double precision,intent(in)   :: Hc(nBas,nBas)
@@ -168,4 +171,21 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max
 
   end if
 
+!------------------------------------------------------------------------
+! Compute complex G0F2 electronic binding energies
+!------------------------------------------------------------------------
+
+  if(doG0F2) then
+
+    call wall_time(start_GF)
+    call cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet, &
+               linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,    &
+               ENuc,ERHF,ERI_MO,CAP,dipole_int_MO,eHF)
+    call wall_time(end_GF)
+
+    t_GF = end_GF - start_GF
+    write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF2 = ',t_GF,' seconds'
+    write(*,*)
+
+  end if
 end subroutine

src/GW/cRG0W0.f90

@@ -135,10 +135,7 @@ subroutine cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TD
 
   ! Linearized or graphical solution?
   Re_eGWlin(:) = eHF(:) + Re_Z(:)*Re_SigC(:) - Im_Z(:)*Im_SigC(:)
-  Im_eGWlin(:) = Re_Z(:)*Im_SigC(:) + Im_Z(:)*Re_SigC(:)
+  Im_eGWlin(:) = e_cap(:) + Re_Z(:)*Im_SigC(:) + Im_Z(:)*Re_SigC(:)
-  do p = 1, nOrb
-        Im_eGWlin(p) = Im_eGWlin(p) + CAP(p,p)  
-  end do
 
   if(linearize) then 
  

src/GW/cRGW_self_energy_diag.f90

@@ -25,8 +25,8 @@ subroutine cRGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Re_S
   integer                       :: i,a,p,m
   double precision              :: num,eps
   double precision              :: eta_tilde
-  double precision              :: Re_DS(nBas)
+  double precision,allocatable  :: Re_DS(:)
-  double precision              :: Im_DS(nBas)
+  double precision,allocatable  :: Im_DS(:)
   
 ! Output variables
 
@@ -37,6 +37,7 @@ subroutine cRGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Re_S
   double precision,intent(out)  :: EcGM
 
 ! Initialize 
+  allocate(Re_DS(nBas),Im_DS(nBas))
   Re_Sig(:) = 0d0
   Im_Sig(:) = 0d0
   Re_DS(:)   = 0d0
@@ -97,5 +98,5 @@ subroutine cRGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Re_S
 ! Compute renormalization factor from derivative 
   Re_Z(:) = (1d0-Re_DS(:))/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
   Im_Z(:) = Im_DS(:)/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
-
+  deallocate(Re_DS,Im_DS)
 end subroutine 

src/HF/cRHF.f90

@@ -146,7 +146,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,ENuc, &
 
     ! Total energy
 
-    ERHF = ET + EV+ EW + EJ + EK
+    ERHF = ET + EV + EW + EJ + EK
 
     ! DIIS extrapolation  !
 
@@ -154,7 +154,6 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,ENuc, &
 
       n_diis = min(n_diis+1,max_diis)
       call complex_DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis,F_diis,err,F)
-      if (rcond<1.0d-10) write(*,*) "!!! DIIS system ill conditioned: rcond = ", rcond," !!!"
     end if
     ! Level shift
     if(level_shift > 0d0 .and. Conv > thresh) call complex_level_shifting(level_shift,nBas,nBas,nO,S,c,F)

src/HF/print_cRHF.f90

@@ -46,7 +46,7 @@ subroutine print_cRHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EW,EJ, EK, ERHF)
   write(*,'(A69)')           '---------------------------------------------------------'
   write(*,'(A33)')           ' Summary               '
   write(*,'(A69)')           '---------------------------------------------------------'
-  write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' One-electron energy = ',real(ET + EV),'+',aimag(ET+EV),'i',' au'
+  write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' One-electron energy = ',real(ET + EV +EW),'+',aimag(ET+EV+EW),'i',' au'
   write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' Kinetic      energy = ',real(ET),'+',aimag(ET),'i',' au'
   write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' Potential    energy = ',real(EV),'+',aimag(EV),'i',' au'
   write(*,'(A33,1X,F16.10,1X,A1,F16.10,A1,A3)') ' CAP          energy = ',real(EW),'+',aimag(EW),'i',' au'
@@ -81,5 +81,10 @@ subroutine print_cRHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EW,EJ, EK, ERHF)
   write(*,'(A37)') '-------------------------------'
   call complex_vecout(nOrb, eHF)
   write(*,*)
+  write(*,'(A37)') '-------------------------------'
+  write(*,'(A37)') ' cRHF orbital energies (eV) '
+  write(*,'(A37)') '-------------------------------'
+  call complex_vecout(nOrb, eHF*HaToeV)
+  write(*,*)
 
 end subroutine 

src/QuAcK/QuAcK.f90

@@ -13,7 +13,7 @@ program QuAcK
   logical                       :: dophRPA,dophRPAx,docrRPA,doppRPA
   logical                       :: doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3
   logical                       :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW
-  logical                       :: docG0W0
+  logical                       :: docG0W0,docG0F2
   logical                       :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh
 
   integer                       :: nNuc
@@ -127,7 +127,7 @@ program QuAcK
                     doG0W0,doevGW,doqsGW,doufG0W0,doufGW,  &
                     doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp, &
                     doG0T0eh,doevGTeh,doqsGTeh,            &
-                    docG0W0,                               &
+                    docG0W0,docG0F2,                       &
                     doRtest,doUtest,doGtest)
 
 !--------------------------!
@@ -258,7 +258,7 @@ program QuAcK
                   dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA,        &
                   doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,                   &
                   doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh,                                       &
-                  docG0W0,                                                                                                &
+                  docG0W0,docG0F2,                                                                                        &
                   nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc,                                                              &
                   S,T,V,Hc,CAP,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift,                       &
                   guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA,                     &

src/QuAcK/RQuAcK.f90

@@ -3,7 +3,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
                   dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA,       & 
                   doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,                  &
                   doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh,                                      & 
-                  docG0W0,                                                                                               & 
+                  docG0W0,docG0F2,                                                                                       & 
                   nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc,                                                             &
                   S,T,V,Hc,CAP_AO,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift,                   &
                   guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,singlet,triplet,TDA,                             &
@@ -34,7 +34,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
   logical,intent(in)            :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW
   logical,intent(in)            :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp
   logical,intent(in)            :: doG0T0eh,doevGTeh,doqsGTeh
-  logical,intent(in)            :: docG0W0
+  logical,intent(in)            :: docG0W0,docG0F2
 
   integer,intent(in)            :: nNuc,nBas,nOrb
   integer,intent(in)            :: nC
@@ -207,7 +207,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
   ! 4-index transform 
   
   call AOtoMO_ERI_RHF(nBas,nOrb,cHF,ERI_AO,ERI_MO)
-  if (docG0W0) call AOtoMO(nBas,nOrb,cHF,CAP_AO,CAP_MO) ! Add the different cases for cGW methods when implemented
+  if (docG0W0 .or. docG0F2) call AOtoMO(nBas,nOrb,cHF,CAP_AO,CAP_MO) ! Add the different cases for cGW methods when implemented
   call wall_time(end_AOtoMO)
 
   t_AOtoMO = end_AOtoMO - start_AOtoMO
@@ -327,15 +327,15 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF,
 ! Green's function module !
 !-------------------------!
 
-  doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3
+doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3 .or. docG0F2
 
   if(doGF) then
 
     call wall_time(start_GF)
-    call RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm_GF,maxSCF_GF, &
+    call RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2,renorm_GF,maxSCF_GF, &
              thresh_GF,max_diis_GF,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,lin_GF, &
              eta_GF,reg_GF,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,            &
-             S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
+             S,X,T,V,Hc,ERI_AO,ERI_MO,CAP_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
     call wall_time(end_GF)
 
     t_GF = end_GF - start_GF

src/QuAcK/read_methods.f90

@@ -1,5 +1,5 @@
 subroutine read_methods(working_dir,                           &
-                        doRHF,doUHF,doGHF,doROHF,doHFB,docRHF,  &
+                        doRHF,doUHF,doGHF,doROHF,doHFB,docRHF, &
                         doMP2,doMP3,                           & 
                         doCCD,dopCCD,doDCD,doCCSD,doCCSDT,     & 
                         do_drCCD,do_rCCD,do_crCCD,do_lCCD,     &
@@ -10,7 +10,7 @@ subroutine read_methods(working_dir,                           &
                         doG0W0,doevGW,doqsGW,doufG0W0,doufGW,  & 
                         doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp, &
                         doG0T0eh,doevGTeh,doqsGTeh,            &
-                        docG0W0,                               &
+                        docG0W0,docG0F2,                       &
                         doRtest,doUtest,doGtest)
 
 ! Read desired methods 
@@ -33,7 +33,7 @@ subroutine read_methods(working_dir,                           &
   logical,intent(out)           :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW
   logical,intent(out)           :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp 
   logical,intent(out)           :: doG0T0eh,doevGTeh,doqsGTeh
-  logical,intent(out)           :: docG0W0
+  logical,intent(out)           :: docG0W0,docG0F2
 
   logical,intent(out)           :: doRtest,doUtest,doGtest
 
@@ -205,12 +205,15 @@ subroutine read_methods(working_dir,                           &
       if(ans2 == 'T') doevGTeh = .true.
       if(ans3 == 'T') doqsGTeh = .true.
       
-      ! Read CAP-GW methods
+      ! Read Complex methods
 
       docG0W0 = .false.
+      doG0F2    = .false.
+
       read(1,*)
-      read(1,*) ans1
+      read(1,*) ans1,ans2
       if(ans1 == 'T') docG0W0    = .true.
+      if(ans2 == 'T') docG0F2    = .true.
       
       ! Read test