diff --git a/src/GW/RG0W0.f90 b/src/GW/RG0W0.f90 index aa24eae..451d31c 100644 --- a/src/GW/RG0W0.f90 +++ b/src/GW/RG0W0.f90 @@ -1,5 +1,5 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & - linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) + linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! Perform G0W0 calculation @@ -25,7 +25,7 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA logical,intent(in) :: triplet logical,intent(in) :: linearize double precision,intent(in) :: eta - logical,intent(in) :: regularize + logical,intent(in) :: doSRG integer,intent(in) :: nBas integer,intent(in) :: nOrb @@ -87,6 +87,15 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA write(*,*) end if +! SRG regularization + + if(doSRG) then + + write(*,*) '*** SRG regularized qsGW scheme ***' + write(*,*) + + end if + ! Memory allocation allocate(Aph(nS,nS),Bph(nS,nS),SigC(nOrb),Z(nOrb),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS), & @@ -113,9 +122,11 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! Compute GW self-energy ! !------------------------! - if(regularize) call GW_regularization(nOrb,nC,nO,nV,nR,nS,eHF,Om,rho) - - call RGW_self_energy_diag(eta,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z) + if(doSRG) then + call RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z) + else + call RGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z) + end if !-----------------------------------! ! Solve the quasi-particle equation ! diff --git a/src/GW/RGW_SRG_self_energy.f90 b/src/GW/RGW_SRG_self_energy.f90 new file mode 100644 index 0000000..4e13078 --- /dev/null +++ b/src/GW/RGW_SRG_self_energy.f90 @@ -0,0 +1,144 @@ +subroutine RGW_SRG_self_energy(nBas,nOrb,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) :: nOrb + 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(nOrb) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: rho(nOrb,nOrb,nS) + +! Local variables + + integer :: i,j,a,b + integer :: p,q + integer :: m + double precision :: Dpim,Dqim,Dpam,Dqam,Diam + double precision :: renorm + double precision :: s + +! Output variables + + double precision,intent(out) :: EcGM + double precision,intent(out) :: SigC(nOrb,nOrb) + double precision,intent(out) :: Z(nOrb) + +!--------------------! +! SRG flow parameter ! +!--------------------! + + s = 500d0 + +!--------------------! +! SRG-GW self-energy ! +!--------------------! + + SigC(:,:) = 0d0 + + ! Occupied part of the correlation self-energy + + !$OMP PARALLEL & + !$OMP SHARED(SigC,rho,s,nS,nC,nO,nOrb,nR,e,Om) & + !$OMP PRIVATE(m,i,q,p,Dpim,Dqim,renorm) & + !$OMP DEFAULT(NONE) + !$OMP DO + do q=nC+1,nOrb-nR + do p=nC+1,nOrb-nR + do m=1,nS + do i=nC+1,nO + + 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 + + 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,nOrb,e,Om) & + !$OMP PRIVATE(m,a,q,p,Dpam,Dqam,renorm) & + !$OMP DEFAULT(NONE) + !$OMP DO + do q=nC+1,nOrb-nR + do p=nC+1,nOrb-nR + do m=1,nS + do a=nO+1,nOrb-nR + + 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 + + end do + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + +!------------------------! +! Renormalization factor ! +!------------------------! + + Z(:) = 0d0 + + ! Occupied part of the renormlization factor + + do p=nC+1,nOrb-nR + do i=nC+1,nO + do m=1,nS + + Dpim = e(p) - e(i) + Om(m) + Z(p) = Z(p) - 2d0*rho(p,i,m)**2*(1d0-exp(-2d0*s*Dpim*Dpim))/Dpim**2 + + end do + end do + end do + + ! Virtual part of the renormlization factor + + do p=nC+1,nOrb-nR + do a=nO+1,nOrb-nR + do m=1,nS + Dpam = e(p) - e(a) - Om(m) + Z(p) = Z(p) - 2d0*rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam**2 + end do + end do + end do + + Z(:) = 1d0/(1d0 - Z(:)) + +!-------------------------------------! +! Galitskii-Migdal correlation energy ! +!-------------------------------------! + + EcGM = 0d0 + do i=nC+1,nO + do a=nO+1,nOrb-nR + do m=1,nS + + Diam = e(a) - e(i) + Om(m) + EcGM = EcGM - 4d0*rho(a,i,m)*rho(a,i,m)*(1d0-exp(-2d0*s*Diam*Diam))/Diam + + end do + end do + end do + +end subroutine diff --git a/src/GW/RGW_SRG_self_energy_diag.f90 b/src/GW/RGW_SRG_self_energy_diag.f90 new file mode 100644 index 0000000..f640e59 --- /dev/null +++ b/src/GW/RGW_SRG_self_energy_diag.f90 @@ -0,0 +1,138 @@ +subroutine RGW_SRG_self_energy_diag(nBas,nOrb,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) :: nOrb + 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(nOrb) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: rho(nOrb,nOrb,nS) + +! Local variables + + integer :: i,j,a,b + integer :: p + integer :: m + double precision :: Dpim,Dqim,Dpam,Dqam,Diam + double precision :: renorm + double precision :: s + +! Output variables + + double precision,intent(out) :: EcGM + double precision,intent(out) :: SigC(nOrb) + double precision,intent(out) :: Z(nOrb) + +!--------------------! +! SRG flow parameter ! +!--------------------! + + s = 500d0 + +!--------------------! +! SRG-GW self-energy ! +!--------------------! + + SigC(:) = 0d0 + + ! Occupied part of the correlation self-energy + + !$OMP PARALLEL & + !$OMP SHARED(SigC,rho,s,nS,nC,nO,nOrb,nR,e,Om) & + !$OMP PRIVATE(m,i,p,Dpim,Dqim,renorm) & + !$OMP DEFAULT(NONE) + !$OMP DO + do p=nC+1,nOrb-nR + do m=1,nS + 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 + + 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,nOrb,e,Om) & + !$OMP PRIVATE(m,a,p,Dpam,Dqam,renorm) & + !$OMP DEFAULT(NONE) + !$OMP DO + do p=nC+1,nOrb-nR + do m=1,nS + 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 + + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + +!------------------------! +! Renormalization factor ! +!------------------------! + + Z(:) = 0d0 + + ! Occupied part of the renormlization factor + + do p=nC+1,nOrb-nR + do i=nC+1,nO + do m=1,nS + + Dpim = e(p) - e(i) + Om(m) + Z(p) = Z(p) - 2d0*rho(p,i,m)**2*(1d0-exp(-2d0*s*Dpim*Dpim))/Dpim**2 + + end do + end do + end do + + ! Virtual part of the renormlization factor + + do p=nC+1,nOrb-nR + do a=nO+1,nOrb-nR + do m=1,nS + Dpam = e(p) - e(a) - Om(m) + Z(p) = Z(p) - 2d0*rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam**2 + end do + end do + end do + + Z(:) = 1d0/(1d0 - Z(:)) + +!-------------------------------------! +! Galitskii-Migdal correlation energy ! +!-------------------------------------! + + EcGM = 0d0 + do i=nC+1,nO + do a=nO+1,nOrb-nR + do m=1,nS + + Diam = e(a) - e(i) + Om(m) + EcGM = EcGM - 4d0*rho(a,i,m)*rho(a,i,m)*(1d0-exp(-2d0*s*Diam*Diam))/Diam + + end do + end do + end do + +end subroutine diff --git a/src/GW/RGW_self_energy.f90 b/src/GW/RGW_self_energy.f90 index 9e368be..0f0ce8a 100644 --- a/src/GW/RGW_self_energy.f90 +++ b/src/GW/RGW_self_energy.f90 @@ -1,4 +1,4 @@ -subroutine RGW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z) +subroutine RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z) ! Compute correlation part of the self-energy and the renormalization factor @@ -9,14 +9,15 @@ subroutine RGW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z) double precision,intent(in) :: eta integer,intent(in) :: nBas + integer,intent(in) :: nOrb 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) :: e(nOrb) double precision,intent(in) :: Om(nS) - double precision,intent(in) :: rho(nBas,nBas,nS) + double precision,intent(in) :: rho(nOrb,nOrb,nS) ! Local variables @@ -27,71 +28,118 @@ subroutine RGW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z) ! Output variables double precision,intent(out) :: EcGM - double precision,intent(out) :: Sig(nBas,nBas) - double precision,intent(out) :: Z(nBas) - -! Initialize - - Sig(:,:) = 0d0 - Z(:) = 0d0 + double precision,intent(out) :: Sig(nOrb,nOrb) + double precision,intent(out) :: Z(nOrb) !----------------! ! GW self-energy ! !----------------! + Sig(:,:) = 0d0 + ! Occupied part of the correlation self-energy -!$OMP PARALLEL & -!$OMP SHARED(Sig,Z,rho,eta,nS,nC,nO,nBas,nR,e,Om) & -!$OMP PRIVATE(m,i,q,p,eps,num) & -!$OMP DEFAULT(NONE) -!$OMP DO - do q=nC+1,nBas-nR - do p=nC+1,nBas-nR - do m=1,nS - do i=nC+1,nO + !$OMP PARALLEL & + !$OMP SHARED(Sig,Z,rho,eta,nS,nC,nO,nOrb,nR,e,Om) & + !$OMP PRIVATE(m,i,q,p,eps,num) & + !$OMP DEFAULT(NONE) + !$OMP DO + do q=nC+1,nOrb-nR + do p=nC+1,nOrb-nR + do m=1,nS + do i=nC+1,nO - eps = e(p) - e(i) + Om(m) - num = 2d0*rho(p,i,m)*rho(q,i,m) - Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) - if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + eps = e(p) - e(i) + Om(m) + num = 2d0*rho(p,i,m)*rho(q,i,m) + Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) - end do end do - 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(Sig,Z,rho,eta,nS,nC,nO,nOrb,nR,e,Om) & + !$OMP PRIVATE(m,a,q,p,eps,num) & + !$OMP DEFAULT(NONE) + !$OMP DO + do q=nC+1,nOrb-nR + do p=nC+1,nOrb-nR + do m=1,nS + do a=nO+1,nOrb-nR + + eps = e(p) - e(a) - Om(m) + num = 2d0*rho(p,a,m)*rho(q,a,m) + Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) + + end do + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + +!------------------------! +! Renormalization factor ! +!------------------------! + + Z(:) = 0d0 + +! Occupied part of the renormalization factor + !$OMP PARALLEL & -!$OMP SHARED(Sig,Z,rho,eta,nS,nC,nO,nBas,nR,e,Om) & -!$OMP PRIVATE(m,a,q,p,eps,num) & +!$OMP SHARED(Sig,Z,rho,eta,nS,nC,nO,nOrb,nR,e,Om) & +!$OMP PRIVATE(m,i,q,p,eps,num) & !$OMP DEFAULT(NONE) -!$OMP DO - do q=nC+1,nBas-nR - do p=nC+1,nBas-nR - do m=1,nS - do a=nO+1,nBas-nR +!$OMP DO + do p=nC+1,nOrb-nR + do m=1,nS + do i=nC+1,nO - eps = e(p) - e(a) - Om(m) - num = 2d0*rho(p,a,m)*rho(q,a,m) - Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) - if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + eps = e(p) - e(i) + Om(m) + num = 2d0*rho(p,i,m)*rho(q,i,m) + Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - end do end do end do end do -!$OMP END DO -!$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL -! Galitskii-Migdal correlation energy +! Virtual part of the renormalization factor + + !$OMP PARALLEL & + !$OMP SHARED(Sig,Z,rho,eta,nS,nC,nO,nOrb,nR,e,Om) & + !$OMP PRIVATE(m,a,q,p,eps,num) & + !$OMP DEFAULT(NONE) + !$OMP DO + do p=nC+1,nOrb-nR + do m=1,nS + do a=nO+1,nOrb-nR + + eps = e(p) - e(a) - Om(m) + num = 2d0*rho(p,a,m)*rho(q,a,m) + Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + + Z(:) = 1d0/(1d0 - Z(:)) + +!-------------------------------------! +! Galitskii-Migdal correlation energy ! +!-------------------------------------! EcGM = 0d0 do m=1,nS - do a=nO+1,nBas-nR + do a=nO+1,nOrb-nR do i=nC+1,nO eps = e(a) - e(i) + Om(m) @@ -102,8 +150,4 @@ subroutine RGW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z) end do end do -! Compute renormalization factor from derivative - - Z(:) = 1d0/(1d0 - Z(:)) - end subroutine diff --git a/src/GW/SRG_qsRGW.f90 b/src/GW/SRG_qsRGW.f90 index 42ee7b7..c664b31 100644 --- a/src/GW/SRG_qsRGW.f90 +++ b/src/GW/SRG_qsRGW.f90 @@ -272,7 +272,6 @@ subroutine SRG_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS ! Save quasiparticles energy for next cycle - Conv = maxval(abs(err)) eOld(:) = eGW(:) !------------------------------------------------------------------------ diff --git a/src/GW/evRGW.f90 b/src/GW/evRGW.f90 index 118ca64..7361fa0 100644 --- a/src/GW/evRGW.f90 +++ b/src/GW/evRGW.f90 @@ -1,5 +1,5 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & - singlet,triplet,linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) + singlet,triplet,linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! Perform self-consistent eigenvalue-only GW calculation @@ -29,7 +29,7 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop logical,intent(in) :: triplet logical,intent(in) :: linearize double precision,intent(in) :: eta - logical,intent(in) :: regularize + logical,intent(in) :: doSRG integer,intent(in) :: nBas integer,intent(in) :: nOrb @@ -125,9 +125,13 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! Compute correlation part of the self-energy - if(regularize) call GW_regularization(nOrb,nC,nO,nV,nR,nS,eGW,Om,rho) + call GW_regularization(nOrb,nC,nO,nV,nR,nS,eGW,Om,rho) - call RGW_self_energy_diag(eta,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) + if(doSRG) then + call RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) + else + call RGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) + end if ! Solve the quasi-particle equation diff --git a/src/GW/qsRGW.f90 b/src/GW/qsRGW.f90 index 61fc434..fe3b9b8 100644 --- a/src/GW/qsRGW.f90 +++ b/src/GW/qsRGW.f90 @@ -1,5 +1,5 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, & - TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc, & + TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet,eta,doSRG,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) @@ -28,7 +28,7 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop logical,intent(in) :: singlet logical,intent(in) :: triplet double precision,intent(in) :: eta - logical,intent(in) :: regularize + logical,intent(in) :: doSRG integer,intent(in) :: nNuc double precision,intent(in) :: ZNuc(nNuc) @@ -123,6 +123,15 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop write(*,*) end if +! SRG regularization + + if(doSRG) then + + write(*,*) '*** SRG regularized qsGW scheme ***' + write(*,*) + + end if + ! Memory allocation allocate(eGW(nOrb)) @@ -195,12 +204,13 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) if(print_W) call print_excitation_energies('phRPA@GW@RHF','singlet',nS,Om) - call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI_MO,XpY,rho) - if(regularize) call GW_regularization(nOrb,nC,nO,nV,nR,nS,eGW,Om,rho) - - call RGW_self_energy(eta,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) + if(doSRG) then + call RGW_SRG_self_energy(nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) + else + call RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) + end if ! Make correlation self-energy Hermitian and transform it back to AO basis