GT regularization

src/GT/G0T0pp.f90

@@ -153,7 +153,10 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
 ! Compute T-matrix version of the self-energy 
 !----------------------------------------------
 
-! if(regularize) call GTpp_regularization(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t)
+  if(regularize) then 
+    call GTpp_regularization(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Om1s,rho1s,Om2s,rho2s)
+    call GTpp_regularization(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Om1t,rho1t,Om2t,rho2t)
+  end if
 
   call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, & 
                              Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z)

src/GT/GTeh_regularization.f90

@@ -0,0 +1,55 @@
+subroutine GTeh_regularization(nBas,nC,nO,nR,nS,e,Om,rhoL,rhoR)
+
+! Regularize GTeh excitation densities via SRG
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  integer,intent(in)            :: e(nBas)
+  integer,intent(in)            :: Om(nS)
+
+! Local variables
+
+  integer                       :: p,i,a,m
+  double precision              :: s
+  double precision              :: kappa
+  double precision              :: Dpim,Dpam
+
+! Output variables
+
+  double precision,intent(inout):: rhoL(nBas,nBas,nS)
+  double precision,intent(inout):: rhoR(nBas,nBas,nS)
+
+! SRG flow parameter 
+
+  s = 100d0
+
+! Regularize excitation densities
+
+  do p=nC+1,nBas-nR
+    do m=1,nS
+
+      do i=nC+1,nO
+        Dpim = e(p) - e(i) + Om(m)
+        kappa  = 1d0 - exp(-Dpim*Dpim*s)
+        rhoL(i,p,m) = kappa*rhoL(i,p,m)
+        rhoR(i,p,m) = kappa*rhoR(i,p,m)
+      enddo
+
+      do a=nO+1,nBas-nR
+        Dpam = e(p) - e(a) - Om(m)
+        kappa  = 1d0 - exp(-Dpam*Dpam*s)
+        rhoL(p,a,m) = kappa*rhoL(p,a,m)
+        rhoR(p,a,m) = kappa*rhoR(p,a,m)
+      enddo
+
+    enddo
+  enddo
+
+end subroutine 

src/GT/GTpp_regularization.f90

@@ -0,0 +1,58 @@
+subroutine GTpp_regularization(nBas,nC,nO,nR,nOO,nVV,e,Om1,rho1,Om2,rho2)
+
+! Regularize GTpp excitation densities via SRG
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nOO
+  integer,intent(in)            :: nVV
+  integer,intent(in)            :: e(nBas)
+
+  double precision,intent(in)   :: Om1(nVV)
+  double precision,intent(in)   :: Om2(nOO)
+
+! Local variables
+
+  integer                       :: p,i,a,m
+  double precision              :: s
+  double precision              :: kappa
+  double precision              :: Dpim,Dpam
+
+! Output variables
+
+  double precision,intent(inout):: rho1(nBas,nBas,nVV)
+  double precision,intent(inout):: rho2(nBas,nBas,nOO)
+
+! SRG flow parameter 
+
+  s = 100d0
+
+! Regularize excitation densities
+
+  do p=nC+1,nBas-nR
+
+    do m=1,nVV
+      do i=nC+1,nO
+        Dpim = e(p) + e(i) - Om1(m)
+        kappa  = 1d0 - exp(-Dpim*Dpim*s)
+        rho1(p,i,m) = kappa*rho1(p,i,m)
+      end do
+    end do
+
+    do m=1,nOO
+      do a=nO+1,nBas-nR
+        Dpam = e(p) + e(a) - Om2(m)
+        kappa  = 1d0 - exp(-Dpam*Dpam*s)
+        rho2(p,a,m) = kappa*rho2(p,a,m)
+      end do
+    end do
+
+  enddo
+
+end subroutine 

src/GT/evGTpp.f90

@@ -171,7 +171,10 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
   ! Compute T-matrix version of the self-energy 
   !----------------------------------------------
 
-  ! if(regularize) call GTpp_regularization(eta,nBas,nC,nO,nV,nR,nOO,nVV,eGT,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t)
+    if(regularize) then 
+      call GTpp_regularization(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s)
+      call GTpp_regularization(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t)
+    end if 
 
     call GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nooS,nVVt,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s, & 
                                Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z)

src/GT/qsGTpp.f90

@@ -223,7 +223,10 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T
     iblock =  4
     call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
 
-  ! if(regularize) call GTpp_regularization(eta,nBas,nC,nO,nV,nR,nOO,nVV,eGT,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t)
+    if(regularize) then 
+      call GTpp_regularization(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s)
+      call GTpp_regularization(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t)
+    end if 
 
     call GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s, & 
                           Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z)