PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-04-24 17:24:53 +02:00
Code Releases Activity

fixing GW regularization

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-08-01 16:30:41 +02:00
parent 0dfa3b0071
commit 01d8e28ca2
13 changed files with 104 additions and 688 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
src/GW/G0W0.f90
View File

@ -114,16 +114,9 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
! Compute GW self-energy !
!------------------------!
if(regularize) then
if(regularize) call GW_regularization(nBas,nC,nO,nR,nS,eHF,Om,rho)
call regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC)
call regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,Z)
else
call GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
end if
call GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
!-----------------------------------!
! Solve the quasi-particle equation !

52
src/GW/GW_regularization.f90 Normal file
View File

@ -0,0 +1,52 @@
subroutine GW_regularization(nBas,nC,nO,nR,nS,e,Om,rho)
! Regularize GW 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):: rho(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)
rho(p,i,m) = kappa*rho(p,i,m)
enddo
do a=nO+1,nBas-nR
Dpam = e(p) - e(a) - Om(m)
kappa = 1d0 - exp(-Dpam*Dpam*s)
rho(p,a,m) = kappa*rho(p,a,m)
enddo
enddo
enddo
end subroutine

35
src/GW/UG0W0.f90
View File

@ -53,10 +53,10 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons
double precision,allocatable :: SigC(:,:)
double precision,allocatable :: Z(:,:)
integer :: nS_aa,nS_bb,nS_sc
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:,:)
double precision,allocatable :: Om(:)
double precision,allocatable :: XpY(:,:)
double precision,allocatable :: XmY(:,:)
double precision,allocatable :: rho(:,:,:,:)
double precision,allocatable :: eGWlin(:,:)
@ -98,7 +98,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons
nS_sc = nS_aa + nS_bb
allocate(SigC(nBas,nspin),Z(nBas,nspin),eGWlin(nBas,nspin), &
OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin))
Om(nS_sc),XpY(nS_sc,nS_sc),XmY(nS_sc,nS_sc),rho(nBas,nBas,nS_sc,nspin))
!-------------------!
! Compute screening !
@ -109,31 +109,28 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons
ispin = 1
call phULR(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,Om,rho,EcRPA,Om,XpY,XmY)
if(print_W) call print_excitation_energies('phRPA@UHF',5,nS_sc,OmRPA)
if(print_W) call print_excitation_energies('phRPA@UHF',5,nS_sc,Om)
!----------------------!
! Excitation densities !
!----------------------!
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho)
!------------------------------------------------!
! Compute self-energy and renormalization factor !
!------------------------------------------------!
if(regularize) then
call unrestricted_regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,EcGM)
call unrestricted_regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z)
else
call UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,Z,EcGM)
do is=1,nspin
call GW_regularization(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eHF(:,is),Om,rho(:,:,:,is))
end do
end if
call UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,Om,rho,SigC,Z,EcGM)
!-----------------------------------!
! Solve the quasi-particle equation !
!-----------------------------------!
@ -153,7 +150,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons
do is=1,nspin
call UGW_QP_graph(eta,nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eHF(:,is), &
OmRPA,rho_RPA(:,:,:,is),eHF(:,is),eGW(:,is),Z(:,is))
Om,rho(:,:,:,is),eHF(:,is),eGW(:,is),Z(:,is))
end do
end if
@ -161,7 +158,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons
! Compute RPA correlation energy
call phULR(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Om,rho,EcRPA,Om,XpY,XmY)
! Dump results
@ -169,7 +166,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons
! Free memory
deallocate(OmRPA,XpY_RPA,XmY_RPA,rho_RPA)
deallocate(Om,XpY,XmY,rho)
! Perform BSE calculation

13
src/GW/evGW.f90
View File

@ -123,22 +123,15 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
! Compute spectral weights
! Compute spectral weights
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho)
! Compute correlation part of the self-energy
if(regularize) then
if(regularize) call GW_regularization(nBas,nC,nO,nR,nS,eGW,Om,rho)
call regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC)
call renormalization_factor_SRG(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,Z)
else
call GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
end if
call GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
! Solve the quasi-particle equation

33
src/GW/evUGW.f90
View File

