diff --git a/src/GW/RGW_SRG_self_energy.f90 b/src/GW/RGW_SRG_self_energy.f90 index 4e13078..4fcc44f 100644 --- a/src/GW/RGW_SRG_self_energy.f90 +++ b/src/GW/RGW_SRG_self_energy.f90 @@ -24,7 +24,6 @@ subroutine RGW_SRG_self_energy(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) integer :: p,q integer :: m double precision :: Dpim,Dqim,Dpam,Dqam,Diam - double precision :: renorm double precision :: s ! Output variables @@ -49,7 +48,7 @@ subroutine RGW_SRG_self_energy(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) !$OMP PARALLEL & !$OMP SHARED(SigC,rho,s,nS,nC,nO,nOrb,nR,e,Om) & - !$OMP PRIVATE(m,i,q,p,Dpim,Dqim,renorm) & + !$OMP PRIVATE(m,i,q,p,Dpim,Dqim) & !$OMP DEFAULT(NONE) !$OMP DO do q=nC+1,nOrb-nR @@ -59,8 +58,9 @@ subroutine RGW_SRG_self_energy(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) Dpim = e(p) - e(i) + Om(m) Dqim = e(q) - e(i) + Om(m) - renorm = (1d0-exp(-s*Dpim*Dpim)*exp(-s*Dqim*Dqim))*(Dpim + Dqim)/(Dpim*Dpim + Dqim*Dqim) - SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,m)*rho(q,i,m)*renorm + SigC(p,q) = SigC(p,q) & + + 2d0*rho(p,i,m)*rho(q,i,m)* & + (1d0-exp(-s*Dpim*Dpim)*exp(-s*Dqim*Dqim))*(Dpim + Dqim)/(Dpim*Dpim + Dqim*Dqim) end do end do @@ -73,7 +73,7 @@ subroutine RGW_SRG_self_energy(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) !$OMP PARALLEL & !$OMP SHARED(SigC,rho,s,nS,nC,nO,nR,nOrb,e,Om) & - !$OMP PRIVATE(m,a,q,p,Dpam,Dqam,renorm) & + !$OMP PRIVATE(m,a,q,p,Dpam,Dqam) & !$OMP DEFAULT(NONE) !$OMP DO do q=nC+1,nOrb-nR @@ -83,8 +83,9 @@ subroutine RGW_SRG_self_energy(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) Dpam = e(p) - e(a) - Om(m) Dqam = e(q) - e(a) - Om(m) - renorm = (1d0-exp(-s*Dpam*Dpam)*exp(-s*Dqam*Dqam))*(Dpam + Dqam)/(Dpam*Dpam + Dqam*Dqam) - SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,m)*rho(q,a,m)*renorm + SigC(p,q) = SigC(p,q) & + + 2d0*rho(p,a,m)*rho(q,a,m)* & + (1d0-exp(-s*Dpam*Dpam)*exp(-s*Dqam*Dqam))*(Dpam + Dqam)/(Dpam*Dpam + Dqam*Dqam) end do end do diff --git a/src/GW/RGW_SRG_self_energy_diag.f90 b/src/GW/RGW_SRG_self_energy_diag.f90 index f640e59..41820e2 100644 --- a/src/GW/RGW_SRG_self_energy_diag.f90 +++ b/src/GW/RGW_SRG_self_energy_diag.f90 @@ -23,8 +23,7 @@ subroutine RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC, integer :: i,j,a,b integer :: p integer :: m - double precision :: Dpim,Dqim,Dpam,Dqam,Diam - double precision :: renorm + double precision :: Dpim,Dpam,Diam double precision :: s ! Output variables @@ -49,7 +48,7 @@ subroutine RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC, !$OMP PARALLEL & !$OMP SHARED(SigC,rho,s,nS,nC,nO,nOrb,nR,e,Om) & - !$OMP PRIVATE(m,i,p,Dpim,Dqim,renorm) & + !$OMP PRIVATE(m,i,p,Dpim) & !$OMP DEFAULT(NONE) !$OMP DO do p=nC+1,nOrb-nR @@ -57,8 +56,7 @@ subroutine RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC, do i=nC+1,nO Dpim = e(p) - e(i) + Om(m) - renorm = (1d0-exp(-2d0*s*Dpim*Dpim))/Dpim - SigC(p) = SigC(p) + 2d0*rho(p,i,m)**2*renorm + SigC(p) = SigC(p) + 2d0*rho(p,i,m)**2*(1d0-exp(-2d0*s*Dpim*Dpim))/Dpim end do end do @@ -70,7 +68,7 @@ subroutine RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC, !$OMP PARALLEL & !$OMP SHARED(SigC,rho,s,nS,nC,nO,nR,nOrb,e,Om) & - !$OMP PRIVATE(m,a,p,Dpam,Dqam,renorm) & + !$OMP PRIVATE(m,a,p,Dpam) & !$OMP DEFAULT(NONE) !$OMP DO do p=nC+1,nOrb-nR @@ -78,8 +76,7 @@ subroutine RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC, do a=nO+1,nOrb-nR Dpam = e(p) - e(a) - Om(m) - renorm = (1d0-exp(-2d0*s*Dpam*Dpam))/Dpam - SigC(p) = SigC(p) + 2d0*rho(p,a,m)**2*renorm + SigC(p) = SigC(p) + 2d0*rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam end do end do diff --git a/src/GW/SRG_qsUGW.f90 b/src/GW/SRG_qsUGW.f90 index 825eebd..04404cc 100644 --- a/src/GW/SRG_qsUGW.f90 +++ b/src/GW/SRG_qsUGW.f90 @@ -217,13 +217,7 @@ subroutine SRG_qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS ! Compute self-energy and renormalization factor ! !------------------------------------------------! - if(regularize) then - do is=1,nspin - call GW_regularization(nBas,nC(is),nO(is),nV(is),nR(is),nSt,eGW(:,is),Om,rho(:,:,:,is)) - end do - end if - - call USRG_self_energy(nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,EcGM,SigC,Z) + call UGW_SRG_self_energy(nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,EcGM,SigC,Z) ! Make correlation self-energy Hermitian and transform it back to AO basis diff --git a/src/GW/UG0W0.f90 b/src/GW/UG0W0.f90 index a32d867..7f8c8f9 100644 --- a/src/GW/UG0W0.f90 +++ b/src/GW/UG0W0.f90 @@ -1,5 +1,5 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, & - linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb, & + linearize,eta,doSRG,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb, & dipole_int_aa,dipole_int_bb,cHF,eHF) ! Perform unrestricted G0W0 calculation @@ -24,7 +24,7 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD logical,intent(in) :: spin_flip logical,intent(in) :: linearize double precision,intent(in) :: eta - logical,intent(in) :: regularize + logical,intent(in) :: doSRG integer,intent(in) :: nBas integer,intent(in) :: nC(nspin) @@ -127,14 +127,12 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD ! Compute self-energy and renormalization factor ! !------------------------------------------------! - if(regularize) then - do is=1,nspin - call GW_regularization(nBas,nC(is),nO(is),nV(is),nR(is),nSt,eHF(:,is),Om,rho(:,:,:,is)) - end do + if(doSRG) then + call UGW_SRG_self_energy_diag(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 - call UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nSt,eHF,Om,rho,SigC,Z,EcGM) - !-----------------------------------! ! Solve the quasi-particle equation ! !-----------------------------------! diff --git a/src/GW/USRG_self_energy.f90 b/src/GW/UGW_SRG_self_energy.f90 similarity index 51% rename from src/GW/USRG_self_energy.f90 rename to src/GW/UGW_SRG_self_energy.f90 index 18820b0..fc636d4 100644 --- a/src/GW/USRG_self_energy.f90 +++ b/src/GW/UGW_SRG_self_energy.f90 @@ -1,4 +1,4 @@ -subroutine USRG_self_energy(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) +subroutine UGW_SRG_self_energy(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! Compute correlation part of the self-energy @@ -21,10 +21,9 @@ subroutine USRG_self_energy(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) integer :: ispin integer :: i,j,a,b - integer :: p,q,r + integer :: p,q integer :: m double precision :: Dpim,Dqim,Dpam,Dqam,Diam - double precision :: t1,t2 double precision :: s ! Output variables @@ -47,97 +46,95 @@ subroutine USRG_self_energy(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! Occupied part of the correlation self-energy -! call wall_time(t1) - !$OMP PARALLEL & !$OMP SHARED(SigC,rho,s,nS,nC,nO,nBas,nR,e,Om) & !$OMP PRIVATE(ispin,m,i,q,p,Dpim,Dqim) & !$OMP DEFAULT(NONE) !$OMP DO do ispin=1,nspin - do q=nC(ispin)+1,nBas-nR(ispin) - do p=nC(ispin)+1,nBas-nR(ispin) + do q=nC(ispin)+1,nBas-nR(ispin) + do p=nC(ispin)+1,nBas-nR(ispin) do m=1,nS - do i=nC(ispin)+1,nO(ispin) - Dpim = e(p,ispin) - e(i,ispin) + Om(m) - Dqim = e(q,ispin) - e(i,ispin) + Om(m) - SigC(p,q,ispin) = SigC(p,q,ispin) & - + rho(p,i,m,ispin)*rho(q,i,m,ispin)*(1d0-dexp(-s*Dpim*Dpim)*dexp(-s*Dqim*Dqim)) & - *(Dpim + Dqim)/(Dpim*Dpim + Dqim*Dqim) - end do + do i=nC(ispin)+1,nO(ispin) + + Dpim = e(p,ispin) - e(i,ispin) + Om(m) + Dqim = e(q,ispin) - e(i,ispin) + Om(m) + SigC(p,q,ispin) = SigC(p,q,ispin) & + + rho(p,i,m,ispin)*rho(q,i,m,ispin)*(1d0-dexp(-s*Dpim*Dpim)*dexp(-s*Dqim*Dqim)) & + *(Dpim + Dqim)/(Dpim*Dpim + Dqim*Dqim) + + end do end do - end do - end do + end do + end do end do !$OMP END DO !$OMP END PARALLEL -! call wall_time(t2) -! print *, "first loop", (t2-t1) + ! Virtual part of the correlation self-energy -! Virtual part of the correlation self-energy - - call wall_time(t1) - !$OMP PARALLEL & - !$OMP SHARED(SigC,rho,s,nS,nC,nO,nR,nBas,e,Om) & - !$OMP PRIVATE(ispin,m,a,q,p,Dpam,Dqam) & - !$OMP DEFAULT(NONE) - !$OMP DO - do ispin=1,nspin - do q=nC(ispin)+1,nBas-nR(ispin) - do p=nC(ispin)+1,nBas-nR(ispin) - do m=1,nS + !$OMP PARALLEL & + !$OMP SHARED(SigC,rho,s,nS,nC,nO,nR,nBas,e,Om) & + !$OMP PRIVATE(ispin,m,a,q,p,Dpam,Dqam) & + !$OMP DEFAULT(NONE) + !$OMP DO + do ispin=1,nspin + do q=nC(ispin)+1,nBas-nR(ispin) + do p=nC(ispin)+1,nBas-nR(ispin) + do m=1,nS do a=nO(ispin)+1,nBas-nR(ispin) - Dpam = e(p,ispin) - e(a,ispin) - Om(m) - Dqam = e(q,ispin) - e(a,ispin) - Om(m) - SigC(p,q,ispin) = SigC(p,q,ispin) & - + rho(p,a,m,ispin)*rho(q,a,m,ispin)*(1d0-exp(-s*Dpam*Dpam)*exp(-s*Dqam*Dqam)) & - *(Dpam + Dqam)/(Dpam*Dpam + Dqam*Dqam) + + Dpam = e(p,ispin) - e(a,ispin) - Om(m) + Dqam = e(q,ispin) - e(a,ispin) - Om(m) + SigC(p,q,ispin) = SigC(p,q,ispin) & + + rho(p,a,m,ispin)*rho(q,a,m,ispin)*(1d0-exp(-s*Dpam*Dpam)*exp(-s*Dqam*Dqam)) & + *(Dpam + Dqam)/(Dpam*Dpam + Dqam*Dqam) + end do - end do + end do + end do end do - end do - end do + end do !$OMP END DO !$OMP END PARALLEL -! call wall_time(t2) -! print *, "second loop", (t2-t1) - - -! Initialize +!------------------------! +! Renormalization factor ! +!------------------------! Z(:,:) = 0d0 - do ispin=1,nspin - do p=nC(ispin)+1,nBas-nR(ispin) - do i=nC(ispin)+1,nO(ispin) - do m=1,nS - Dpim = e(p,ispin) - e(i,ispin) + Om(m) - Z(p,ispin) = Z(p,ispin) - rho(p,i,m,ispin)**2*(1d0-dexp(-2d0*s*Dpim*Dpim))/Dpim**2 - end do - end do - end do - end do - - ! Virtual part of the correlation self-energy + ! Occupied part of the renormalization factor do ispin=1,nspin - do p=nC(ispin)+1,nBas-nR(ispin) - do a=nO(ispin)+1,nBas-nR(ispin) + do p=nC(ispin)+1,nBas-nR(ispin) + do i=nC(ispin)+1,nO(ispin) do m=1,nS - Dpam = e(p,ispin) - e(a,ispin) - Om(m) - Z(p,ispin) = Z(p,ispin) - rho(p,a,m,ispin)**2*(1d0-dexp(-2d0*s*Dpam*Dpam))/Dpam**2 + Dpim = e(p,ispin) - e(i,ispin) + Om(m) + Z(p,ispin) = Z(p,ispin) - rho(p,i,m,ispin)**2*(1d0-dexp(-2d0*s*Dpim*Dpim))/Dpim**2 end do - end do - end do + end do + end do end do -! Compute renormalization factor from derivative of SigC + ! Virtual part of the renormalization factor + + do ispin=1,nspin + do p=nC(ispin)+1,nBas-nR(ispin) + do a=nO(ispin)+1,nBas-nR(ispin) + do m=1,nS + Dpam = e(p,ispin) - e(a,ispin) - Om(m) + Z(p,ispin) = Z(p,ispin) - rho(p,a,m,ispin)**2*(1d0-dexp(-2d0*s*Dpam*Dpam))/Dpam**2 + end do + end do + end do + end do Z(:,:) = 1d0/(1d0 - Z(:,:)) -! Galitskii-Migdal correlation energy +!-------------------------------------! +! Galitskii-Migdal correlation energy ! +!-------------------------------------! EcGM = 0d0 do ispin=1,nspin diff --git a/src/GW/UGW_SRG_self_energy_diag.f90 b/src/GW/UGW_SRG_self_energy_diag.f90 new file mode 100644 index 0000000..ee4b725 --- /dev/null +++ b/src/GW/UGW_SRG_self_energy_diag.f90 @@ -0,0 +1,141 @@ +subroutine UGW_SRG_self_energy_diag(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) + +! Compute correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + 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) :: nS + double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: rho(nBas,nBas,nS,nspin) + +! Local variables + + integer :: ispin + integer :: i,j,a,b + integer :: p + integer :: m + double precision :: Dpim,Dpam,Diam + double precision :: s + +! Output variables + + double precision,intent(out) :: EcGM + double precision,intent(out) :: SigC(nBas,nspin) + double precision,intent(out) :: Z(nBas,nspin) + +! SRG flow parameter + + s = 500d0 + +! Initialize + + SigC(:,:) = 0d0 + +!--------------------! +! SRG-GW self-energy ! +!--------------------! + + ! Occupied part of the correlation self-energy + + !$OMP PARALLEL & + !$OMP SHARED(SigC,rho,s,nS,nC,nO,nBas,nR,e,Om) & + !$OMP PRIVATE(ispin,m,i,p,Dpim) & + !$OMP DEFAULT(NONE) + !$OMP DO + do ispin=1,nspin + do p=nC(ispin)+1,nBas-nR(ispin) + do m=1,nS + do i=nC(ispin)+1,nO(ispin) + + Dpim = e(p,ispin) - e(i,ispin) + Om(m) + SigC(p,ispin) = SigC(p,ispin) + rho(p,i,m,ispin)**2*(1d0-dexp(-2d0*s*Dpim*Dpim))/Dpim + + end do + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + + ! Virtual part of the correlation self-energy + + !$OMP PARALLEL & + !$OMP SHARED(SigC,rho,s,nS,nC,nO,nR,nBas,e,Om) & + !$OMP PRIVATE(ispin,m,a,p,Dpam) & + !$OMP DEFAULT(NONE) + !$OMP DO + do ispin=1,nspin + do p=nC(ispin)+1,nBas-nR(ispin) + do m=1,nS + do a=nO(ispin)+1,nBas-nR(ispin) + + Dpam = e(p,ispin) - e(a,ispin) - Om(m) + SigC(p,ispin) = SigC(p,ispin) + rho(p,a,m,ispin)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam + + end do + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + +!------------------------! +! Renormalization factor ! +!------------------------! + + Z(:,:) = 0d0 + + ! Occupied part of the renormalization factor + + do ispin=1,nspin + do p=nC(ispin)+1,nBas-nR(ispin) + do i=nC(ispin)+1,nO(ispin) + do m=1,nS + Dpim = e(p,ispin) - e(i,ispin) + Om(m) + Z(p,ispin) = Z(p,ispin) - rho(p,i,m,ispin)**2*(1d0-dexp(-2d0*s*Dpim*Dpim))/Dpim**2 + end do + end do + end do + end do + + ! Virtual part of the renormalization factor + + do ispin=1,nspin + do p=nC(ispin)+1,nBas-nR(ispin) + do a=nO(ispin)+1,nBas-nR(ispin) + do m=1,nS + Dpam = e(p,ispin) - e(a,ispin) - Om(m) + Z(p,ispin) = Z(p,ispin) - rho(p,a,m,ispin)**2*(1d0-dexp(-2d0*s*Dpam*Dpam))/Dpam**2 + end do + end do + end do + end do + + Z(:,:) = 1d0/(1d0 - Z(:,:)) + +!-------------------------------------! +! Galitskii-Migdal correlation energy ! +!-------------------------------------! + + EcGM = 0d0 + do ispin=1,nspin + do i=nC(ispin)+1,nO(ispin) + do a=nO(ispin)+1,nBas-nR(ispin) + do m=1,nS + Diam = e(a,ispin) - e(i,ispin) + Om(m) + EcGM = EcGM - rho(a,i,m,ispin)*rho(a,i,m,ispin)*(1d0-exp(-2d0*s*Diam*Diam))/Diam + end do + end do + end do + end do + +end subroutine diff --git a/src/GW/evRGW.f90 b/src/GW/evRGW.f90 index 7361fa0..5229e96 100644 --- a/src/GW/evRGW.f90 +++ b/src/GW/evRGW.f90 @@ -125,8 +125,6 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! Compute correlation part of the self-energy - call GW_regularization(nOrb,nC,nO,nV,nR,nS,eGW,Om,rho) - if(doSRG) then call RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) else diff --git a/src/GW/evUGW.f90 b/src/GW/evUGW.f90 index 63e323c..f073194 100644 --- a/src/GW/evUGW.f90 +++ b/src/GW/evUGW.f90 @@ -1,5 +1,5 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA, & - spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc, & + spin_conserved,spin_flip,linearize,eta,doSRG,nBas,nC,nO,nV,nR,nS,ENuc, & EUHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF) ! Perform self-consistent eigenvalue-only GW calculation @@ -28,7 +28,7 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE logical,intent(in) :: spin_flip logical,intent(in) :: linearize double precision,intent(in) :: eta - logical,intent(in) :: regularize + logical,intent(in) :: doSRG integer,intent(in) :: nBas integer,intent(in) :: nC(nspin) @@ -153,13 +153,12 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE ! Compute self-energy and renormalization factor ! !------------------------------------------------! - if(regularize) then - do is=1,nspin - call GW_regularization(nBas,nC(is),nO(is),nV(is),nR(is),nSt,eGW(:,is),Om,rho(:,:,:,is)) - end do + if(doSRG) then + call UGW_SRG_self_energy_diag(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 - call UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,SigC,Z,EcGM) !-----------------------------------! ! Solve the quasi-particle equation ! diff --git a/src/GW/qsUGW.f90 b/src/GW/qsUGW.f90 index 69e3f11..0894911 100644 --- a/src/GW/qsUGW.f90 +++ b/src/GW/qsUGW.f90 @@ -1,5 +1,5 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,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, & + 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) ! Perform a quasiparticle self-consistent GW calculation @@ -25,7 +25,7 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE logical,intent(in) :: spin_conserved logical,intent(in) :: spin_flip double precision,intent(in) :: eta - logical,intent(in) :: regularize + logical,intent(in) :: doSRG integer,intent(in) :: nNuc double precision,intent(in) :: ZNuc(nNuc) @@ -217,13 +217,12 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE ! Compute self-energy and renormalization factor ! !------------------------------------------------! - if(regularize) then - do is=1,nspin - call GW_regularization(nBas,nC(is),nO(is),nV(is),nR(is),nSt,eGW(:,is),Om,rho(:,:,:,is)) - end do + if(doSRG) then + call UGW_SRG_self_energy(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 - call UGW_self_energy(eta,nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,SigC,Z,EcGM) ! Make correlation self-energy Hermitian and transform it back to AO basis