From 04a42701e2001a7b1d832244653a405fddf1f42d Mon Sep 17 00:00:00 2001 From: pfloos Date: Wed, 11 Sep 2024 18:41:27 +0200 Subject: [PATCH] working on SRG for GW --- src/GW/GG0W0.f90 | 37 ++-------------- src/GW/GGW_SRG_self_energy.f90 | 5 ++- src/GW/GGW_SRG_self_energy_diag.f90 | 5 ++- src/GW/RG0W0.f90 | 18 ++++---- src/GW/RGW_SRG_ReSigC.f90 | 52 ++++++++++++++++++++++ src/GW/RGW_SRG_self_energy.f90 | 5 ++- src/GW/RGW_SRG_self_energy_diag.f90 | 5 ++- src/GW/RGW_phACFDT.f90 | 6 +-- src/GW/RGW_phBSE.f90 | 1 - src/GW/UG0W0.f90 | 30 +++---------- src/GW/UGW_SRG_self_energy.f90 | 5 ++- src/GW/UGW_SRG_self_energy_diag.f90 | 5 ++- src/GW/UGW_phACFDT.f90 | 17 ++++++++ src/GW/UGW_phBSE.f90 | 17 +++++++- src/GW/evGGW.f90 | 67 +++-------------------------- src/GW/evRGW.f90 | 25 +++-------- src/GW/evUGW.f90 | 46 +++----------------- src/GW/qsGGW.f90 | 49 ++------------------- src/GW/qsRGW.f90 | 6 ++- src/GW/qsUGW.f90 | 38 ++++------------ 20 files changed, 160 insertions(+), 279 deletions(-) create mode 100644 src/GW/RGW_SRG_ReSigC.f90 diff --git a/src/GW/GG0W0.f90 b/src/GW/GG0W0.f90 index 038df13..109b2b0 100644 --- a/src/GW/GG0W0.f90 +++ b/src/GW/GG0W0.f90 @@ -40,6 +40,7 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA logical :: print_W = .true. logical :: dRPA + double precision :: flow double precision :: EcRPA double precision :: EcBSE double precision :: EcGM @@ -68,7 +69,6 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! Initialization dRPA = .true. - EcRPA = 0d0 ! TDA for W @@ -79,6 +79,8 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! SRG regularization + flow = 500d0 + if(doSRG) then write(*,*) '*** SRG regularized qsGW scheme ***' @@ -113,7 +115,7 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA !------------------------! if(doSRG) then - call GGW_SRG_self_energy_diag(nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z) + call GGW_SRG_self_energy_diag(flow,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z) else call GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z) end if @@ -142,8 +144,6 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA end if -! call GW_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eHF,Om,rho) - ! Compute the RPA correlation energy call phGLR_A(dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph) @@ -174,35 +174,6 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA write(*,*)'-------------------------------------------------------------------------------' write(*,*) -! Compute the BSE correlation energy via the adiabatic connection - -! if(doACFDT) then - -! write(*,*) '-------------------------------------------------------------' -! write(*,*) ' Adiabatic connection version of BSE@G0W0 correlation energy ' -! write(*,*) '-------------------------------------------------------------' -! write(*,*) - -! if(doXBS) then - -! write(*,*) '*** scaled screening version (XBS) ***' -! write(*,*) - -! end if - -! call GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,EcBSE) - -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy (singlet) =',EcBSE(1),' au' -! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy (triplet) =',EcBSE(2),' au' -! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2),' au' -! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 total energy =',ENuc + EGHF + EcBSE(1) + EcBSE(2),' au' -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) - -! end if - end if if(doppBSE) then diff --git a/src/GW/GGW_SRG_self_energy.f90 b/src/GW/GGW_SRG_self_energy.f90 index a0c7c70..a77fbf8 100644 --- a/src/GW/GGW_SRG_self_energy.f90 +++ b/src/GW/GGW_SRG_self_energy.f90 @@ -1,4 +1,4 @@ -subroutine GGW_SRG_self_energy(nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) +subroutine GGW_SRG_self_energy(flow,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! Compute correlation part of the self-energy @@ -7,6 +7,7 @@ subroutine GGW_SRG_self_energy(nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! Input variables + double precision,intent(in) :: flow integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO @@ -35,7 +36,7 @@ subroutine GGW_SRG_self_energy(nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! SRG flow parameter ! !--------------------! - s = 500d0 + s = flow !--------------------! ! SRG-GW self-energy ! diff --git a/src/GW/GGW_SRG_self_energy_diag.f90 b/src/GW/GGW_SRG_self_energy_diag.f90 index c63b2b3..e07199e 100644 --- a/src/GW/GGW_SRG_self_energy_diag.f90 +++ b/src/GW/GGW_SRG_self_energy_diag.f90 @@ -1,4 +1,4 @@ -subroutine GGW_SRG_self_energy_diag(nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) +subroutine GGW_SRG_self_energy_diag(flow,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! Compute correlation part of the self-energy @@ -7,6 +7,7 @@ subroutine GGW_SRG_self_energy_diag(nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! Input variables + double precision,intent(in) :: flow integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO @@ -35,7 +36,7 @@ subroutine GGW_SRG_self_energy_diag(nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! SRG flow parameter ! !--------------------! - s = 500d0 + s = flow !--------------------! ! SRG-GW self-energy ! diff --git a/src/GW/RG0W0.f90 b/src/GW/RG0W0.f90 index 451d31c..32c2a61 100644 --- a/src/GW/RG0W0.f90 +++ b/src/GW/RG0W0.f90 @@ -46,7 +46,7 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA logical :: plot_self = .false. logical :: dRPA_W integer :: isp_W - double precision :: lambda + double precision :: flow double precision :: EcRPA double precision :: EcBSE(nspin) double precision :: EcGM @@ -73,10 +73,6 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA write(*,*)'*******************************' write(*,*) -! Initialization - - lambda = 1d0 - ! Spin manifold and TDA for dynamical screening isp_W = 1 @@ -89,6 +85,8 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! SRG regularization + flow = 500d0 + if(doSRG) then write(*,*) '*** SRG regularized qsGW scheme ***' @@ -105,8 +103,8 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! Compute screening ! !-------------------! - call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,lambda,eHF,ERI,Aph) - if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,lambda,ERI,Bph) + call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) + if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) @@ -123,7 +121,7 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA !------------------------! if(doSRG) then - call RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z) + call RGW_SRG_self_energy_diag(flow,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 @@ -162,8 +160,8 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! Compute the RPA correlation energy - call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,lambda,eGW,ERI,Aph) - if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,lambda,ERI,Bph) + call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph) + if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) diff --git a/src/GW/RGW_SRG_ReSigC.f90 b/src/GW/RGW_SRG_ReSigC.f90 new file mode 100644 index 0000000..e15f607 --- /dev/null +++ b/src/GW/RGW_SRG_ReSigC.f90 @@ -0,0 +1,52 @@ +double precision function RGW_SRG_ReSigC(p,w,s,nBas,nC,nO,nV,nR,nS,e,Om,rho) + +! Compute diagonal of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: p + double precision,intent(in) :: w + double precision,intent(in) :: s + 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) :: Om(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + integer :: i,a,m + double precision :: Dpim,Dpam + +! Initialize + + RGW_SRG_ReSigC = 0d0 + +! Occupied part of the correlation self-energy + + do i=nC+1,nO + do m=1,nS + Dpim = w - e(i) + Om(m) + RGW_SRG_ReSigC = RGW_SRG_ReSigC & + + 2d0*rho(p,i,m)**2*(1d0-exp(-2d0*s*Dpim*Dpim))/Dpim + end do + end do + +! Virtual part of the correlation self-energy + + do a=nO+1,nBas-nR + do m=1,nS + Dpam = w - e(a) - Om(m) + RGW_SRG_ReSigC = RGW_SRG_ReSigC & + + 2d0*rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam + end do + end do + +end function diff --git a/src/GW/RGW_SRG_self_energy.f90 b/src/GW/RGW_SRG_self_energy.f90 index 4fcc44f..6ad1c90 100644 --- a/src/GW/RGW_SRG_self_energy.f90 +++ b/src/GW/RGW_SRG_self_energy.f90 @@ -1,4 +1,4 @@ -subroutine RGW_SRG_self_energy(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) +subroutine RGW_SRG_self_energy(flow,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! Compute correlation part of the self-energy @@ -7,6 +7,7 @@ subroutine RGW_SRG_self_energy(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! Input variables + double precision,intent(in) :: flow integer,intent(in) :: nBas integer,intent(in) :: nOrb integer,intent(in) :: nC @@ -36,7 +37,7 @@ subroutine RGW_SRG_self_energy(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! SRG flow parameter ! !--------------------! - s = 500d0 + s = flow !--------------------! ! SRG-GW self-energy ! diff --git a/src/GW/RGW_SRG_self_energy_diag.f90 b/src/GW/RGW_SRG_self_energy_diag.f90 index 41820e2..eef9f29 100644 --- a/src/GW/RGW_SRG_self_energy_diag.f90 +++ b/src/GW/RGW_SRG_self_energy_diag.f90 @@ -1,4 +1,4 @@ -subroutine RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) +subroutine RGW_SRG_self_energy_diag(flow,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! Compute correlation part of the self-energy @@ -7,6 +7,7 @@ subroutine RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC, ! Input variables + double precision,intent(in) :: flow integer,intent(in) :: nBas integer,intent(in) :: nOrb integer,intent(in) :: nC @@ -36,7 +37,7 @@ subroutine RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC, ! SRG flow parameter ! !--------------------! - s = 500d0 + s = flow !--------------------! ! SRG-GW self-energy ! diff --git a/src/GW/RGW_phACFDT.f90 b/src/GW/RGW_phACFDT.f90 index 3c7c18b..e2b2254 100644 --- a/src/GW/RGW_phACFDT.f90 +++ b/src/GW/RGW_phACFDT.f90 @@ -62,9 +62,9 @@ subroutine RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas, ! Hello World - write(*,*) '-------------------------------------------------------------' - write(*,*) ' Adiabatic connection version of BSE@G0W0 correlation energy ' - write(*,*) '-------------------------------------------------------------' + write(*,*) '-----------------------------------------------------------' + write(*,*) ' Adiabatic connection version of BSE@GW correlation energy ' + write(*,*) '-----------------------------------------------------------' write(*,*) ! eXtended BSE diff --git a/src/GW/RGW_phBSE.f90 b/src/GW/RGW_phBSE.f90 index 5953591..edb4111 100644 --- a/src/GW/RGW_phBSE.f90 +++ b/src/GW/RGW_phBSE.f90 @@ -190,5 +190,4 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple end if - end subroutine diff --git a/src/GW/UG0W0.f90 b/src/GW/UG0W0.f90 index 3643c89..0975f97 100644 --- a/src/GW/UG0W0.f90 +++ b/src/GW/UG0W0.f90 @@ -49,6 +49,7 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD logical :: dRPA integer :: is integer :: ispin + double precision :: flow double precision :: EcRPA double precision :: EcGM(nspin) double precision :: EcBSE(nspin) @@ -88,6 +89,8 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD ! SRG regularization + flow = 500d0 + if(doSRG) then write(*,*) '*** SRG regularized qsGW scheme ***' @@ -130,7 +133,7 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD !------------------------------------------------! if(doSRG) then - call UGW_SRG_self_energy_diag(nBas,nC,nO,nV,nR,nSt,eHF,Om,rho,SigC,Z,EcGM) + call UGW_SRG_self_energy_diag(flow,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 @@ -186,20 +189,9 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD if(dophBSE) then - call UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,S, & + call UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,S, & ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF,eGW,EcBSE) - if(exchange_kernel) then - - EcBSE(1) = 0.5d0*EcBSE(1) - EcBSE(2) = 0.5d0*EcBSE(2) - - else - - EcBSE(2) = 0.0d0 - - end if - write(*,*) write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy (spin-conserved) = ',EcBSE(1),' au' @@ -213,18 +205,6 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD if(doACFDT) then - write(*,*) '------------------------------------------------------------' - write(*,*) 'Adiabatic connection version of BSE@UG0W0 correlation energy' - write(*,*) '------------------------------------------------------------' - write(*,*) - - if(doXBS) then - - write(*,*) '*** scaled screening version (XBS) ***' - write(*,*) - - end if - call UGW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,dophBSE,spin_conserved,spin_flip,eta, & nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGW,EcBSE) diff --git a/src/GW/UGW_SRG_self_energy.f90 b/src/GW/UGW_SRG_self_energy.f90 index fc636d4..6dcd8d4 100644 --- a/src/GW/UGW_SRG_self_energy.f90 +++ b/src/GW/UGW_SRG_self_energy.f90 @@ -1,4 +1,4 @@ -subroutine UGW_SRG_self_energy(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) +subroutine UGW_SRG_self_energy(flow,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! Compute correlation part of the self-energy @@ -7,6 +7,7 @@ subroutine UGW_SRG_self_energy(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! Input variables + double precision,intent(in) :: flow integer,intent(in) :: nBas integer,intent(in) :: nC(nspin) integer,intent(in) :: nO(nspin) @@ -34,7 +35,7 @@ subroutine UGW_SRG_self_energy(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! SRG flow parameter - s = 500d0 + s = flow ! Initialize diff --git a/src/GW/UGW_SRG_self_energy_diag.f90 b/src/GW/UGW_SRG_self_energy_diag.f90 index ee4b725..e45148b 100644 --- a/src/GW/UGW_SRG_self_energy_diag.f90 +++ b/src/GW/UGW_SRG_self_energy_diag.f90 @@ -1,4 +1,4 @@ -subroutine UGW_SRG_self_energy_diag(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) +subroutine UGW_SRG_self_energy_diag(flow,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! Compute correlation part of the self-energy @@ -7,6 +7,7 @@ subroutine UGW_SRG_self_energy_diag(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! Input variables + double precision,intent(in) :: flow integer,intent(in) :: nBas integer,intent(in) :: nC(nspin) integer,intent(in) :: nO(nspin) @@ -34,7 +35,7 @@ subroutine UGW_SRG_self_energy_diag(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! SRG flow parameter - s = 500d0 + s = flow ! Initialize diff --git a/src/GW/UGW_phACFDT.f90 b/src/GW/UGW_phACFDT.f90 index 87e681d..ec7303a 100644 --- a/src/GW/UGW_phACFDT.f90 +++ b/src/GW/UGW_phACFDT.f90 @@ -63,6 +63,23 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s allocate(Ec(nAC,nspin)) + +! Hello World + + write(*,*) '-----------------------------------------------------------' + write(*,*) ' Adiabatic connection version of BSE@GW correlation energy ' + write(*,*) '-----------------------------------------------------------' + write(*,*) + +! eXtended BSE + + if(doXBS) then + + write(*,*) '*** scaled screening version (XBS) ***' + write(*,*) + + end if + ! Antisymmetrized kernel version if(exchange_kernel) then diff --git a/src/GW/UGW_phBSE.f90 b/src/GW/UGW_phBSE.f90 index 302a19b..85ca02b 100644 --- a/src/GW/UGW_phBSE.f90 +++ b/src/GW/UGW_phBSE.f90 @@ -1,4 +1,4 @@ -subroutine UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, & +subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, & S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cW,eW,eGW,EcBSE) ! Compute the Bethe-Salpeter excitation energies @@ -8,6 +8,7 @@ subroutine UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO ! Input variables + logical,intent(in) :: exchange_kernel logical,intent(in) :: TDA_W logical,intent(in) :: TDA logical,intent(in) :: dBSE @@ -161,4 +162,18 @@ subroutine UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO end if + +! Scale properly correlation energy if exchange is included in interaction kernel + + if(exchange_kernel) then + + EcBSE(1) = 0.5d0*EcBSE(1) + EcBSE(2) = 0.5d0*EcBSE(2) + + else + + EcBSE(2) = 0.0d0 + + end if + end subroutine diff --git a/src/GW/evGGW.f90 b/src/GW/evGGW.f90 index cc9717e..f8c59f0 100644 --- a/src/GW/evGGW.f90 +++ b/src/GW/evGGW.f90 @@ -41,10 +41,10 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! Local variables - logical :: linear_mixing logical :: dRPA = .true. integer :: nSCF integer :: n_diis + double precision :: flow double precision :: rcond double precision :: Conv double precision :: EcRPA @@ -81,6 +81,8 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! SRG regularization + flow = 500d0 + if(doSRG) then write(*,*) '*** SRG regularized qsGW scheme ***' @@ -88,11 +90,6 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop end if -! Linear mixing - - linear_mixing = .false. - alpha = 0.2d0 - ! Memory allocation allocate(Aph(nS,nS),Bph(nS,nS),eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas), & @@ -130,7 +127,7 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! Compute correlation part of the self-energy if(doSRG) then - call GGW_SRG_self_energy_diag(nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) + call GGW_SRG_self_energy_diag(flow,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) else call GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) end if @@ -163,18 +160,10 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! Linear mixing or DIIS extrapolation - if(linear_mixing) then + if(max_diis > 1) then - eGW(:) = alpha*eGW(:) + (1d0 - alpha)*eOld(:) - - else - n_diis = min(n_diis+1,max_diis) - if(abs(rcond) > 1d-7) then - call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGW-eOld,eGW) - else - n_diis = 0 - end if + call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGW-eOld,eGW) end if @@ -222,52 +211,8 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop write(*,*)'-------------------------------------------------------------------------------' write(*,*) -! Compute the BSE correlation energy via the adiabatic connection - -! if(doACFDT) then - -! write(*,*) '------------------------------------------------------' -! write(*,*) 'Adiabatic connection version of BSE correlation energy' -! write(*,*) '------------------------------------------------------' -! write(*,*) - -! if(doXBS) then - -! write(*,*) '*** scaled screening version (XBS) ***' -! write(*,*) - -! end if - -! call GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcBSE) - -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (singlet) =',EcBSE(1) -! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (triplet) =',EcBSE(2) -! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy =',EcBSE(1) + EcBSE(2) -! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW total energy =',ENuc + EGHF + EcBSE(1) + EcBSE(2) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) - -! end if - end if -! if(doppBSE) then - -! call GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) - -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy (singlet) =',EcBSE(1) -! write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy (triplet) =',3d0*EcBSE(2) -! write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy =',EcBSE(1) + 3d0*EcBSE(2) -! write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW total energy =',ENuc + EGHF + EcBSE(1) + 3d0*EcBSE(2) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) - -! end if - ! Testing zone if(dotest) then diff --git a/src/GW/evRGW.f90 b/src/GW/evRGW.f90 index f7e9258..438ffb4 100644 --- a/src/GW/evRGW.f90 +++ b/src/GW/evRGW.f90 @@ -44,11 +44,11 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! Local variables - logical :: linear_mixing logical :: dRPA = .true. integer :: ispin integer :: nSCF integer :: n_diis + double precision :: flow double precision :: rcond double precision :: Conv double precision :: EcRPA @@ -85,6 +85,8 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! SRG regularization + flow = 500d0 + if(doSRG) then write(*,*) '*** SRG regularized qsGW scheme ***' @@ -92,11 +94,6 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop end if -! Linear mixing - - linear_mixing = .false. - alpha = 0.2d0 - ! Memory allocation allocate(Aph(nS,nS),Bph(nS,nS),eGW(nOrb),eOld(nOrb),Z(nOrb),SigC(nOrb), & @@ -135,7 +132,7 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! Compute correlation part of the self-energy if(doSRG) then - call RGW_SRG_self_energy_diag(nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) + call RGW_SRG_self_energy_diag(flow,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 @@ -166,20 +163,12 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop call print_evRGW(nOrb,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM) - ! Linear mixing or DIIS extrapolation + ! DIIS extrapolation - if(linear_mixing) then - - eGW(:) = alpha*eGW(:) + (1d0 - alpha)*eOld(:) - - else + if(max_diis > 1) then n_diis = min(n_diis+1,max_diis) - if(abs(rcond) > 1d-7) then - call DIIS_extrapolation(rcond,nOrb,nOrb,n_diis,error_diis,e_diis,eGW-eOld,eGW) - else - n_diis = 0 - end if + call DIIS_extrapolation(rcond,nOrb,nOrb,n_diis,error_diis,e_diis,eGW-eOld,eGW) end if diff --git a/src/GW/evUGW.f90 b/src/GW/evUGW.f90 index b8a9085..80f4ce0 100644 --- a/src/GW/evUGW.f90 +++ b/src/GW/evUGW.f90 @@ -49,11 +49,11 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE ! Local variables logical :: dRPA - logical :: linear_mixing integer :: is integer :: ispin integer :: nSCF integer :: n_diis + double precision :: flow double precision :: rcond(nspin) double precision :: Conv double precision :: EcRPA(nspin) @@ -92,6 +92,8 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE ! SRG regularization + flow = 500d0 + if(doSRG) then write(*,*) '*** SRG regularized qsGW scheme ***' @@ -104,11 +106,6 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE EcRPA(:) = 0d0 dRPA = .true. -! Linear mixing - - linear_mixing = .false. - alpha = 0.2d0 - ! Memory allocation nSa = nS(1) @@ -156,7 +153,7 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE !------------------------------------------------! if(doSRG) then - call UGW_SRG_self_energy_diag(nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,SigC,Z,EcGM) + call UGW_SRG_self_energy_diag(flow,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 @@ -200,11 +197,7 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE ! Linear mixing or DIIS extrapolation - if(linear_mixing) then - - eGW(:,:) = alpha*eGW(:,:) + (1d0 - alpha)*eOld(:,:) - - else + if(max_diis > 1) then n_diis = min(n_diis+1,max_diis) do is=1,nspin @@ -212,10 +205,6 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE e_diis(:,1:n_diis,is),eGW(:,is)-eOld(:,is),eGW(:,is)) end do -! Reset DIIS if required - - if(minval(rcond(:)) < 1d-15) n_diis = 0 - end if ! Save quasiparticles energy for next cycle @@ -253,20 +242,9 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE if(BSE) then - call UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, & + call UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, & S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eGW,eGW,EcBSE) - if(exchange_kernel) then - - EcBSE(1) = 0.5d0*EcBSE(1) - EcBSE(2) = 0.5d0*EcBSE(2) - - else - - EcBSE(2) = 0.0d0 - - end if - write(*,*) write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@evGW@UHF correlation energy (spin-conserved) =',EcBSE(1),' au' @@ -280,18 +258,6 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE if(doACFDT) then - write(*,*) '--------------------------------------------------------------' - write(*,*) ' Adiabatic connection version of BSE@evUGW correlation energy ' - write(*,*) '--------------------------------------------------------------' - write(*,*) - - if(doXBS) then - - write(*,*) '*** scaled screening version (XBS) ***' - write(*,*) - - end if - call UGW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,spin_conserved,spin_flip, & eta,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eGW,eGW,EcRPA) diff --git a/src/GW/qsGGW.f90 b/src/GW/qsGGW.f90 index 0f32c34..49a214a 100644 --- a/src/GW/qsGGW.f90 +++ b/src/GW/qsGGW.f90 @@ -60,6 +60,7 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop integer :: nBas2Sq integer :: ixyz integer :: n_diis + double precision :: flow double precision :: ET,ETaa,ETbb double precision :: EV,EVaa,EVbb double precision :: EJ,EJaaaa,EJaabb,EJbbaa,EJbbbb @@ -132,6 +133,8 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! SRG regularization + flow = 500d0 + if(doSRG) then write(*,*) '*** SRG regularized qsGW scheme ***' @@ -264,7 +267,7 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop call GGW_excitation_density(nBas2,nC,nO,nR,nS,ERI_MO,XpY,rho) if(doSRG) then - call GGW_SRG_self_energy(nBas2,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) + call GGW_SRG_self_energy(flow,nBas2,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) else call GGW_self_energy(eta,nBas2,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) end if @@ -395,52 +398,8 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop write(*,*)'-------------------------------------------------------------------------------' write(*,*) -! Compute the BSE correlation energy via the adiabatic connection - -! if(doACFDT) then - -! write(*,*) '------------------------------------------------------' -! write(*,*) 'Adiabatic connection version of BSE correlation energy' -! write(*,*) '------------------------------------------------------' -! write(*,*) - -! if(doXBS) then - -! write(*,*) '*** scaled screening version (XBS) ***' -! write(*,*) - -! end if - -! call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcBSE) - -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy (singlet) =',EcBSE(1) -! write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy (triplet) =',EcBSE(2) -! write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy =',EcBSE(1) + EcBSE(2) -! write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW total energy =',ENuc + EqsGW + EcBSE(1) + EcBSE(2) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) - -! end if - end if -! if(doppBSE) then -! -! call GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eHF,eGW,EcBSE) -! -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@qsGW correlation energy (singlet) =',EcBSE(1) -! write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@qsGW correlation energy (triplet) =',3d0*EcBSE(2) -! write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@qsGW correlation energy =',EcBSE(1) + 3d0*EcBSE(2) -! write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@qsGW total energy =',ENuc + EGHF + EcBSE(1) + 3d0*EcBSE(2) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) - -! end if - ! Testing zone if(dotest) then diff --git a/src/GW/qsRGW.f90 b/src/GW/qsRGW.f90 index fe3b9b8..a8e4ce2 100644 --- a/src/GW/qsRGW.f90 +++ b/src/GW/qsRGW.f90 @@ -76,6 +76,8 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop double precision,external :: trace_matrix double precision :: dipole(ncart) + double precision :: flow + logical :: dRPA_W = .true. logical :: print_W = .false. double precision,allocatable :: err_diis(:,:) @@ -125,6 +127,8 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! SRG regularization + flow = 500d0 + if(doSRG) then write(*,*) '*** SRG regularized qsGW scheme ***' @@ -207,7 +211,7 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI_MO,XpY,rho) if(doSRG) then - call RGW_SRG_self_energy(nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) + call RGW_SRG_self_energy(flow,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 diff --git a/src/GW/qsUGW.f90 b/src/GW/qsUGW.f90 index f9ebb50..4fff640 100644 --- a/src/GW/qsUGW.f90 +++ b/src/GW/qsUGW.f90 @@ -1,4 +1,4 @@ -subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, & +subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, & 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) @@ -17,7 +17,7 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE logical,intent(in) :: doACFDT logical,intent(in) :: exchange_kernel logical,intent(in) :: doXBS - logical,intent(in) :: BSE + logical,intent(in) :: dophBSE logical,intent(in) :: TDA_W logical,intent(in) :: TDA logical,intent(in) :: dBSE @@ -66,6 +66,7 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE integer :: is integer :: n_diis integer :: nSa,nSb,nSt + double precision :: flow double precision :: dipole(ncart) double precision :: ET(nspin) @@ -128,6 +129,8 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE ! SRG regularization + flow = 500d0 + if(doSRG) then write(*,*) '*** SRG regularized qsGW scheme ***' @@ -220,7 +223,7 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE !------------------------------------------------! if(doSRG) then - call UGW_SRG_self_energy(nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,SigC,Z,EcGM) + call UGW_SRG_self_energy(flow,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 @@ -358,22 +361,11 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE ! Perform BSE calculation - if(BSE) then + if(dophBSE) then - call UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, & + call UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, & S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,c,eGW,eGW,EcBSE) - if(exchange_kernel) then - - EcBSE(1) = 0.5d0*EcBSE(1) - EcBSE(2) = 0.5d0*EcBSE(2) - - else - - EcBSE(2) = 0.0d0 - - end if - write(*,*) write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@qsGW@UHF correlation energy (spin-conserved) = ',EcBSE(1),' au' @@ -387,19 +379,7 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE if(doACFDT) then - write(*,*) '--------------------------------------------------------------' - write(*,*) ' Adiabatic connection version of BSE@qsUGW correlation energy ' - write(*,*) '--------------------------------------------------------------' - write(*,*) - - if(doXBS) then - - write(*,*) '*** scaled screening version (XBS) ***' - write(*,*) - - end if - - call UGW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,spin_conserved,spin_flip, & + call UGW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,dophBSE,spin_conserved,spin_flip, & eta,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eGW,eGW,EcRPA) write(*,*)