@ -65,10 +65,10 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
double precision,allocatable :: Z(:,:)
integer :: nS_aa,nS_bb,nS_sc
double precision,allocatable :: SigC(:,:)
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:,:)
double precision,allocatable :: Om(:)
double precision,allocatable :: XpY(:,:)
double precision,allocatable :: XmY(:,:)
double precision,allocatable :: rho(:,:,:,:)
! Hello world
@ -104,7 +104,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
nS_sc = nS_aa + nS_bb
allocate(eGW(nBas,nspin),eOld(nBas,nspin),Z(nBas,nspin),SigC(nBas,nspin), &
OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin), &
Om(nS_sc),XpY(nS_sc,nS_sc),XmY(nS_sc,nS_sc),rho(nBas,nBas,nS_sc,nspin), &
error_diis(nBas,max_diis,nspin),e_diis(nBas,max_diis,nspin))
! Initialization
@ -129,29 +129,26 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
! Compute screening
call phULR(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Om,rho,EcRPA,Om,XpY,XmY)
!----------------------!
! Excitation densities !
!----------------------!
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho)
!------------------------------------------------!
! Compute self-energy and renormalization factor !
!------------------------------------------------!
if(regularize) then
call unrestricted_regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,EcGM)
call unrestricted_regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,Z)
else
call UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,Z,EcGM)
if(regularize) then
do is=1,nspin
call GW_regularization(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eGW(:,is),Om,rho(:,:,:,is))
end do
end if
call UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,Om,rho,SigC,Z,EcGM)
!-----------------------------------!
! Solve the quasi-particle equation !
!-----------------------------------!
@ -170,7 +167,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
do is=1,nspin
call UGW_QP_graph(eta,nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eHF(:,is), &
OmRPA,rho_RPA(:,:,:,is),eOld(:,is),eGW(:,is),Z(:,is))
Om,rho(:,:,:,is),eOld(:,is),eGW(:,is),Z(:,is))
end do
end if
@ -232,7 +229,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
! Deallocate memory
deallocate(eOld,Z,SigC,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error_diis,e_diis)
deallocate(eOld,Z,SigC,Om,XpY,XmY,rho,error_diis,e_diis)
! Perform BSE calculation

11
src/GW/qsGW.f90
View File

@ -182,16 +182,9 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,rho)
if(regularize) then
if(regularize) call GW_regularization(nBas,nC,nO,nR,nS,eGW,Om,rho)
call regularized_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC)
call regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,Z)
else
call GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
endif
call GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
! Make correlation self-energy Hermitian and transform it back to AO basis

31
src/GW/qsUGW.f90
View File

@ -81,10 +81,10 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
double precision,external :: trace_matrix
double precision,allocatable :: error_diis(:,:,:)
double precision,allocatable :: F_diis(:,:,:)
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:,:)
double precision,allocatable :: Om(:)
double precision,allocatable :: XpY(:,:)
double precision,allocatable :: XmY(:,:)
double precision,allocatable :: rho(:,:,:,:)
double precision,allocatable :: c(:,:,:)
double precision,allocatable :: cp(:,:,:)
double precision,allocatable :: eGW(:,:)
@ -139,8 +139,8 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
allocate(eGW(nBas,nspin),eOld(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin),F(nBas,nBas,nspin), &
Fp(nBas,nBas,nspin),J(nBas,nBas,nspin),K(nBas,nBas,nspin),SigC(nBas,nBas,nspin),SigCp(nBas,nBas,nspin), &
SigCm(nBas,nBas,nspin),Z(nBas,nspin),OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc), &
rho_RPA(nBas,nBas,nS_sc,nspin),error(nBas,nBas,nspin),error_diis(nBasSq,max_diis,nspin), &
SigCm(nBas,nBas,nspin),Z(nBas,nspin),Om(nS_sc),XpY(nS_sc,nS_sc),XmY(nS_sc,nS_sc), &
rho(nBas,nBas,nS_sc,nspin),error(nBas,nBas,nspin),error_diis(nBasSq,max_diis,nspin), &
F_diis(nBasSq,max_diis,nspin))
! Initialization
@ -198,29 +198,26 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
! Compute linear response
call phULR(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Om,rho,EcRPA,Om,XpY,XmY)
!----------------------!
! Excitation densities !
!----------------------!
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho)
!------------------------------------------------!
! Compute self-energy and renormalization factor !
!------------------------------------------------!
if(regularize) then
call unrestricted_regularized_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,EcGM)
call unrestricted_regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,Z)
else
call UGW_self_energy(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,Z,EcGM)
do is=1,nspin
call GW_regularization(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eGW(:,is),Om,rho(:,:,:,is))
end do
end if
call UGW_self_energy(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,Om,rho,SigC,Z,EcGM)
! Make correlation self-energy Hermitian and transform it back to AO basis
do is=1,nspin
@ -355,7 +352,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
! Deallocate memory
deallocate(cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error,error_diis,F_diis)
deallocate(cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,Om,XpY,XmY,rho,error,error_diis,F_diis)
! Perform BSE calculation

