regularized GW

input/methods

@@ -13,7 +13,7 @@
 # G0F2* evGF2* qsGF2* G0F3 evGF3
   F     F       F      F    F
 # G0W0* evGW* qsGW* ufG0W0 ufGW
-  F    F    T       F      F
+  F    T    F       F      F
 # G0T0 evGT qsGT 
   F    F    F
 # MCMP2

input/options

@@ -1,15 +1,17 @@
 # HF: maxSCF thresh   DIIS n_diis guess_type ortho_type mix_guess stability
-       1024  0.00001  T    5      1          1            T         F
+      1024   0.00001  T    5      1          1          T         F
 # MP: 
       
 # CC: maxSCF thresh   DIIS n_diis
       64     0.00001  T    5
 # spin: TDA singlet triplet spin_conserved spin_flip 
-         F   T       T           T          T 
+         F  T       T       T              T 
-# GF: maxSCF thresh  DIIS n_diis lin eta renorm 
+# GF: maxSCF thresh  DIIS n_diis lin eta renorm reg
-       256  0.00001 T    5      T    0.00367493  3      
+      256  0.00001   T    5      T   0.0 3      F
-# GW/GT: maxSCF thresh  DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 
+# GW: maxSCF thresh  DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 reg
-         256   0.00001   T    5     T   0.0 F      F     F     F   F  
+      256    0.00001  T    5     T   0.0 F      F     F     F   F   T
+# GT: maxSCF thresh  DIIS n_diis lin eta TDA_T reg
+      256    0.00001 T    5      T   0.0 F     F  
 # ACFDT: AC Kx  XBS
          F  F   F
 # BSE: BSE dBSE dTDA evDyn

mol/h2.xyz

@@ -1,4 +1,4 @@
 2
  
 H 0. 0. 0.
-H 0. 0. 0.5
+H 0. 0. 0.7

src/GF/G0F2.f90

@@ -1,4 +1,5 @@
-subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
+subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,regularize, &
+                nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
 
 ! Perform a one-shot second-order Green function calculation
 
@@ -16,6 +17,7 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,nBas,nC,nO
   logical,intent(in)            :: triplet
   logical,intent(in)            :: linearize
   double precision,intent(in)   :: eta
+  logical,intent(in)            :: regularize
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nO
   integer,intent(in)            :: nC

src/GF/UG0F2.f90

@@ -1,4 +1,4 @@
-subroutine UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
+subroutine UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
                  S,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,eHF)
 
 ! Perform unrestricted G0W0 calculation
@@ -18,6 +18,7 @@ subroutine UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linearize,eta,
   logical,intent(in)            :: spin_flip
   logical,intent(in)            :: linearize
   double precision,intent(in)   :: eta
+  logical,intent(in)            :: regularize
 
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC(nspin)

src/GF/evGF2.f90

@@ -1,5 +1,5 @@
 subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet, &
-                 linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
+                 linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
 
 ! Perform eigenvalue self-consistent second-order Green function calculation
 
@@ -20,6 +20,8 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet,
   logical,intent(in)            :: triplet
   logical,intent(in)            :: linearize
   double precision,intent(in)   :: eta
+  logical,intent(in)            :: regularize
+
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nO
   integer,intent(in)            :: nC

src/GF/evUGF2.f90

@@ -1,5 +1,5 @@
 subroutine evUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, &
-                  eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,   & 
+                  eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,   & 
                   dipole_int_aa,dipole_int_bb,cHF,eHF)
 
 ! Perform self-consistent eigenvalue-only GW calculation
@@ -22,6 +22,7 @@ subroutine evUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,
   logical,intent(in)            :: spin_conserved
   logical,intent(in)            :: spin_flip
   double precision,intent(in)   :: eta
+  logical,intent(in)            :: regularize
 
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC(nspin)

src/GF/qsGF2.f90

@@ -1,5 +1,5 @@
 subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, & 
-                eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
+                eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
                 S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
 
 ! Perform a quasiparticle self-consistent GF2 calculation
@@ -20,6 +20,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,
   logical,intent(in)            :: singlet
   logical,intent(in)            :: triplet
   double precision,intent(in)   :: eta
+  logical,intent(in)            :: regularize
 
   integer,intent(in)            :: nNuc
   double precision,intent(in)   :: ZNuc(nNuc)

src/GF/qsUGF2.f90

@@ -1,4 +1,4 @@
-subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta, &
+subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,regularize, &
                  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)
 
@@ -20,6 +20,7 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,
   logical,intent(in)            :: spin_conserved
   logical,intent(in)            :: spin_flip
   double precision,intent(in)   :: eta
+  logical,intent(in)            :: regularize
 
   integer,intent(in)            :: nNuc
   double precision,intent(in)   :: ZNuc(nNuc)

src/GT/G0T0.f90

@@ -1,5 +1,5 @@
 subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,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,eG0T0)
+                linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0T0)
 
 ! Perform one-shot calculation with a T-matrix self-energy (G0T0)
 
@@ -21,6 +21,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
   logical,intent(in)            :: triplet
   logical,intent(in)            :: linearize
   double precision,intent(in)   :: eta
+  logical,intent(in)            :: regularize
 
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC

src/GT/evGT.f90

@@ -1,5 +1,5 @@
 subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & 
-                BSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas, & 
+                BSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,regularize,nBas, & 
                 nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0T0)
 
 ! Perform eigenvalue self-consistent calculation with a T-matrix self-energy (evGT)
@@ -24,6 +24,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
   logical,intent(in)            :: singlet
   logical,intent(in)            :: triplet
   double precision,intent(in)   :: eta
+  logical,intent(in)            :: regularize
 
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC

src/GT/qsGT.f90

@@ -1,5 +1,5 @@
 subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,          & 
-                dBSE,dTDA,evDyn,singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
+                dBSE,dTDA,evDyn,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
                 S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
 
 ! Perform a quasiparticle self-consistent GT calculation
@@ -24,6 +24,7 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T
   logical,intent(in)            :: singlet
   logical,intent(in)            :: triplet
   double precision,intent(in)   :: eta
+  logical,intent(in)            :: regularize
 
   integer,intent(in)            :: nNuc
   double precision,intent(in)   :: ZNuc(nNuc)

src/GT/static_Tmatrix_TA.f90

