From 32998c8a3256b630223d14576aac508acca9dccf Mon Sep 17 00:00:00 2001 From: pfloos Date: Sat, 14 Sep 2024 23:25:38 +0200 Subject: [PATCH] symmetrized upfolded phBSE scheme --- src/GT/RGTpp_phBSE_dynamic_perturbation.f90 | 5 - src/GW/RGW_phBSE.f90 | 2 + src/GW/RGW_phBSE_dynamic_perturbation.f90 | 5 - src/GW/RGW_phBSE_upfolded.f90 | 20 +- src/GW/RGW_phBSE_upfolded_sym.f90 | 259 ++++++++++++++++++++ src/GW/UGW_phBSE_dynamic_perturbation.f90 | 5 - 6 files changed, 271 insertions(+), 25 deletions(-) create mode 100644 src/GW/RGW_phBSE_upfolded_sym.f90 diff --git a/src/GT/RGTpp_phBSE_dynamic_perturbation.f90 b/src/GT/RGTpp_phBSE_dynamic_perturbation.f90 index bb8dace..04835e9 100644 --- a/src/GT/RGTpp_phBSE_dynamic_perturbation.f90 +++ b/src/GT/RGTpp_phBSE_dynamic_perturbation.f90 @@ -47,7 +47,6 @@ subroutine RGTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,n integer :: ia integer :: maxS = 10 - double precision :: gapGT double precision,allocatable :: OmDyn(:) double precision,allocatable :: ZDyn(:) @@ -111,13 +110,9 @@ subroutine RGTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,n ! Dump results ! !--------------! - gapGT = eGT(nO+1) - eGT(nO) - write(*,*) '---------------------------------------------------------------------------------------------------' write(*,*) ' First-order dynamical correction to static Bethe-Salpeter excitation energies ' write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,'(A57,F10.6,A3)') ' BSE neutral excitation must be lower than the GT gap = ',gapGT*HaToeV,' eV' - write(*,*) '---------------------------------------------------------------------------------------------------' write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)' write(*,*) '---------------------------------------------------------------------------------------------------' diff --git a/src/GW/RGW_phBSE.f90 b/src/GW/RGW_phBSE.f90 index 1438194..24b7eed 100644 --- a/src/GW/RGW_phBSE.f90 +++ b/src/GW/RGW_phBSE.f90 @@ -146,6 +146,8 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple ! Upfolded phBSE ! !----------------! + call RGW_phBSE_upfolded_sym(ispin,nOrb,nOrb,nC,nO,nV,nR,nS,ERI,rho_RPA,OmRPA,eW) + call RGW_phBSE_upfolded(ispin,nOrb,nOrb,nC,nO,nV,nR,nS,ERI,rho_RPA,OmRPA,eGW) end if diff --git a/src/GW/RGW_phBSE_dynamic_perturbation.f90 b/src/GW/RGW_phBSE_dynamic_perturbation.f90 index 57342c8..b94897e 100644 --- a/src/GW/RGW_phBSE_dynamic_perturbation.f90 +++ b/src/GW/RGW_phBSE_dynamic_perturbation.f90 @@ -35,7 +35,6 @@ subroutine RGW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS, integer :: ia integer :: maxS = 10 - double precision :: gapGW double precision,allocatable :: Om_dyn(:) double precision,allocatable :: Z_dyn(:) @@ -73,13 +72,9 @@ subroutine RGW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS, end if - gapGW = eGW(nO+1) - eGW(nO) - write(*,*) '---------------------------------------------------------------------------------------------------' write(*,*) ' First-order dynamical correction to static Bethe-Salpeter excitation energies ' write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,'(A57,F10.6,A3)') ' BSE neutral excitation must be lower than the GW gap = ',gapGW*HaToeV,' eV' - write(*,*) '---------------------------------------------------------------------------------------------------' write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)' write(*,*) '---------------------------------------------------------------------------------------------------' diff --git a/src/GW/RGW_phBSE_upfolded.f90 b/src/GW/RGW_phBSE_upfolded.f90 index 6bb1342..91d48af 100644 --- a/src/GW/RGW_phBSE_upfolded.f90 +++ b/src/GW/RGW_phBSE_upfolded.f90 @@ -27,7 +27,7 @@ subroutine RGW_phBSE_upfolded(ispin,nBas,nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW) integer :: a,b,c,d integer :: m,ia,jb,iam,kcm integer,parameter :: maxH = 20 - double precision :: tmp + double precision :: Jph,Kph,C2h2p integer :: n1h,n1p,n1h1p,n2h2p,nH double precision,external :: Kronecker_delta @@ -153,15 +153,15 @@ subroutine RGW_phBSE_upfolded(ispin,nBas,nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW) ! Jph - tmp = - sqrt(2d0)*Kronecker_delta(a,c)*rho(i,k,m) - H(ia ,n1h1p+kcm) = tmp - H(n1h1p+n2h2p+kcm,ia) = tmp + Jph = - sqrt(2d0)*Kronecker_delta(a,c)*rho(i,k,m) + H(ia ,n1h1p+kcm) = Jph + H(n1h1p+n2h2p+kcm,ia) = Jph ! Kph - tmp = sqrt(2d0)*Kronecker_delta(i,k)*rho(a,c,m) - H(ia ,n1h1p+n2h2p+kcm) = tmp - H(n1h1p+kcm,ia) = tmp + Kph = + sqrt(2d0)*Kronecker_delta(i,k)*rho(a,c,m) + H(ia ,n1h1p+n2h2p+kcm) = Kph + H(n1h1p+kcm,ia) = Kph end do end do @@ -180,9 +180,9 @@ subroutine RGW_phBSE_upfolded(ispin,nBas,nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW) do m=1,nS iam = iam + 1 - tmp = Om(m) + eGW(a) - eGW(i) - H(n1h1p +iam,n1h1p +iam) = tmp - H(n1h1p+n2h2p+iam,n1h1p+n2h2p+iam) = tmp + C2h2p = Om(m) + eGW(a) - eGW(i) + H(n1h1p +iam,n1h1p +iam) = C2h2p + H(n1h1p+n2h2p+iam,n1h1p+n2h2p+iam) = C2h2p end do end do diff --git a/src/GW/RGW_phBSE_upfolded_sym.f90 b/src/GW/RGW_phBSE_upfolded_sym.f90 new file mode 100644 index 0000000..362d17b --- /dev/null +++ b/src/GW/RGW_phBSE_upfolded_sym.f90 @@ -0,0 +1,259 @@ +subroutine RGW_phBSE_upfolded_sym(ispin,nBas,nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eHF) + +! Upfolded phBSE@GW (TDA only!) + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + 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) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: rho(nOrb,nOrb,nS) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: eHF(nOrb) + +! Local variables + + integer :: s + integer :: i,j,k,l + integer :: a,b,c,d + integer :: m,ia,jb,iam,kcm + integer,parameter :: maxH = 20 + double precision :: Jph,Kph,C2h2p + + integer :: n1h,n1p,n1h1p,n2h2p,nH + double precision,external :: Kronecker_delta + + integer,allocatable :: order(:) + double precision,allocatable :: H(:,:) + double precision,allocatable :: OmBSE(:) + double precision,allocatable :: Z(:) + +! Output variables + +! Hello world + + write(*,*) + write(*,*)'*********************************************' + write(*,*)'* Symmetrized Upfolded phBSE@GW Calculation *' + write(*,*)'*********************************************' + write(*,*) + +! TDA for W + + write(*,*) 'Tamm-Dancoff approximation by default!' + write(*,*) + +! Dimension of the supermatrix + + n1h = nO + n1p = nV + + n1h1p = n1h*n1p + n2h2p = n1h1p*n1h1p + nH = n1h1p + n2h2p + +! Memory allocation + + allocate(order(nH),H(nH,nH),OmBSE(nH),Z(nH)) + +! Initialization + + H(:,:) = 0d0 + +!-------------------------------! +! Compute BSE supermatrix ! +!-------------------------------! +! ! +! | A Jph + Kph | ! +! H = | | ! +! | Jph + Kph C2h2p | ! +! ! +!-------------------------------! + + !----------------------! + ! Block A for singlets ! + !----------------------! + + if(ispin == 1) then + + ia = 0 + do i=nC+1,nO + do a=nO+1,nOrb-nR + ia = ia + 1 + + jb = 0 + do j=nC+1,nO + do b=nO+1,nOrb-nR + jb = jb + 1 + + H(ia,jb) = (eHF(a) - eHF(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + + 2d0*ERI(i,b,a,j) - ERI(i,b,j,a) + + do k=nC+1,nO + do m=1,nS + + H(ia,jb) = H(ia,jb) & + + Kronecker_delta(i,j)*rho(a,k,m)*rho(b,k,m)/(eHF(a) - eHF(k) + Om(m)) & + + Kronecker_delta(i,j)*rho(a,k,m)*rho(b,k,m)/(eHF(b) - eHF(k) + Om(m)) + + end do + end do + + do c=nO+1,nOrb-nR + do m=1,nS + + H(ia,jb) = H(ia,jb) & + - Kronecker_delta(a,b)*rho(i,c,m)*rho(j,c,m)/(eHF(i) - eHF(c) - Om(m)) & + - Kronecker_delta(a,b)*rho(i,c,m)*rho(j,c,m)/(eHF(j) - eHF(c) - Om(m)) + + end do + end do + + end do + end do + + end do + end do + + end if + + !----------------------! + ! Block A for triplets ! + !----------------------! + + if(ispin == 2) then + + ia = 0 + do i=nC+1,nO + do a=nO+1,nOrb-nR + ia = ia + 1 + + jb = 0 + do j=nC+1,nO + do b=nO+1,nOrb-nR + jb = jb + 1 + + H(ia,jb) = (eHF(a) - eHF(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + - ERI(i,b,j,a) + + do k=nC+1,nO + do m=1,nS + + H(ia,jb) = H(ia,jb) & + + Kronecker_delta(i,j)*rho(a,k,m)*rho(b,k,m)/(eHF(a) - eHF(k) + Om(m)) & + + Kronecker_delta(i,j)*rho(a,k,m)*rho(b,k,m)/(eHF(b) - eHF(k) + Om(m)) + + end do + end do + + do c=nO+1,nOrb-nR + do m=1,nS + + H(ia,jb) = H(ia,jb) & + - Kronecker_delta(a,b)*rho(i,c,m)*rho(j,c,m)/(eHF(i) - eHF(c) - Om(m)) & + - Kronecker_delta(a,b)*rho(i,c,m)*rho(j,c,m)/(eHF(j) - eHF(c) - Om(m)) + + end do + end do + + end do + end do + + end do + end do + + end if + + !------------------! + ! Blocks Jph & Kph ! + !------------------! + + ia = 0 + do i=nC+1,nO + do a=nO+1,nOrb-nR + ia = ia + 1 + + kcm = 0 + do k=nC+1,nO + do c=nO+1,nOrb-nR + do m=1,nS + kcm = kcm + 1 + + ! Jph + Kph + + Jph = - sqrt(2d0)*Kronecker_delta(a,c)*rho(i,k,m) + Kph = + sqrt(2d0)*Kronecker_delta(i,k)*rho(a,c,m) + H(ia ,n1h1p+kcm) = Jph + Kph + H(n1h1p+kcm,ia) = Jph + Kph + + end do + end do + end do + + end do + end do + + !-------------! + ! Block 2h2p ! + !-------------! + + iam = 0 + do i=nC+1,nO + do a=nO+1,nOrb-nR + do m=1,nS + iam = iam + 1 + + C2h2p = Om(m) + eHF(a) - eHF(i) + H(n1h1p+iam,n1h1p+iam) = C2h2p + + end do + end do + end do + +!-------------------------! +! Diagonalize supermatrix ! +!-------------------------! + + call diagonalize_matrix(nH,H,OmBSE) + +!-----------------! +! Compute weights ! +!-----------------! + + Z(:) = 0d0 + do s=1,nH + do ia=1,n1h1p + Z(s) = Z(s) + H(ia,s)**2 + end do + end do + +!--------------! +! Dump results ! +!--------------! + + write(*,*)'-------------------------------------------' + write(*,*)' Upfolded phBSE excitation energies (eV) ' + write(*,*)'-------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') & + '|','#','|','OmBSE (eV)','|','Z','|' + write(*,*)'-------------------------------------------' + + do s=1,min(nH,maxH) + if(Z(s) > 1d-7) & + write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & + '|',s,'|',OmBSE(s)*HaToeV,'|',Z(s),'|' + end do + + write(*,*)'-------------------------------------------' + write(*,*) + +end subroutine diff --git a/src/GW/UGW_phBSE_dynamic_perturbation.f90 b/src/GW/UGW_phBSE_dynamic_perturbation.f90 index 75655c9..65bf043 100644 --- a/src/GW/UGW_phBSE_dynamic_perturbation.f90 +++ b/src/GW/UGW_phBSE_dynamic_perturbation.f90 @@ -41,7 +41,6 @@ subroutine UGW_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nSa integer :: ia integer,parameter :: maxS = 10 - double precision :: gapGW double precision,allocatable :: OmDyn(:) double precision,allocatable :: ZDyn(:) @@ -63,13 +62,9 @@ subroutine UGW_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nSa write(*,*) end if - gapGW = min(eGW(nO(1)+1,1),eGW(nO(2)+1,2)) - max(eGW(nO(1),1),eGW(nO(2),2)) - write(*,*) '---------------------------------------------------------------------------------------------------' write(*,*) ' First-order dynamical correction to static Bethe-Salpeter excitation energies ' write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,'(A57,F10.6,A3)') ' BSE neutral excitation must be lower than the GW gap = ',gapGW*HaToeV,' eV' - write(*,*) '---------------------------------------------------------------------------------------------------' write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)' write(*,*) '---------------------------------------------------------------------------------------------------'