73
src/GW/regularized_renormalization_factor.f90
View File

@ -1,73 +0,0 @@
subroutine regularized_renormalization_factor(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
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 = 1d0
! 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(-2d0*eps**2/kappa**2))/eps
dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2)
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(-2d0*eps**2/kappa**2))/eps
dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2)
Z(p) = Z(p) - 2d0*rho(p,a,jb)**2*dfk
end do
end do
end do
! Compute renormalization factor from derivative of SigC
Z(:) = 1d0/(1d0 - Z(:))
end subroutine

91
src/GW/regularized_self_energy_correlation.f90
View File

@ -1,91 +0,0 @@
subroutine regularized_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC)
! Compute correlation part of the regularized self-energy
implicit none
include 'parameters.h'
! Input variables
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,j,a,b
integer :: p,q,r
integer :: jb
double precision :: eps
double precision :: kappa
double precision :: fk
! Output variables
double precision,intent(out) :: SigC(nBas,nBas)
double precision,intent(out) :: EcGM
! Initialize
SigC(:,:) = 0d0
!---------------------------------------------!
! Parameters for regularized MP2 calculations !
!---------------------------------------------!
kappa = 1d0
!----------------!
! GW self-energy !
!----------------!
! Occupied part of the correlation self-energy
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do i=nC+1,nO
do jb=1,nS
eps = e(p) - e(i) + Omega(jb)
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,jb)*rho(q,i,jb)*fk
end do
end do
end do
end do
! Virtual part of the correlation self-energy
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(p) - e(a) - Omega(jb)
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,jb)*rho(q,a,jb)*fk
end do
end do
end do
end do
! GM correlation energy
EcGM = 0d0
do i=nC+1,nO
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(a) - e(i) + Omega(jb)
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
EcGM = EcGM - 4d0*rho(a,i,jb)**2*fk
end do
end do
end do
end subroutine

72
src/GW/regularized_self_energy_correlation_diag.f90
View File

@ -1,72 +0,0 @@
subroutine regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC)
! Compute diagonal of the correlation part of the regularized self-energy
implicit none
include 'parameters.h'
! Input variables
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,q,jb
double precision :: Dpijb,Dpajb
! Output variables
double precision,intent(out) :: SigC(nBas)
double precision,intent(out) :: EcGM
! Initialize
SigC(:) = 0d0
!-----------------------------
! GW self-energy
!-----------------------------
! Occupied part of the correlation self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
do jb=1,nS
Dpijb = e(p) - e(i) + Omega(jb)
SigC(p) = SigC(p) + 2d0*rho(p,i,jb)**2*(1d0 - exp(-2d0*eta*Dpijb*Dpijb))/Dpijb
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
Dpajb = e(p) - e(a) - Omega(jb)
SigC(p) = SigC(p) + 2d0*rho(p,a,jb)**2*(1d0 - exp(-2d0*eta*Dpajb*Dpajb))/Dpajb
end do
end do
end do
! Galitskii-Migdal correlation energy
EcGM = 0d0
do i=nC+1,nO
do a=nO+1,nBas-nR
do jb=1,nS
EcGM = EcGM - 4d0*rho(a,i,jb)**2
end do
end do
end do
end subroutine

111
src/GW/unrestricted_regularized_renormalization_factor.f90
View File

@ -1,111 +0,0 @@
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

133
src/GW/unrestricted_regularized_self_energy_correlation.f90
View File

@ -1,133 +0,0 @@
subroutine unrestricted_regularized_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_regularized_self_energy_correlation

126
src/GW/unrestricted_regularized_self_energy_correlation_diag.f90
View File

@ -1,126 +0,0 @@
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