@@ -50,7 +50,7 @@ subroutine static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r
 
           do kl=1,nOO
 !           chi = chi + lambda*rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2)
-            chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2)
+            chi = chi - rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2)
           enddo
 
           TA(ia,jb) = TA(ia,jb) - 1d0*lambda*chi

src/GT/static_Tmatrix_TB.f90

@@ -50,7 +50,7 @@ subroutine static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r
 
           do kl=1,nOO
 !           chi = chi + lambda*rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2
-            chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2
+            chi = chi - rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2
           enddo
 
           TB(ia,jb) = TB(ia,jb) - 1d0*lambda*chi

src/GW/G0W0.f90

@@ -1,5 +1,5 @@
 subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & 
-                dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,            & 
+                dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,regularize, & 
                 nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0W0)
 
 ! Perform G0W0 calculation
@@ -24,6 +24,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, &
   logical,intent(in)            :: triplet
   logical,intent(in)            :: linearize
   double precision,intent(in)   :: eta
+  logical,intent(in)            :: regularize
 
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC
@@ -124,13 +125,18 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, &
 !------------------------!
 
   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)
 
-!--------------------------------!
+  if(regularize) then 
-! Compute renormalization factor !
-!--------------------------------!
 
-  call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
+    call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
+    call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
+
+  else
+
+    call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
+    call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
+
+  end if
 
 !-----------------------------------!
 ! Solve the quasi-particle equation !

src/GW/UG0W0.f90

@@ -1,5 +1,5 @@
 subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, &
-                 linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb,                & 
+                 linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb,     & 
                  dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eGW)
 
 ! Perform unrestricted G0W0 calculation
@@ -24,6 +24,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
   logical,intent(in)            :: spin_flip
   logical,intent(in)            :: linearize
   double precision,intent(in)   :: eta
+  logical,intent(in)            :: regularize
 
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC(nspin)

src/GW/evGW.f90

@@ -1,5 +1,5 @@
 subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,  & 
-                G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,        &
+                G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,        &
                 ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0W0)
 
 ! Perform self-consistent eigenvalue-only GW calculation
@@ -29,6 +29,8 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
   logical,intent(in)            :: singlet
   logical,intent(in)            :: triplet
   double precision,intent(in)   :: eta
+  logical,intent(in)            :: regularize
+
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC
   integer,intent(in)            :: nO
@@ -151,7 +153,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
       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)
 
-    endif
+    end if
 
 !   Compute spectral weights
 
@@ -161,17 +163,33 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
 
     if(G0W) then
 
-!     call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
+      if(regularize) then 
-      call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
+
-      call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
+        call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
+        call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
+
+      else
+
+        call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
+        call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
+
+      end if
 
     else 
 
-!     call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
+      if(regularize) then 
-      call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
-      call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z)
 
-    endif
+        call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
+        call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z)
+
+      else
+
+        call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
+        call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z)
+
+      end if
+
+    end if
 
     ! Solve the quasi-particle equation
 
@@ -198,9 +216,9 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
         call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGW-eOld,eGW)
       else
         n_diis = 0
-      endif
+      end if
 
-    endif
+    end if
 
     ! Save quasiparticles energy for next cycle
 
@@ -210,7 +228,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
 
     nSCF = nSCF + 1
 
-  enddo
+  end do
 !------------------------------------------------------------------------
 ! End main loop
 !------------------------------------------------------------------------
@@ -231,7 +249,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
 
     stop
 
-  endif
+  end if
 
 ! Deallocate memory
 
@@ -288,6 +306,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
 
     end if
 
-  endif
+  end if
 
 end subroutine evGW

src/GW/evUGW.f90

@@ -1,5 +1,5 @@
 subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,    & 
-                G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc, &
+                G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc, &
                 EUHF,S,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eG0W0)
 
 ! Perform self-consistent eigenvalue-only GW calculation
@@ -29,6 +29,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE
   logical,intent(in)            :: spin_conserved
   logical,intent(in)            :: spin_flip
   double precision,intent(in)   :: eta
+  logical,intent(in)            :: regularize
 
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC(nspin)

src/GW/qsGW.f90

@@ -1,5 +1,5 @@
 subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,          & 
-                G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
+                G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
                 S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
 
 ! Perform a quasiparticle self-consistent GW calculation
@@ -27,6 +27,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
   logical,intent(in)            :: singlet
   logical,intent(in)            :: triplet
   double precision,intent(in)   :: eta
+  logical,intent(in)            :: regularize
 
   integer,intent(in)            :: nNuc
   double precision,intent(in)   :: ZNuc(nNuc)
@@ -192,15 +193,31 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
 
     if(G0W) then
 
-!     call regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
+      if(regularize) then
-      call self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
+
-      call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
+        call regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
+        call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
+
+      else
+
+        call self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
+        call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
+
+      end if
 
      else
 
-!     call regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
+       if(regularize) then
-      call self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
+
-      call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z)
+         call regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
+         call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z)
+
+       else
+
+         call self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
+         call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z)
+
+       endif
 
      endif
 

src/GW/qsUGW.f90

@@ -1,5 +1,5 @@
 subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,    & 
-                G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO, & 
+                G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,regularize,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)
 
@@ -28,6 +28,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE
   logical,intent(in)            :: spin_conserved
   logical,intent(in)            :: spin_flip
   double precision,intent(in)   :: eta
+  logical,intent(in)            :: regularize
 
   integer,intent(in)            :: nNuc
   double precision,intent(in)   :: ZNuc(nNuc)

src/GW/regularized_renormalization_factor.f90

