From 179a0835b30305ae0cc991d6ae25c0f9a4519386 Mon Sep 17 00:00:00 2001 From: Loris Burth Date: Tue, 22 Apr 2025 19:58:53 +0200 Subject: [PATCH] Added SRG for complex GW methods --- src/GW/RGW_SRG_self_energy.f90 | 2 +- src/GW/RGW_self_energy.f90 | 3 +- src/GW/complex_RGW_QP_graph.f90 | 9 +- src/GW/complex_RGW_SRG_SigC_dSigC.f90 | 88 +++++++++ src/GW/complex_RGW_SRG_self_energy.f90 | 193 ++++++++++++++++++++ src/GW/complex_RGW_SRG_self_energy_diag.f90 | 8 +- src/GW/complex_RGW_self_energy.f90 | 2 +- src/GW/complex_cRG0W0.f90 | 7 +- src/GW/complex_evRGW.f90 | 11 +- src/GW/complex_qsRGW.f90 | 6 +- 10 files changed, 310 insertions(+), 19 deletions(-) create mode 100644 src/GW/complex_RGW_SRG_SigC_dSigC.f90 create mode 100644 src/GW/complex_RGW_SRG_self_energy.f90 diff --git a/src/GW/RGW_SRG_self_energy.f90 b/src/GW/RGW_SRG_self_energy.f90 index 6ad1c90..53b4521 100644 --- a/src/GW/RGW_SRG_self_energy.f90 +++ b/src/GW/RGW_SRG_self_energy.f90 @@ -71,7 +71,7 @@ subroutine RGW_SRG_self_energy(flow,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC, !$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) & diff --git a/src/GW/RGW_self_energy.f90 b/src/GW/RGW_self_energy.f90 index b3d9c2c..4888c23 100644 --- a/src/GW/RGW_self_energy.f90 +++ b/src/GW/RGW_self_energy.f90 @@ -130,7 +130,7 @@ subroutine RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z) end do !$OMP END DO !$OMP END PARALLEL - + Z(:) = 1d0/(1d0 - Z(:)) !-------------------------------------! @@ -149,5 +149,4 @@ subroutine RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z) end do end do end do - end subroutine diff --git a/src/GW/complex_RGW_QP_graph.f90 b/src/GW/complex_RGW_QP_graph.f90 index a4af64d..46144b9 100644 --- a/src/GW/complex_RGW_QP_graph.f90 +++ b/src/GW/complex_RGW_QP_graph.f90 @@ -64,11 +64,16 @@ subroutine complex_RGW_QP_graph(doSRG,eta,flow,nBas,nC,nO,nV,nR,nS,Re_eHF,Im_eHF nIt = nIt + 1 + if(doSRG) then - call complex_RGW_SigC_dSigC(p,eta,nBas,nC,nO,nV,nR,nS,& + call complex_RGW_SRG_SigC_dSigC(flow,p,eta,nBas,nC,nO,nV,nR,nS,& Re_w,Im_w,Re_eOld,Im_eOld,Om,rho,& Re_SigC,Im_SigC,Re_dSigC,Im_dSigC) - + else + call complex_RGW_SigC_dSigC(p,eta,nBas,nC,nO,nV,nR,nS,& + Re_w,Im_w,Re_eOld,Im_eOld,Om,rho,& + Re_SigC,Im_SigC,Re_dSigC,Im_dSigC) + end if Re_f = Re_w - Re_eHF(p) - Re_SigC Im_f = Im_w - Im_eHF(p) - Im_SigC Re_df = (1d0 - Re_dSigC)/((1d0 - Re_dSigC)**2 + Im_dSigC**2) diff --git a/src/GW/complex_RGW_SRG_SigC_dSigC.f90 b/src/GW/complex_RGW_SRG_SigC_dSigC.f90 new file mode 100644 index 0000000..8c2fa38 --- /dev/null +++ b/src/GW/complex_RGW_SRG_SigC_dSigC.f90 @@ -0,0 +1,88 @@ +subroutine complex_RGW_SRG_SigC_dSigC(flow,p,eta,nBas,nC,nO,nV,nR,nS,Re_w,Im_w,Re_e,Im_e,Om,rho,Re_SigC,Im_SigC,Re_DS,Im_DS) + + ! Complute diagonal of the correlation part of the self-energy and its derivative fully complex + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: p + double precision,intent(in) :: eta + double precision,intent(in) :: flow + 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) :: Re_e(nBas) + double precision,intent(in) :: Im_e(nBas) + double precision,intent(in) :: Re_w + double precision,intent(in) :: Im_w + complex*16,intent(in) :: Om(nS) + complex*16,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + integer :: i,a,m + double precision :: eps,s + double precision :: eta_tilde + complex*16 :: num + complex*16 :: tmp + +! Output variables + + double precision,intent(out) :: Re_SigC + double precision,intent(out) :: Im_SigC + double precision,intent(out) :: Re_DS + double precision,intent(out) :: Im_DS + +! Initialize + Re_SigC = 0d0 + Im_SigC = 0d0 + Re_DS = 0d0 + Im_DS = 0d0 + s = flow + +! Compute self energy and its derivative + +! Occupied part + + do i=nC+1,nO + do m=1,nS + eps = Re_w - Re_e(i) + real(Om(m)) + eta_tilde = eta - Im_w + Im_e(i) - aimag(Om(m)) + num = 2d0*rho(p,i,m)**2*(1d0 - exp(-2d0*s*(eps**2 + eta_tilde**2))) + + tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),& + eta_tilde/(eps**2 + eta_tilde**2),kind=8) + Re_SigC = Re_SigC + real(tmp) + Im_SigC = Im_SigC + aimag(tmp) + + tmp = num*cmplx(-(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2,& + -2*eta_tilde*eps/(eps**2 + eta_tilde**2)**2,kind=8) + Re_DS = Re_DS + real(tmp) + Im_DS = Im_DS + aimag(tmp) + end do + end do + +! Virtual part + + do a=nO+1,nBas-nR + do m=1,nS + eps = Re_w - Re_e(a) - real(Om(m)) + eta_tilde = eta + Im_w - Im_e(a) - aimag(Om(m)) + num = 2d0*rho(p,a,m)**2*(1d0 - exp(-2d0*s*(eps**2 + eta_tilde**2))) + + tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),-eta_tilde/(eps**2 + eta_tilde**2),kind=8) + Re_SigC = Re_SigC + real(tmp) + Im_SigC = Im_SigC + aimag(tmp) + + tmp = num*cmplx(-(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2,& + 2*eta_tilde*eps/(eps**2 + eta_tilde**2)**2,kind=8) + Re_DS = Re_DS + real(tmp) + Im_DS = Im_DS + aimag(tmp) + end do + end do +end subroutine diff --git a/src/GW/complex_RGW_SRG_self_energy.f90 b/src/GW/complex_RGW_SRG_self_energy.f90 new file mode 100644 index 0000000..ae3a03e --- /dev/null +++ b/src/GW/complex_RGW_SRG_self_energy.f90 @@ -0,0 +1,193 @@ +subroutine complex_RGW_SRG_self_energy(flow,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 + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + double precision,intent(in) :: flow + 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 + complex*16,intent(in) :: e(nOrb) + complex*16,intent(in) :: Om(nS) + complex*16,intent(in) :: rho(nOrb,nOrb,nS) + +! Local variables + + integer :: i,j,a,b + integer :: p,q,m + double precision :: eps_p,eps_q,eta_tilde_p,eta_tilde_q + double precision :: eps,eta_tilde,s + complex*16 :: num,tmp + double precision,allocatable :: Re_DS(:) + double precision,allocatable :: Im_DS(:) + double precision,allocatable :: Re_Sig(:,:) + double precision,allocatable :: Im_Sig(:,:) + double precision,allocatable :: Re_Z(:) + double precision,allocatable :: Im_Z(:) + +! Output variables + + complex*16,intent(out) :: EcGM + complex*16,intent(out) :: Sig(nOrb,nOrb) + complex*16,intent(out) :: Z(nOrb) + +!----------------! +! GW self-energy ! +!----------------! + allocate(Re_DS(nOrb),Im_DS(nOrb),Re_Z(nOrb),Im_Z(nOrb),Re_Sig(nOrb,nOrb),Im_Sig(nOrb,nOrb)) + + Re_Sig(:,:) = 0d0 + Im_Sig(:,:) = 0d0 + Re_DS(:) = 0d0 + Im_DS(:) = 0d0 + s = flow + +! Occupied part of the correlation self-energy + + !$OMP PARALLEL & + !$OMP SHARED(Re_Sig,Im_Sig,rho,eta,nS,nC,nO,nOrb,nR,e,Om,s) & + !$OMP PRIVATE(m,i,q,p,eps_p,eps_q,num,eta_tilde_p,eta_tilde_q,tmp) & + !$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_p = real(e(p)) - real(e(i)) + real(Om(m)) + eps_q = real(e(q)) - real(e(i)) + real(Om(m)) + eta_tilde_p = eta - aimag(e(p)) + aimag(e(i)) - aimag(Om(m)) + eta_tilde_q = eta - aimag(e(q)) + aimag(e(i)) - aimag(Om(m)) + num = 2d0*rho(p,i,m)*rho(q,i,m)*(1d0 - exp(-s*(eps_p**2+eta_tilde_p**2 + eps_q**2 + eta_tilde_q**2))) + tmp = num*cmplx((eps_p + eps_q)/(eps_p**2 + eps_q**2 + eta_tilde_p**2 + eta_tilde_q**2),& + (eta_tilde_p + eta_tilde_q)/(eps_p**2 + eps_q**2 + eta_tilde_p**2 + eta_tilde_q**2),kind=8) + Re_Sig(p,q) = Re_Sig(p,q) + real(tmp) + Im_Sig(p,q) = Im_Sig(p,q) + aimag(tmp) + + 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(Re_Sig,Im_Sig,rho,eta,nS,nC,nO,nOrb,nR,e,Om,s) & + !$OMP PRIVATE(m,a,q,p,eps_p,eps_q,num,eta_tilde_p,eta_tilde_q,tmp) & + !$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_p = real(e(p)) - real(e(a)) - real(Om(m)) + eps_q = real(e(q)) - real(e(a)) - real(Om(m)) + eta_tilde_p = eta + aimag(e(p)) - aimag(e(a)) - aimag(Om(m)) + eta_tilde_q = eta + aimag(e(q)) - aimag(e(a)) - aimag(Om(m)) + num = 2d0*rho(p,a,m)*rho(q,a,m)*(1d0 - exp(-s*(eps_p**2+eta_tilde_p**2 + eps_q**2 + eta_tilde_q**2))) + tmp = num*cmplx((eps_p +eps_q)/(eps_p**2+eta_tilde_p**2 + eps_q**2 + eta_tilde_q**2),& + -(eta_tilde_p + eta_tilde_q)/(eps_p**2+eta_tilde_p**2 + eps_q**2 + eta_tilde_q**2),kind=8) + Re_Sig(p,q) = Re_Sig(p,q) + real(tmp) + Im_Sig(p,q) = Im_Sig(p,q) + aimag(tmp) + + end do + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + +!------------------------! +! Renormalization factor ! +!------------------------! + + !Occupied part of the renormalization factor + + !$OMP PARALLEL & + !$OMP SHARED(Re_DS,Im_DS,rho,eta,nS,nC,nO,nOrb,nR,e,Om,s) & + !$OMP PRIVATE(m,i,p,eps,num,eta_tilde,tmp) & + !$OMP DEFAULT(NONE) + !$OMP DO + do p=nC+1,nOrb-nR + do m=1,nS + do i=nC+1,nO + eps = real(e(p)) - real(e(i)) + real(Om(m)) + eta_tilde = eta - aimag(e(p)) + aimag(e(i)) - aimag(Om(m)) + num = 2d0*rho(p,i,m)*rho(p,i,m)*(1d0-exp(-2d0*s*(eps**2+eta_tilde**2))) + tmp = num*cmplx(-(eps**2-eta_tilde**2)/(eps**2 + eta_tilde**2)**2,& + -2*eta_tilde*eps/(eps**2 + eta_tilde**2)**2,kind=8) + Re_DS(p) = Re_DS(p) + real(tmp) + Im_DS(p) = Im_DS(p) + aimag(tmp) + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + +! Virtual part of the renormalization factor + + !$OMP PARALLEL & + !$OMP SHARED(Re_DS,Im_DS,rho,eta,nS,nC,nO,nOrb,nR,e,Om,s) & + !$OMP PRIVATE(m,a,p,eps,num,eta_tilde,tmp) & + !$OMP DEFAULT(NONE) + !$OMP DO + do p=nC+1,nOrb-nR + do m=1,nS + do a=nO+1,nOrb-nR + + eps = real(e(p)) - real(e(a)) - real(Om(m)) + eta_tilde = eta + aimag(e(p)) - aimag(e(a)) - aimag(Om(m)) + num = 2d0*rho(p,a,m)*rho(p,a,m)*(1d0-exp(-2d0*s*(eps**2+eta_tilde**2))) + tmp = num*cmplx(-(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2,& + 2*eta_tilde*eps/eps/(eps**2 + eta_tilde**2)**2,kind=8) + Re_DS(p) = Re_DS(p) + real(tmp) + Im_DS(p) = Im_DS(p) + aimag(tmp) + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + +! Compute renormalization factor from derivative + Re_Z(:) = (1d0-Re_DS(:))/((1d0 - Re_DS(:))**2 + Im_DS(:)**2) + Im_Z(:) = Im_DS(:)/((1d0 - Re_DS(:))**2 + Im_DS(:)**2) + + Z = cmplx(Re_Z,Im_Z,kind=8) + Sig = cmplx(Re_Sig,Im_Sig,kind=8) + + deallocate(Re_DS) + deallocate(Im_DS) + deallocate(Re_Z) + deallocate(Im_Z) + deallocate(Re_Sig) + deallocate(Im_Sig) + +!!-------------------------------------! +!! Galitskii-Migdal correlation energy ! +!!-------------------------------------! +! +! EcGM = 0d0 +! do m=1,nS +! do a=nO+1,nOrb-nR +! do i=nC+1,nO +! +! eps = e(a) - e(i) + Om(m) +! num = 4d0*rho(a,i,m)*rho(a,i,m) +! EcGM = EcGM - num*eps/(eps**2 + eta**2) +! +! end do +! end do +! end do +! +end subroutine diff --git a/src/GW/complex_RGW_SRG_self_energy_diag.f90 b/src/GW/complex_RGW_SRG_self_energy_diag.f90 index c9e9b5f..223466f 100644 --- a/src/GW/complex_RGW_SRG_self_energy_diag.f90 +++ b/src/GW/complex_RGW_SRG_self_energy_diag.f90 @@ -24,7 +24,7 @@ subroutine complex_RGW_SRG_self_energy_diag(flow,eta,nBas,nOrb,nC,nO,nV,nR,nS,Re ! Local variables integer :: i,a,p,m - double precision :: eps + double precision :: eps,s complex*16 :: num double precision :: eta_tilde double precision,allocatable :: Re_DS(:) @@ -46,6 +46,8 @@ subroutine complex_RGW_SRG_self_energy_diag(flow,eta,nBas,nOrb,nC,nO,nV,nR,nS,Re Re_DS(:) = 0d0 Im_DS(:) = 0d0 + s = flow + !----------------! ! GW self-energy ! !----------------! @@ -56,7 +58,7 @@ subroutine complex_RGW_SRG_self_energy_diag(flow,eta,nBas,nOrb,nC,nO,nV,nR,nS,Re do m=1,nS eps = Re_e(p) - Re_e(i) + real(Om(m)) eta_tilde = eta - Im_e(p) + Im_e(i) - aimag(Om(m)) - num = 2d0*rho(p,i,m)**2 + num = 2d0*rho(p,i,m)**2*(1d0 - exp(-2d0*s*(eps**2 + eta_tilde**2))) tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),& eta_tilde/(eps**2+eta_tilde**2),kind=8) Re_Sig(p) = Re_Sig(p) + real(tmp) @@ -78,7 +80,7 @@ subroutine complex_RGW_SRG_self_energy_diag(flow,eta,nBas,nOrb,nC,nO,nV,nR,nS,Re eps = Re_e(p) - Re_e(a) - real(Om(m)) eta_tilde = eta + Im_e(p) - Im_e(a) - aimag(Om(m)) - num = 2d0*rho(p,a,m)**2 + num = 2d0*rho(p,a,m)**2*(1d0 - exp(-2d0*s*(eps**2 + eta_tilde**2))) tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),& -eta_tilde/(eps**2 + eta_tilde**2),kind=8) Re_Sig(p) = Re_Sig(p) + real(tmp) diff --git a/src/GW/complex_RGW_self_energy.f90 b/src/GW/complex_RGW_self_energy.f90 index 8dc0447..6fcae20 100644 --- a/src/GW/complex_RGW_self_energy.f90 +++ b/src/GW/complex_RGW_self_energy.f90 @@ -41,7 +41,7 @@ subroutine complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Si !----------------! ! GW self-energy ! !----------------! - allocate(Re_DS(nBas),Im_DS(nBas),Re_Z(nOrb),Im_Z(nOrb),Re_Sig(nOrb,nOrb),Im_Sig(nOrb,nOrb)) + allocate(Re_DS(nOrb),Im_DS(nOrb),Re_Z(nOrb),Im_Z(nOrb),Re_Sig(nOrb,nOrb),Im_Sig(nOrb,nOrb)) Re_Sig(:,:) = 0d0 Im_Sig(:,:) = 0d0 diff --git a/src/GW/complex_cRG0W0.f90 b/src/GW/complex_cRG0W0.f90 index f7f129d..df5e28a 100644 --- a/src/GW/complex_cRG0W0.f90 +++ b/src/GW/complex_cRG0W0.f90 @@ -94,9 +94,7 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, flow = 500d0 if(doSRG) then - ! Not implemented write(*,*) '*** SRG regularized G0W0 scheme ***' - write(*,*) '!!! No SRG with cRG0W0 !!!' write(*,*) end if @@ -126,8 +124,11 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, !------------------------! ! Compute GW self-energy ! !------------------------! +if(doSRG) then + call complex_RGW_SRG_self_energy_diag(flow,eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_eHF,Im_eHF,Om,rho,EcGM,Re_SigC,Im_SigC,Re_Z,Im_Z) +else call complex_RGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_eHF,Im_eHF,Om,rho,EcGM,Re_SigC,Im_SigC,Re_Z,Im_Z) - +end if !-----------------------------------! ! Solve the quasi-particle equation ! !-----------------------------------! diff --git a/src/GW/complex_evRGW.f90 b/src/GW/complex_evRGW.f90 index f7ca17e..9865516 100644 --- a/src/GW/complex_evRGW.f90 +++ b/src/GW/complex_evRGW.f90 @@ -146,11 +146,14 @@ subroutine complex_evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,d call complex_RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho) - ! Compute correlation part of the self-energy - ! Implement here the srg version if(doSRG) .. complex_RGW_SRG_self_energy_diag - call complex_RGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_eGW,Im_eGW,Om,rho,& + if(doSRG) then + call complex_RGW_SRG_self_energy_diag(flow,eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_eGW,Im_eGW,Om,rho,& EcGM,Re_SigC,Im_SigC,Re_Z,Im_Z) - + else + call complex_RGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_eGW,Im_eGW,Om,rho,& + EcGM,Re_SigC,Im_SigC,Re_Z,Im_Z) + end if + ! Solve the quasi-particle equation if(linearize) then diff --git a/src/GW/complex_qsRGW.f90 b/src/GW/complex_qsRGW.f90 index c4a6b28..ef68f7b 100644 --- a/src/GW/complex_qsRGW.f90 +++ b/src/GW/complex_qsRGW.f90 @@ -213,10 +213,10 @@ subroutine complex_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,d if(print_W) call print_excitation_energies('phRPA@GW@RHF','singlet',nS,Om) call complex_RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI_MO,XpY,rho) - + if(doSRG) then - write(*,*) "SRG not implemented" - !call complex_RGW_SRG_self_energy(flow,nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) + call complex_RGW_SRG_self_energy(flow,eta,nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,& + rho,EcGM,SigC,Z) else call complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,& EcGM,SigC,Z)