@@ -0,0 +1,87 @@
+subroutine regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z)
+
+! Compute the regularized version of the GW renormalization factor
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  logical,intent(in)            :: COHSEX
+  double precision,intent(in)   :: eta
+  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)   :: Omega(nS)
+  double precision,intent(in)   :: rho(nBas,nBas,nS)
+
+! Local variables
+
+  integer                       :: i,a,p,jb
+  double precision              :: eps
+
+  double precision              :: kappa
+  double precision              :: fk,dfk
+
+! Output variables
+
+  double precision,intent(out)  :: Z(nBas)
+
+! Initialize
+
+  Z(:)  = 0d0
+
+!-----------------------------------------!
+! Parameters for regularized calculations !
+!-----------------------------------------!
+
+  kappa = 1.1d0
+
+! static COHSEX approximation
+
+  if(COHSEX) then
+    
+    Z(:) = 1d0
+    return
+ 
+  else
+
+    ! Occupied part of the correlation self-energy
+
+    do p=nC+1,nBas-nR
+      do i=nC+1,nO
+        do jb=1,nS
+          eps = e(p) - e(i) + Omega(jb) 
+          fk  = (1d0 - exp(-kappa*abs(eps)))**2/eps
+          dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
+          dfk = dfk*fk
+          Z(p) = Z(p)  - 2d0*rho(p,i,jb)**2*dfk
+        end do
+      end do
+    end do
+
+    ! Virtual part of the correlation self-energy
+
+    do p=nC+1,nBas-nR
+      do a=nO+1,nBas-nR
+        do jb=1,nS
+          eps = e(p) - e(a) - Omega(jb) 
+          fk  = (1d0 - exp(-kappa*abs(eps)))**2/eps
+          dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
+          dfk = dfk*fk
+          Z(p) = Z(p)  - 2d0*rho(p,a,jb)**2*(eps/(eps**2 + eta**2))**2
+        end do
+      end do
+    end do
+
+  end if
+
+! Compute renormalization factor from derivative of SigC
+ 
+  Z(:) = 1d0/(1d0 - Z(:))
+
+end subroutine regularized_renormalization_factor

src/GW/regularized_self_energy_correlation.f90

@@ -9,7 +9,12 @@ subroutine regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,
 
   logical,intent(in)            :: COHSEX
   double precision,intent(in)   :: eta
-  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
+  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)   :: Omega(nS)
   double precision,intent(in)   :: rho(nBas,nBas,nS)

src/GW/regularized_self_energy_correlation_diag.f90

@@ -23,7 +23,6 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,
 
   integer                       :: i,a,p,q,jb
   double precision              :: eps
-  double precision,external     :: SigC_dcgw
 
   double precision              :: kappa
   double precision              :: fk
@@ -37,9 +36,9 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,
 
   SigC(:) = 0d0
 
-!---------------------------------------------!
+!-----------------------------------------!
-! Parameters for regularized MP2 calculations !
+! Parameters for regularized calculations !
-!---------------------------------------------!
+!-----------------------------------------!
 
   kappa = 1.1d0
 

src/GW/renormalization_factor.f90

@@ -9,14 +9,19 @@ subroutine renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z)
 
   logical,intent(in)            :: COHSEX
   double precision,intent(in)   :: eta
-  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
+  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)   :: Omega(nS)
   double precision,intent(in)   :: rho(nBas,nBas,nS)
 
 ! Local variables
 
-  integer                       :: i,j,a,b,x,jb
+  integer                       :: p,i,a,jb
   double precision              :: eps
 
 ! Output variables
@@ -38,30 +43,22 @@ subroutine renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z)
 
     ! Occupied part of the correlation self-energy
 
-    do x=nC+1,nBas-nR
+    do p=nC+1,nBas-nR
       do i=nC+1,nO
-        jb = 0
+        do jb=1,nS
-        do j=nC+1,nO
+          eps = e(p) - e(i) + Omega(jb) 
-          do b=nO+1,nBas-nR
+          Z(p) = Z(p)  - 2d0*rho(p,i,jb)**2*(eps/(eps**2 + eta**2))**2
-            jb = jb + 1
-            eps = e(x) - e(i) + Omega(jb) 
-            Z(x) = Z(x)  - 2d0*rho(x,i,jb)**2*(eps/(eps**2 + eta**2))**2
-          end do
         end do
       end do
     end do
 
     ! Virtual part of the correlation self-energy
 
-    do x=nC+1,nBas-nR
+    do p=nC+1,nBas-nR
       do a=nO+1,nBas-nR
-        jb = 0
+        do jb=1,nS
-        do j=nC+1,nO
+          eps = e(p) - e(a) - Omega(jb) 
-          do b=nO+1,nBas-nR
+          Z(p) = Z(p)  - 2d0*rho(p,a,jb)**2*(eps/(eps**2 + eta**2))**2
-            jb = jb + 1
-            eps = e(x) - e(a) - Omega(jb) 
-            Z(x) = Z(x)  - 2d0*rho(x,a,jb)**2*(eps/(eps**2 + eta**2))**2
-          end do
         end do
       end do
     end do

src/GW/self_energy_correlation.f90

@@ -9,7 +9,12 @@ subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Ec
 
   logical,intent(in)            :: COHSEX
   double precision,intent(in)   :: eta
-  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
+  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)   :: Omega(nS)
   double precision,intent(in)   :: rho(nBas,nBas,nS)

src/GW/self_energy_correlation_diag.f90

@@ -23,7 +23,6 @@ subroutine self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,r
 
   integer                       :: i,a,p,q,jb
   double precision              :: eps
-  double precision,external     :: SigC_dcgw
 
 ! Output variables
 

src/GW/ufBSE.f90

@@ -1,4 +1,4 @@
-subroutine ufBSE(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
+subroutine ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
 
 ! Unfold BSE@GW equations
 
@@ -7,7 +7,6 @@ subroutine ufBSE(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
 
 ! Input variables
 
-  double precision,intent(in)   :: eta
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC
   integer,intent(in)            :: nO

src/GW/ufG0W0.f90

@@ -1,4 +1,4 @@
-subroutine ufG0W0(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
+subroutine ufG0W0(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
 
 ! Unfold G0W0 equations
 
@@ -7,7 +7,6 @@ subroutine ufG0W0(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
 
 ! Input variables
 
-  double precision,intent(in)   :: eta
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC
   integer,intent(in)            :: nO

src/GW/ufGW.f90

@@ -1,4 +1,4 @@
-subroutine ufGW(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
+subroutine ufGW(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
 
 ! Unfold GW equations
 
@@ -7,7 +7,6 @@ subroutine ufGW(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
 
 ! Input variables
 
-  double precision,intent(in)   :: eta
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC
   integer,intent(in)            :: nO

src/GW/unrestricted_regularized_renormalization_factor.f90

@@ -0,0 +1,111 @@
+subroutine unrestricted_regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,Z)
+
+! Compute the renormalization factor in the unrestricted formalism
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC(nspin)
+  integer,intent(in)            :: nO(nspin)
+  integer,intent(in)            :: nV(nspin)
+  integer,intent(in)            :: nR(nspin)
+  integer,intent(in)            :: nSt
+  double precision,intent(in)   :: e(nBas,nspin)
+  double precision,intent(in)   :: Omega(nSt)
+  double precision,intent(in)   :: rho(nBas,nBas,nSt,nspin)
+
+! Local variables
+
+  integer                       :: i,a,p,jb
+  double precision              :: eps
+
+  double precision              :: kappa
+  double precision              :: fk,dfk
+
+! Output variables
+
+  double precision,intent(out)  :: Z(nBas,nspin)
+
+! Initialize 
+
+  Z(:,:) = 0d0
+
+!-----------------------------------------!
+! Parameters for regularized calculations !
+!-----------------------------------------!
+
+  kappa = 1.1d0
+
+!--------------!
+! Spin-up part !
+!--------------!
+
+  ! Occupied part of the renormalization factor
+
+  do p=nC(1)+1,nBas-nR(1)
+    do i=nC(1)+1,nO(1)
+      do jb=1,nSt
+        eps = e(p,1) - e(i,1) + Omega(jb)
+        fk  = (1d0 - exp(-kappa*abs(eps)))**2/eps
+        dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
+        dfk = dfk*fk
+        Z(p,1) = Z(p,1) + rho(p,i,jb,1)**2*dfk
+      end do
+    end do
+  end do
+
+  ! Virtual part of the correlation self-energy
+
+  do p=nC(1)+1,nBas-nR(1)
+    do a=nO(1)+1,nBas-nR(1)
+      do jb=1,nSt
+        eps = e(p,1) - e(a,1) - Omega(jb)
+        fk  = (1d0 - exp(-kappa*abs(eps)))**2/eps
+        dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
+        dfk = dfk*fk
+        Z(p,1) = Z(p,1) + rho(p,a,jb,1)**2*dfk
+      end do
+    end do
+  end do
+
+!----------------!
+! Spin-down part !
+!----------------!
+
+  ! Occupied part of the correlation self-energy
+
+  do p=nC(2)+1,nBas-nR(2)
+    do i=nC(2)+1,nO(2)
+      do jb=1,nSt
+        eps = e(p,2) - e(i,2) + Omega(jb)
+        fk  = (1d0 - exp(-kappa*abs(eps)))**2/eps
+        dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
+        dfk = dfk*fk
+        Z(p,2) = Z(p,2) + rho(p,i,jb,2)**2*dfk
+      end do
+    end do
+  end do
+
+  ! Virtual part of the correlation self-energy
+
+  do p=nC(2)+1,nBas-nR(2)
+    do a=nO(2)+1,nBas-nR(2)
+      do jb=1,nSt
+        eps = e(p,2) - e(a,2) - Omega(jb)
+        fk  = (1d0 - exp(-kappa*abs(eps)))**2/eps
+        dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
+        dfk = dfk*fk
+        Z(p,2) = Z(p,2) + rho(p,a,jb,2)**2*dfk
+      end do
+    end do
+  end do
+
+! Final rescaling
+
+  Z(:,:) = 1d0/(1d0 + Z(:,:))
+
+end subroutine unrestricted_regularized_renormalization_factor

src/GW/unrestricted_regularized_self_energy_correlation.f90

@@ -0,0 +1,133 @@
+subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC,EcGM)
+
+! Compute diagonal of the correlation part of the self-energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC(nspin)
+  integer,intent(in)            :: nO(nspin)
+  integer,intent(in)            :: nV(nspin)
+  integer,intent(in)            :: nR(nspin)
+  integer,intent(in)            :: nSt
+  double precision,intent(in)   :: e(nBas,nspin)
+  double precision,intent(in)   :: Omega(nSt)
+  double precision,intent(in)   :: rho(nBas,nBas,nSt,nspin)
+
+! Local variables
+
+  integer                       :: i,a,p,q,jb
+  double precision              :: eps
+
+  double precision              :: kappa
+  double precision              :: fk
+
+! Output variables
+
+  double precision,intent(out)  :: SigC(nBas,nBas,nspin)
+  double precision              :: EcGM(nspin)
+
+! Initialize 
+
+  SigC(:,:,:) = 0d0
+  EcGM(:)   = 0d0
+
+!-----------------------------------------!
+! Parameters for regularized calculations !
+!-----------------------------------------!
+
+  kappa = 1.1d0
+
+!--------------!
+! Spin-up part !
+!--------------!
+
+  ! Occupied part of the correlation self-energy
+
+  do p=nC(1)+1,nBas-nR(1)
+    do q=nC(1)+1,nBas-nR(1)
+      do i=nC(1)+1,nO(1)
+        do jb=1,nSt
+          eps = e(p,1) - e(i,1) + Omega(jb)
+          fk  = (1d0 - exp(-kappa*abs(eps)))**2/eps
+          SigC(p,q,1) = SigC(p,q,1) + rho(p,i,jb,1)*rho(q,i,jb,1)*eps/(eps**2 + eta**2)
+        end do
+      end do
+    end do
+  end do
+
+  ! Virtual part of the correlation self-energy
+
+  do p=nC(1)+1,nBas-nR(1)
+    do q=nC(1)+1,nBas-nR(1)
+      do a=nO(1)+1,nBas-nR(1)
+        do jb=1,nSt
+          eps = e(p,1) - e(a,1) - Omega(jb)
+          fk  = (1d0 - exp(-kappa*abs(eps)))**2/eps
+          SigC(p,q,1) = SigC(p,q,1) + rho(p,a,jb,1)*rho(q,a,jb,1)*eps/(eps**2 + eta**2)
+        end do
+      end do
+    end do
+  end do
+
+  ! GM correlation energy
+
+  do i=nC(1)+1,nO(1)
+    do a=nO(1)+1,nBas-nR(1)
+      do jb=1,nSt
+        eps = e(a,1) - e(i,1) + Omega(jb)
+        fk  = (1d0 - exp(-kappa*abs(eps)))**2/eps
+        EcGM(1) = EcGM(1) - rho(a,i,jb,1)**2*eps/(eps**2 + eta**2)
+      end do
+    end do
+  end do
+
+!----------------!
+! Spin-down part !
+!----------------!
+
+  ! Occupied part of the correlation self-energy
+
+  do p=nC(2)+1,nBas-nR(2)
+    do q=nC(2)+1,nBas-nR(2)
+      do i=nC(2)+1,nO(2)
+        do jb=1,nSt
+          eps = e(p,2) - e(i,2) + Omega(jb)
+          fk  = (1d0 - exp(-kappa*abs(eps)))**2/eps
+          SigC(p,q,2) = SigC(p,q,2) + rho(p,i,jb,2)*rho(q,i,jb,2)*eps/(eps**2 + eta**2)
+        end do
+      end do
+    end do
+  end do
+
+  ! Virtual part of the correlation self-energy
+
+  do p=nC(2)+1,nBas-nR(2)
+    do q=nC(2)+1,nBas-nR(2)
+      do a=nO(2)+1,nBas-nR(2)
+        do jb=1,nSt
+          eps = e(p,2) - e(a,2) - Omega(jb)
+          fk  = (1d0 - exp(-kappa*abs(eps)))**2/eps
+          SigC(p,q,2) = SigC(p,q,2) + rho(p,a,jb,2)*rho(q,a,jb,2)*eps/(eps**2 + eta**2)
+        end do
+      end do
+    end do
+  end do
+
+  ! GM correlation energy
+
+  do i=nC(2)+1,nO(2)
+    do a=nO(2)+1,nBas-nR(2)
+      do jb=1,nSt
+        eps = e(a,2) - e(i,2) + Omega(jb)
+        fk  = (1d0 - exp(-kappa*abs(eps)))**2/eps
+        EcGM(2) = EcGM(2) - rho(a,i,jb,2)**2*eps/(eps**2 + eta**2)
+      end do
+    end do
+  end do
+
+end subroutine unrestricted_self_energy_correlation

src/GW/unrestricted_regularized_self_energy_correlation_diag.f90

@@ -0,0 +1,126 @@
+subroutine unrestricted_regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC,EcGM)
+
+! Compute diagonal of the correlation part of the self-energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC(nspin)
+  integer,intent(in)            :: nO(nspin)
+  integer,intent(in)            :: nV(nspin)
+  integer,intent(in)            :: nR(nspin)
+  integer,intent(in)            :: nSt
+  double precision,intent(in)   :: e(nBas,nspin)
+  double precision,intent(in)   :: Omega(nSt)
+  double precision,intent(in)   :: rho(nBas,nBas,nSt,nspin)
+
+! Local variables
+
+  integer                       :: i,a,p,q,jb
+  double precision              :: eps
+
+  double precision              :: kappa
+  double precision              :: fk
+
+! Output variables
+
+  double precision,intent(out)  :: SigC(nBas,nspin)
+  double precision              :: EcGM(nspin)
+
+! Initialize 
+
+  SigC(:,:) = 0d0
+  EcGM(:)   = 0d0
+
+!-----------------------------------------!
+! Parameters for regularized calculations !
+!-----------------------------------------!
+
+  kappa = 1.1d0
+
+!--------------!
+! Spin-up part !
+!--------------!
+
+  ! Occupied part of the correlation self-energy
+
+  do p=nC(1)+1,nBas-nR(1)
+    do i=nC(1)+1,nO(1)
+      do jb=1,nSt
+        eps = e(p,1) - e(i,1) + Omega(jb)
+        fk  = (1d0 - exp(-kappa*abs(eps)))**2/eps
+        SigC(p,1) = SigC(p,1) + rho(p,i,jb,1)**2*fk
+      end do
+    end do
+  end do
+
+  ! Virtual part of the correlation self-energy
+
+  do p=nC(1)+1,nBas-nR(1)
+    do a=nO(1)+1,nBas-nR(1)
+      do jb=1,nSt
+        eps = e(p,1) - e(a,1) - Omega(jb)
+        fk  = (1d0 - exp(-kappa*abs(eps)))**2/eps
+        SigC(p,1) = SigC(p,1) + rho(p,a,jb,1)**2*fk
+      end do
+    end do
+  end do
+
+  ! GM correlation energy
+
+  do i=nC(1)+1,nO(1)
+    do a=nO(1)+1,nBas-nR(1)
+      do jb=1,nSt
+        eps = e(a,1) - e(i,1) + Omega(jb)
+        fk  = (1d0 - exp(-kappa*abs(eps)))**2/eps
+        EcGM(1) = EcGM(1) - rho(a,i,jb,1)**2*fk
+      end do
+    end do
+  end do
+
+!----------------!
+! Spin-down part !
+!----------------!
+
+  ! Occupied part of the correlation self-energy
+
+  do p=nC(2)+1,nBas-nR(2)
+    do i=nC(2)+1,nO(2)
+      do jb=1,nSt
+        eps = e(p,2) - e(i,2) + Omega(jb)
+        fk  = (1d0 - exp(-kappa*abs(eps)))**2/eps
+        SigC(p,2) = SigC(p,2) + rho(p,i,jb,2)**2*fk
+      end do
+    end do
+  end do
+
+  ! Virtual part of the correlation self-energy
+
+  do p=nC(2)+1,nBas-nR(2)
+    do a=nO(2)+1,nBas-nR(2)
+      do jb=1,nSt
+        eps = e(p,2) - e(a,2) - Omega(jb)
+        fk  = (1d0 - exp(-kappa*abs(eps)))**2/eps
+        SigC(p,2) = SigC(p,2) + rho(p,a,jb,2)**2*fk
+      end do
+    end do
+  end do
+
+  ! GM correlation energy
+
+  do i=nC(2)+1,nO(2)
+    do a=nO(2)+1,nBas-nR(2)
+      do jb=1,nSt
+        eps = e(a,2) - e(i,2) + Omega(jb)
+        fk  = (1d0 - exp(-kappa*abs(eps)))**2/eps
+        EcGM(2) = EcGM(2) - rho(a,i,jb,2)**2*fk
+      end do
+    end do
+  end do
+
+
+end subroutine unrestricted_regularized_self_energy_correlation_diag

src/GW/unrestricted_renormalization_factor.f90

@@ -20,7 +20,7 @@ subroutine unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,
 
 ! Local variables
 
-  integer                       :: i,j,a,b,p,q,jb
+  integer                       :: i,a,p,jb
   double precision              :: eps
 
 ! Output variables

src/GW/unrestricted_self_energy_correlation.f90

@@ -20,7 +20,7 @@ subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega
 
 ! Local variables
 
-  integer                       :: i,j,a,b,p,q,jb
+  integer                       :: i,a,p,q,jb
   double precision              :: eps
 
 ! Output variables

src/GW/unrestricted_self_energy_correlation_diag.f90

@@ -20,7 +20,7 @@ subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSt,e,
 
 ! Local variables
 
-  integer                       :: i,j,a,b,p,q,jb
+  integer                       :: i,a,p,q,jb
   double precision              :: eps
 
 ! Output variables

src/MBPT/Makefile

@@ -1,10 +0,0 @@
-default:
-	ninja
-	make -C ..
-
-clean:
-	ninja -t clean
-
-debug:
-	ninja -t clean
-	make -C .. debug

src/MBPT/Sangalli_A_matrix_dynamic.f90.x

@@ -1,78 +0,0 @@
-subroutine Sangalli_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,A_dyn)
-
-! Compute the dynamic part of the Bethe-Salpeter equation matrices
-
-  implicit none
-  include 'parameters.h'
-
-! Input variables
-
-  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
-  double precision,intent(in)   :: eta
-  double precision,intent(in)   :: lambda
-  double precision,intent(in)   :: eGW(nBas)
-  double precision,intent(in)   :: OmRPA(nS)
-  double precision,intent(in)   :: OmBSE
-  double precision,intent(in)   :: rho(nBas,nBas,nS)
-  
-! Local variables
-
-  integer                       :: maxS
-  double precision              :: chi
-  double precision              :: eps
-  integer                       :: i,j,a,b,ia,jb,kc
-
-! Output variables
-
-  double precision,intent(out)  :: A_dyn(nS,nS)
-
-! Initialization
-
-  A_dyn(:,:) = 0d0
-
-! Number of poles taken into account 
-
-  maxS = nS
-
-! Build dynamic A matrix
-
-  do 
-
-
-  (ERI(i,k,j,c)*KroneckerDelta(a,b) + ERI(k,a,c,b)*KroneckerDelta(i,j))*(X)
-
-  ia = 0
-  do i=nC+1,nO
-    do a=nO+1,nBas-nR
-      ia = ia + 1
-      jb = 0
-      do j=nC+1,nO
-        do b=nO+1,nBas-nR
-          jb = jb + 1
- 
-          chi = 0d0
-          do kc=1,maxS
-
-            chi = chi + rho(i,j,kc)*rho(a,b,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2)
-
-          enddo
-
-          A_dyn(ia,jb) = A_dyn(ia,jb) - 4d0*lambda*chi
-
-          chi = 0d0
-          do kc=1,maxS
-            do ld=1,maxS
-
-            eps = OmBSE - (OmRPA(kc) + OmRPA(ld))
-            chi = chi + cRPA(ia,kc,ld)*cRPA(jb,kc,ld)*eps/(eps**2 + (2d0*eta)**2)
-
-          enddo
-
-          A_dyn(ia,jb) = A_dyn(ia,jb) - 2d0*lambda*chi
-
-        enddo
-      enddo
-    enddo
-  enddo
-
-end subroutine Sangalli_A_matrix_dynamic

src/MBPT/excitation_density_RI.f90

@@ -1,65 +0,0 @@
-subroutine excitation_density_RI(nBas,nC,nO,nR,nS,c,G,XpY,rho)
-
-! Compute excitation densities in the RI approximation
-
-  implicit none
-
-! Input variables
-
-  integer,intent(in)            :: nBas,nC,nO,nR,nS
-  double precision,intent(in)   :: c(nBas,nBas),G(nBas,nBas,nBas,nBas),XpY(nS,nS)
-
-! Local variables
-
-  double precision,allocatable  :: scr(:,:,:)
-  integer                       :: mu,nu,la,si,ia,jb,x,y,j,b
-
-! Output variables
-
-  double precision,intent(out)  :: rho(nBas,nBas,nS)
-
-! Memory allocation
-  allocate(scr(nBas,nBas,nS))
-
-  rho(:,:,:) = 0d0
-  do nu=1,nBas
-    do si=1,nBas
-      do ia=1,nS
-        jb = 0
-        do j=nC+1,nO
-          do b=nO+1,nBas-nR
-            jb = jb + 1
-            rho(nu,si,ia) = rho(nu,si,ia) + c(nu,j)*XpY(ia,jb)*c(si,b)
-          enddo
-        enddo
-      enddo
-    enddo
-  enddo
- 
-  scr(:,:,:) = 0d0   
-  do mu=1,nBas
-    do la=1,nBas
-      do ia=1,nS
-        do nu=1,nBas
-          do si=1,nBas
-            scr(mu,la,ia) = scr(mu,la,ia) + G(mu,nu,la,si)*rho(nu,si,ia)
-          enddo
-        enddo
-      enddo
-    enddo
-  enddo
-
-  rho(:,:,:) = 0d0   
-  do ia=1,nS
-    do x=nC+1,nBas-nR
-      do y=nC+1,nBas-nR
-        do mu=1,nBas
-          do la=1,nBas
-            rho(x,y,ia) = rho(x,y,ia) + c(mu,x)*scr(mu,la,ia)*c(la,y)
-          enddo
-        enddo
-      enddo
-    enddo
-  enddo
-
-end subroutine excitation_density_RI

src/MBPT/excitation_density_SOSEX_RI.f90

@@ -1,65 +0,0 @@
-subroutine excitation_density_SOSEX_RI(nBas,nC,nO,nR,nS,c,G,XpY,rho)
-
-! Compute excitation densities
-
-  implicit none
-
-! Input variables
-
-  integer,intent(in)            :: nBas,nC,nO,nR,nS
-  double precision,intent(in)   :: c(nBas,nBas),G(nBas,nBas,nBas,nBas),XpY(nS,nS)
-
-! Local variables
-
-  double precision,allocatable  :: scr(:,:,:)
-  integer                       :: mu,nu,la,si,ia,jb,x,y,j,b
-
-! Output variables
-
-  double precision,intent(out)  :: rho(nBas,nBas,nS)
-
-! Memory allocation
-  allocate(scr(nBas,nBas,nS))
-
-  rho(:,:,:) = 0d0
-  do nu=1,nBas
-    do si=1,nBas
-      do ia=1,nS
-        jb = 0
-        do j=nC+1,nO
-          do b=nO+1,nBas-nR
-            jb = jb + 1
-            rho(nu,si,ia) = rho(nu,si,ia) + c(nu,j)*XpY(ia,jb)*c(si,b)
-          enddo
-        enddo
-      enddo
-    enddo
-  enddo
- 
-  scr(:,:,:) = 0d0   
-  do mu=1,nBas
-    do la=1,nBas
-      do ia=1,nS
-        do nu=1,nBas
-          do si=1,nBas
-            scr(mu,la,ia) = scr(mu,la,ia) + G(mu,nu,la,si)*rho(nu,si,ia)
-          enddo
-        enddo
-      enddo
-    enddo
-  enddo
-
-  rho(:,:,:) = 0d0   
-  do ia=1,nS
-    do x=nC+1,nBas-nR
-      do y=nC+1,nBas-nR
-        do mu=1,nBas
-          do la=1,nBas
-            rho(x,y,ia) = rho(x,y,ia) + c(mu,x)*scr(mu,la,ia)*c(la,y)
-          enddo
-        enddo
-      enddo
-    enddo
-  enddo
-
-end subroutine excitation_density_SOSEX_RI

src/QuAcK/QuAcK.f90

@@ -124,14 +124,19 @@ program QuAcK
 
   integer                       :: maxSCF_GF,n_diis_GF,renormGF
   double precision              :: thresh_GF
-  logical                       :: DIIS_GF,linGF
+  logical                       :: DIIS_GF,linGF,regGF
   double precision              :: eta_GF
 
   integer                       :: maxSCF_GW,n_diis_GW
   double precision              :: thresh_GW
-  logical                       :: DIIS_GW,COHSEX,SOSEX,TDA_W,G0W,GW0,linGW
+  logical                       :: DIIS_GW,COHSEX,SOSEX,TDA_W,G0W,GW0,linGW,regGW
   double precision              :: eta_GW
 
+  integer                       :: maxSCF_GT,n_diis_GT
+  double precision              :: thresh_GT
+  logical                       :: DIIS_GT,TDA_T,linGT,regGT
+  double precision              :: eta_GT
+
   logical                       :: BSE,dBSE,dTDA,evDyn
 
   integer                       :: nMC,nEq,nWalk,nPrint,iSeed
@@ -176,9 +181,10 @@ program QuAcK
   call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,dostab, &
                     maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC,                                  &
                     TDA,singlet,triplet,spin_conserved,spin_flip,                           &
-                    maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,            &
+                    maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,regGF,      &
-                    maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,                     & 
+                    maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,regGW,               & 
                     COHSEX,SOSEX,TDA_W,G0W,GW0,                                             &  
+                    maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T,         & 
                     doACFDT,exchange_kernel,doXBS,                                          &
                     BSE,dBSE,dTDA,evDyn,                                                    &
                     nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift)
@@ -868,13 +874,14 @@ program QuAcK
 
     if(unrestricted) then
 
-      call UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linGF,eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
+      call UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linGF,eta_GF,regGF,        &
-                 S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb,eHF)
+                 nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & 
+                 dipole_int_aa,dipole_int_bb,eHF)
 
     else
 
-      call G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linGF, & 
+      call G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linGF,eta_GF,regGF, & 
-                eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
+                nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
 
     end if
 
@@ -897,13 +904,13 @@ program QuAcK
     if(unrestricted) then
 
       call evUGF2(maxSCF_GF,thresh_GF,n_diis_GF,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,    &
-                  eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & 
+                  eta_GF,regGF,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & 
                   dipole_int_aa,dipole_int_bb,cHF,eHF)
 
     else
 
       call evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF_GF,thresh_GF,n_diis_GF, & 
-                 singlet,triplet,linGF,eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & 
+                 singlet,triplet,linGF,eta_GF,regGF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & 
                  ERI_MO,dipole_int_MO,eHF)
 
     end if
@@ -926,13 +933,13 @@ program QuAcK
 
     if(unrestricted) then
 
-      call qsUGF2(maxSCF_GF,thresh_GF,n_diis_GF,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GF, &
+      call qsUGF2(maxSCF_GF,thresh_GF,n_diis_GF,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GF,regGF, &
                   nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO,                        & 
                   ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF)
 
     else
 
-      call qsGF2(maxSCF_GF,thresh_GF,n_diis_GF,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GF,nNuc,ZNuc,rNuc,ENuc, & 
+      call qsGF2(maxSCF_GF,thresh_GF,n_diis_GF,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GF,regGF,nNuc,ZNuc,rNuc,ENuc, & 
                  nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
 
     end if
@@ -989,7 +996,7 @@ program QuAcK
     if(unrestricted) then 
 
       call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,   & 
-                 linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & 
+                 linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & 
                  dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eG0W0)
     else
 
@@ -997,10 +1004,10 @@ program QuAcK
 
       if(SOSEX) then
         call G0W0_SOSEX(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, &
-                        eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0)
+                        eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0)
       else
         call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, &
-                  linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0)
+                  linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0)
       end if
 
     end if
@@ -1023,14 +1030,14 @@ program QuAcK
     if(unrestricted) then 
 
       call evUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,   &
-                G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,    &
+                G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,    &
                 EUHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb, & 
                 PHF,cHF,eHF,Vxc,eG0W0)
 
     else
 
       call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX, &
-                BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, &
+                BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW,regGW, &
                 nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0)
     end if
     call cpu_time(end_evGW)
@@ -1052,14 +1059,14 @@ program QuAcK
     if(unrestricted) then 
     
       call qsUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, &
-                 G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO, &
+                 G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,regGW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO, &
                  nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_AO,      &
                  dipole_int_aa,dipole_int_bb,PHF,cHF,eHF)
 
     else
  
       call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,               &
-                BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW,nNuc,ZNuc,rNuc,ENuc, & 
+                BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW,regGW,nNuc,ZNuc,rNuc,ENuc, & 
                 nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
 
     end if
@@ -1079,7 +1086,7 @@ program QuAcK
   if(doufG0W0) then
     
     call cpu_time(start_ufGW)
-    call ufG0W0(eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
+    call ufG0W0(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
     call cpu_time(end_ufGW)
   
     t_ufGW = end_ufGW - start_ufGW
@@ -1095,14 +1102,14 @@ program QuAcK
   if(doufGW) then
     
     call cpu_time(start_ufGW)
-    call ufGW(eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
+    call ufGW(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
     call cpu_time(end_ufGW)
   
     t_ufGW = end_ufGW - start_ufGW
     write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufGW = ',t_ufGW,' seconds'
     write(*,*)
 
-!   call ufBSE(eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0)
+!   call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0)
 
   end if
 
@@ -1115,8 +1122,8 @@ program QuAcK
   if(doG0T0) then
     
     call cpu_time(start_G0T0)
-    call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, &
+    call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,  &
-              linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,      &  
+              linGT,eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO, &  
               PHF,cHF,eHF,Vxc,eG0T0)
     call cpu_time(end_G0T0)
   
@@ -1133,8 +1140,8 @@ program QuAcK
   if(doevGT) then
     
     call cpu_time(start_evGT)
-    call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, &
+    call evGT(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, &
-              BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GW,        & 
+              BSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GT,regGT,  & 
               nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,   &
               PHF,cHF,eHF,Vxc,eG0T0)
     call cpu_time(end_evGT)
@@ -1153,9 +1160,10 @@ program QuAcK
 
     call cpu_time(start_qsGT)
 
-    call qsGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,         &
+    call qsGT(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, &
-              BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GW,nNuc,ZNuc,rNuc,ENuc, & 
+              BSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GT,regGT,  &
-              nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
+              nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,     & 
+              ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
 
     call cpu_time(end_qsGT)
 

src/QuAcK/read_options.f90

@@ -1,9 +1,10 @@
 subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,dostab, &
                         maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC,                                  &
                         TDA,singlet,triplet,spin_conserved,spin_flip,                           &
-                        maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,            &
+                        maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,regGF,      &
-                        maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,                     &
+                        maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,regGW,               &
                         COHSEX,SOSEX,TDA_W,G0W,GW0,                                             &
+                        maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T,         &
                         doACFDT,exchange_kernel,doXBS,                                          &
                         BSE,dBSE,dTDA,evDyn,                                                    &
                         nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift)
@@ -41,6 +42,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
   logical,intent(out)           :: linGF
   integer,intent(out)           :: renormGF
   double precision,intent(out)  :: eta_GF
+  logical,intent(out)           :: regGF
 
   integer,intent(out)           :: maxSCF_GW
   double precision,intent(out)  :: thresh_GW
@@ -53,6 +55,16 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
   logical,intent(out)           :: GW0
   logical,intent(out)           :: linGW
   double precision,intent(out)  :: eta_GW
+  logical,intent(out)           :: regGW
+
+  integer,intent(out)           :: maxSCF_GT
+  double precision,intent(out)  :: thresh_GT
+  logical,intent(out)           :: DIIS_GT
+  integer,intent(out)           :: n_diis_GT
+  logical,intent(out)           :: TDA_T
+  logical,intent(out)           :: linGT
+  double precision,intent(out)  :: eta_GT
+  logical,intent(out)           :: regGT
 
   logical,intent(out)           :: doACFDT
   logical,intent(out)           :: exchange_kernel
@@ -73,7 +85,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
   
 ! Local variables
 
-  character(len=1)              :: answer1,answer2,answer3,answer4,answer5,answer6,answer7
+  character(len=1)              :: answer1,answer2,answer3,answer4,answer5,answer6,answer7,answer8
 
 ! Open file with method specification
 
@@ -135,7 +147,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
   if(answer4 == 'T') spin_conserved = .true.
   if(answer5 == 'T') spin_flip      = .true.
 
-! Read Green function options
+! Read GF options
 
   maxSCF_GF = 64
   thresh_GF = 1d-5
@@ -144,12 +156,14 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
   linGF     = .false.
   eta_GF    = 0d0
   renormGF  = 0
+  regGF     = .false.
 
   read(1,*) 
-  read(1,*) maxSCF_GF,thresh_GF,answer1,n_diis_GF,answer2,eta_GF,renormGF
+  read(1,*) maxSCF_GF,thresh_GF,answer1,n_diis_GF,answer2,eta_GF,renormGF,answer3
 
   if(answer1 == 'T') DIIS_GF = .true.
   if(answer2 == 'T') linGF   = .true.
+  if(answer3 == 'T') regGF   = .true.
   if(.not.DIIS_GF) n_diis_GF = 1
 
 ! Read GW options
@@ -160,6 +174,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
   n_diis_GW = 5
   linGW     = .false.
   eta_GW    = 0d0
+  regGW     = .false.
   COHSEX    = .false.
   SOSEX     = .false.
   TDA_W     = .false.
@@ -168,7 +183,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
 
   read(1,*) 
   read(1,*) maxSCF_GW,thresh_GW,answer1,n_diis_GW,answer2,eta_GW, &
-            answer3,answer4,answer5,answer6,answer7
+            answer3,answer4,answer5,answer6,answer7,answer8
 
   if(answer1 == 'T') DIIS_GW = .true.
   if(answer2 == 'T') linGW   = .true.
@@ -177,8 +192,30 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
   if(answer5 == 'T') TDA_W   = .true.
   if(answer6 == 'T') G0W     = .true.
   if(answer7 == 'T') GW0     = .true.
+  if(answer8 == 'T') regGW   = .true.
   if(.not.DIIS_GW) n_diis_GW = 1
 
+! Read GF options
+
+  maxSCF_GF = 64
+  thresh_GF = 1d-5
+  DIIS_GF   = .false.
+  n_diis_GF = 5
+  linGF     = .false.
+  eta_GF    = 0d0
+  regGF     = .false.
+  TDA_T     = .false.
+
+  read(1,*) 
+  read(1,*) maxSCF_GT,thresh_GT,answer1,n_diis_GT,answer2,eta_GT, &
+            answer3,answer4
+
+  if(answer1 == 'T') DIIS_GT = .true.
+  if(answer2 == 'T') linGT   = .true.
+  if(answer3 == 'T') TDA_T   = .true.
+  if(answer4 == 'T') regGT   = .true.
+  if(.not.DIIS_GT) n_diis_GT = 1
+
 ! Options for adiabatic connection
 
   doACFDT         = .false.