mirror of
https://github.com/pfloos/quack
synced 2024-12-22 04:14:26 +01:00
GGW ppBSE dyn
This commit is contained in:
Pierre-Francois Loos
parent
dc0f77d525
commit
88b0424cd3
@ -111,9 +111,9 @@ subroutine GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,
|
||||
! Compute the dynamical screening at the ppBSE level !
|
||||
!----------------------------------------------------!
|
||||
|
||||
! if(dBSE) &
|
||||
! call GGW_ppBSE_dynamic_perturbation(dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, &
|
||||
! Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta)
|
||||
if(dBSE) &
|
||||
call GGW_ppBSE_dynamic_perturbation(dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, &
|
||||
Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta)
|
||||
|
||||
|
||||
deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta)
|
||||
|
||||
82
src/GW/GGW_ppBSE_dynamic_kernel_B.f90
Normal file
82
src/GW/GGW_ppBSE_dynamic_kernel_B.f90
Normal file
@ -0,0 +1,82 @@
|
||||
subroutine GGW_ppBSE_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGW,Om,rho,KB_dyn)
|
||||
|
||||
! Compute the dynamic part of the Bethe-Salpeter equation matrices
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nS
|
||||
integer,intent(in) :: nOO
|
||||
integer,intent(in) :: nVV
|
||||
double precision,intent(in) :: eta
|
||||
double precision,intent(in) :: lambda
|
||||
double precision,intent(in) :: eGW(nBas)
|
||||
double precision,intent(in) :: Om(nS)
|
||||
double precision,intent(in) :: rho(nBas,nBas,nS)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision,external :: Kronecker_delta
|
||||
double precision :: dem,num
|
||||
integer :: m
|
||||
integer :: a,b,i,j
|
||||
integer :: ab,ij
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: KB_dyn(nVV,nOO)
|
||||
|
||||
! Initialization
|
||||
|
||||
KB_dyn(:,:) = 0d0
|
||||
|
||||
! Build dynamic B matrix
|
||||
|
||||
ab = 0
|
||||
do a=nO+1,nBas-nR
|
||||
do b=a+1,nBas-nR
|
||||
ab = ab + 1
|
||||
|
||||
ij = 0
|
||||
do i=nC+1,nO
|
||||
do j=i+1,nO
|
||||
ij = ij + 1
|
||||
|
||||
do m=1,nS
|
||||
|
||||
dem = eGW(j) - Om(m) - eGW(b)
|
||||
num = rho(a,i,m)*rho(b,j,m)
|
||||
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
dem = eGW(j) - Om(m) - eGW(a)
|
||||
num = rho(b,i,m)*rho(a,j,m)
|
||||
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij) - num*dem/(dem**2 + eta**2)
|
||||
|
||||
dem = eGW(i) - Om(m) - eGW(a)
|
||||
num = rho(a,i,m)*rho(b,j,m)
|
||||
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2)
|
||||
|
||||
dem = eGW(i) - Om(m) - eGW(b)
|
||||
num = rho(b,i,m)*rho(a,j,m)
|
||||
|
||||
KB_dyn(ab,ij) = KB_dyn(ab,ij) - num*dem/(dem**2 + eta**2)
|
||||
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
end subroutine
|
||||
88
src/GW/GGW_ppBSE_dynamic_kernel_C.f90
Normal file
88
src/GW/GGW_ppBSE_dynamic_kernel_C.f90
Normal file
@ -0,0 +1,88 @@
|
||||
subroutine GGW_ppBSE_dynamic_kernel_C(eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,eGW,Om,rho,OmBSE,KC_dyn,ZC_dyn)
|
||||
|
||||
! Compute the dynamic part of the Bethe-Salpeter equation matrices
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nS
|
||||
integer,intent(in) :: nVV
|
||||
double precision,intent(in) :: eta
|
||||
double precision,intent(in) :: lambda
|
||||
double precision,intent(in) :: eGW(nBas)
|
||||
double precision,intent(in) :: Om(nS)
|
||||
double precision,intent(in) :: rho(nBas,nBas,nS)
|
||||
double precision,intent(in) :: OmBSE
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision,external :: Kronecker_delta
|
||||
double precision :: dem,num
|
||||
integer :: m
|
||||
integer :: a,b,c,d
|
||||
integer :: ab,cd
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: KC_dyn(nVV,nVV)
|
||||
double precision,intent(out) :: ZC_dyn(nVV,nVV)
|
||||
|
||||
! Initialization
|
||||
|
||||
KC_dyn(:,:) = 0d0
|
||||
ZC_dyn(:,:) = 0d0
|
||||
|
||||
! Build dynamic C matrix
|
||||
|
||||
ab = 0
|
||||
do a=nO+1,nBas-nR
|
||||
do b=a+1,nBas-nR
|
||||
ab = ab + 1
|
||||
|
||||
cd = 0
|
||||
do c=nO+1,nBas-nR
|
||||
do d=c+1,nBas-nR
|
||||
cd = cd + 1
|
||||
|
||||
do m=1,nS
|
||||
|
||||
dem = OmBSE - eGW(c) - Om(m) - eGW(b)
|
||||
num = rho(a,c,m)*rho(b,d,m)
|
||||
|
||||
KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2)
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
dem = OmBSE - eGW(c) - Om(m) - eGW(a)
|
||||
num = rho(b,c,m)*rho(a,d,m)
|
||||
|
||||
KC_dyn(ab,cd) = KC_dyn(ab,cd) - num*dem/(dem**2 + eta**2)
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
dem = OmBSE - eGW(d) - Om(m) - eGW(a)
|
||||
num = rho(a,c,m)*rho(b,d,m)
|
||||
|
||||
KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2)
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
dem = OmBSE - eGW(d) - Om(m) - eGW(b)
|
||||
num = rho(b,c,m)*rho(a,d,m)
|
||||
|
||||
KC_dyn(ab,cd) = KC_dyn(ab,cd) - num*dem/(dem**2 + eta**2)
|
||||
ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
end subroutine
|
||||
88
src/GW/GGW_ppBSE_dynamic_kernel_D.f90
Normal file
88
src/GW/GGW_ppBSE_dynamic_kernel_D.f90
Normal file
@ -0,0 +1,88 @@
|
||||
subroutine GGW_ppBSE_dynamic_kernel_D(eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,eGW,Om,rho,OmBSE,KD_dyn,ZD_dyn)
|
||||
|
||||
! Compute the dynamic part of the Bethe-Salpeter equation matrices
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nS
|
||||
integer,intent(in) :: nOO
|
||||
double precision,intent(in) :: eta
|
||||
double precision,intent(in) :: lambda
|
||||
double precision,intent(in) :: eGW(nBas)
|
||||
double precision,intent(in) :: Om(nS)
|
||||
double precision,intent(in) :: rho(nBas,nBas,nS)
|
||||
double precision,intent(in) :: OmBSE
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision,external :: Kronecker_delta
|
||||
double precision :: dem,num
|
||||
integer :: m
|
||||
integer :: i,j,k,l
|
||||
integer :: ij,kl
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: KD_dyn(nOO,nOO)
|
||||
double precision,intent(out) :: ZD_dyn(nOO,nOO)
|
||||
|
||||
! Initialization
|
||||
|
||||
KD_dyn(:,:) = 0d0
|
||||
ZD_dyn(:,:) = 0d0
|
||||
|
||||
! Build dynamic D matrix
|
||||
|
||||
ij = 0
|
||||
do i=nC+1,nO
|
||||
do j=i+1,nO
|
||||
ij = ij + 1
|
||||
|
||||
kl = 0
|
||||
do k=nC+1,nO
|
||||
do l=k+1,nO
|
||||
kl = kl + 1
|
||||
|
||||
do m=1,nS
|
||||
|
||||
dem = OmBSE - eGW(k) + Om(m) - eGW(i)
|
||||
num = rho(i,k,m)*rho(j,l,m)
|
||||
|
||||
KD_dyn(ij,kl) = KD_dyn(ij,kl) - num*dem/(dem**2 + eta**2)
|
||||
ZD_dyn(ij,kl) = ZD_dyn(ij,kl) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
dem = OmBSE - eGW(k) + Om(m) - eGW(j)
|
||||
num = rho(j,k,m)*rho(i,l,m)
|
||||
|
||||
KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2)
|
||||
ZD_dyn(ij,kl) = ZD_dyn(ij,kl) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
dem = OmBSE - eGW(l) + Om(m) - eGW(j)
|
||||
num = rho(i,k,m)*rho(j,l,m)
|
||||
|
||||
KD_dyn(ij,kl) = KD_dyn(ij,kl) - num*dem/(dem**2 + eta**2)
|
||||
ZD_dyn(ij,kl) = ZD_dyn(ij,kl) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
dem = OmBSE - eGW(l) + Om(m) - eGW(i)
|
||||
num = rho(j,k,m)*rho(i,l,m)
|
||||
|
||||
KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2)
|
||||
ZD_dyn(ij,kl) = ZD_dyn(ij,kl) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
|
||||
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
end subroutine
|
||||
152
src/GW/GGW_ppBSE_dynamic_perturbation.f90
Normal file
152
src/GW/GGW_ppBSE_dynamic_perturbation.f90
Normal file
@ -0,0 +1,152 @@
|
||||
subroutine GGW_ppBSE_dynamic_perturbation(dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int, &
|
||||
OmRPA,rho_RPA,Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta)
|
||||
|
||||
! Compute dynamical effects via perturbation theory for BSE
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
logical,intent(in) :: dTDA
|
||||
double precision,intent(in) :: eta
|
||||
integer,intent(in) :: nOrb
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nS
|
||||
integer,intent(in) :: nOO
|
||||
integer,intent(in) :: nVV
|
||||
|
||||
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
|
||||
double precision,intent(in) :: eW(nOrb)
|
||||
double precision,intent(in) :: eGW(nOrb)
|
||||
double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart)
|
||||
double precision,intent(in) :: OmRPA(nS)
|
||||
double precision,intent(in) :: rho_RPA(nOrb,nOrb,nS)
|
||||
double precision,intent(in) :: Om1(nVV)
|
||||
double precision,intent(in) :: X1(nVV,nVV)
|
||||
double precision,intent(in) :: Y1(nOO,nVV)
|
||||
double precision,intent(in) :: Om2(nOO)
|
||||
double precision,intent(in) :: X2(nVV,nOO)
|
||||
double precision,intent(in) :: Y2(nOO,nOO)
|
||||
double precision,intent(in) :: KB_sta(nVV,nOO)
|
||||
double precision,intent(in) :: KC_sta(nVV,nVV)
|
||||
double precision,intent(in) :: KD_sta(nOO,nOO)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: ab,ij,kl
|
||||
|
||||
integer :: maxOO = 10
|
||||
integer :: maxVV = 0
|
||||
|
||||
double precision,allocatable :: Om1_dyn(:)
|
||||
double precision,allocatable :: Om2_dyn(:)
|
||||
double precision,allocatable :: Z1_dyn(:)
|
||||
double precision,allocatable :: Z2_dyn(:)
|
||||
|
||||
double precision,allocatable :: KB_dyn(:,:)
|
||||
double precision,allocatable :: KC_dyn(:,:)
|
||||
double precision,allocatable :: KD_dyn(:,:)
|
||||
double precision,allocatable :: ZC_dyn(:,:)
|
||||
double precision,allocatable :: ZD_dyn(:,:)
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(Om1_dyn(maxVV),Om2_dyn(maxOO),Z1_dyn(maxVV),Z2_dyn(maxOO), &
|
||||
KB_dyn(nVV,nOO),KC_dyn(nVV,nVV),KD_dyn(nOO,nOO), &
|
||||
ZC_dyn(nVV,nVV),ZD_dyn(nOO,nOO))
|
||||
|
||||
if(dTDA) then
|
||||
write(*,*)
|
||||
write(*,*) '*** dynamical TDA activated ***'
|
||||
write(*,*)
|
||||
end if
|
||||
|
||||
write(*,*) '---------------------------------------------------------------------------------------------------'
|
||||
write(*,*) ' First-order dynamical correction to static ppBSE double electron attachment energies '
|
||||
write(*,*) '---------------------------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)'
|
||||
write(*,*) '---------------------------------------------------------------------------------------------------'
|
||||
|
||||
do ab=1,min(nVV,maxVV)
|
||||
|
||||
if(dTDA) then
|
||||
|
||||
call GGW_ppBSE_dynamic_kernel_C(eta,nOrb,nC,nO,nV,nR,nS,nVV,1d0,eGW,OmRPA,rho_RPA,Om1(ab),KC_dyn,ZC_dyn)
|
||||
|
||||
Z1_dyn(ab) = + dot_product(X1(:,ab),matmul(ZC_dyn,X1(:,ab)))
|
||||
Om1_dyn(ab) = + dot_product(X1(:,ab),matmul(KC_dyn - KC_sta,X1(:,ab)))
|
||||
|
||||
else
|
||||
|
||||
call GGW_ppBSE_dynamic_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGW,OmRPA,rho_RPA,KB_dyn)
|
||||
call GGW_ppBSE_dynamic_kernel_C(eta,nOrb,nC,nO,nV,nR,nS,nVV,1d0,eGW,OmRPA,rho_RPA,Om1(ab),KC_dyn,ZC_dyn)
|
||||
call GGW_ppBSE_dynamic_kernel_D(eta,nOrb,nC,nO,nV,nR,nS,nOO,1d0,eGW,OmRPA,rho_RPA,Om1(ab),KD_dyn,ZD_dyn)
|
||||
|
||||
Z1_dyn(ab) = dot_product(X1(:,ab),matmul(ZC_dyn,X1(:,ab))) &
|
||||
+ dot_product(Y1(:,ab),matmul(ZD_dyn,Y1(:,ab)))
|
||||
|
||||
Om1_dyn(ab) = dot_product(X1(:,ab),matmul(KC_dyn - KC_sta,X1(:,ab))) &
|
||||
- dot_product(Y1(:,ab),matmul(KD_dyn - KD_sta,Y1(:,ab))) &
|
||||
+ dot_product(X1(:,ab),matmul(KB_dyn - KB_sta,Y1(:,ab))) &
|
||||
- dot_product(Y1(:,ab),matmul(transpose(KB_dyn - KB_sta),X1(:,ab)))
|
||||
|
||||
end if
|
||||
|
||||
Z1_dyn(ab) = 1d0/(1d0 - Z1_dyn(ab))
|
||||
Om1_dyn(ab) = Z1_dyn(ab)*Om1_dyn(ab)
|
||||
|
||||
write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') &
|
||||
ab,Om1(ab)*HaToeV,(Om1(ab)+Om1_dyn(ab))*HaToeV,Om1_dyn(ab)*HaToeV,Z1_dyn(ab)
|
||||
|
||||
end do
|
||||
write(*,*) '---------------------------------------------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
write(*,*) '---------------------------------------------------------------------------------------------------'
|
||||
write(*,*) ' First-order dynamical correction to static ppBSE double electron detachment energies '
|
||||
write(*,*) '---------------------------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)'
|
||||
write(*,*) '---------------------------------------------------------------------------------------------------'
|
||||
|
||||
kl = 0
|
||||
do ij=nOO,max(1,nOO+1-maxOO),-1
|
||||
kl = kl + 1
|
||||
|
||||
if(dTDA) then
|
||||
|
||||
call GGW_ppBSE_dynamic_kernel_D(eta,nOrb,nC,nO,nV,nR,nS,nOO,1d0,eGW,OmRPA,rho_RPA,Om2(ij),KD_dyn,ZD_dyn)
|
||||
|
||||
Z2_dyn(kl) = - dot_product(Y2(:,ij),matmul(ZD_dyn,Y2(:,ij)))
|
||||
Om2_dyn(kl) = - dot_product(Y2(:,ij),matmul(KD_dyn - KD_sta,Y2(:,ij)))
|
||||
|
||||
else
|
||||
|
||||
call GGW_ppBSE_dynamic_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGW,OmRPA,rho_RPA,KB_dyn)
|
||||
call GGW_ppBSE_dynamic_kernel_C(eta,nOrb,nC,nO,nV,nR,nS,nVV,1d0,eGW,OmRPA,rho_RPA,Om2(ij),KC_dyn,ZC_dyn)
|
||||
call GGW_ppBSE_dynamic_kernel_D(eta,nOrb,nC,nO,nV,nR,nS,nOO,1d0,eGW,OmRPA,rho_RPA,Om2(ij),KD_dyn,ZD_dyn)
|
||||
|
||||
Z2_dyn(kl) = dot_product(X2(:,ij),matmul(ZC_dyn,X2(:,ij))) &
|
||||
+ dot_product(Y2(:,ij),matmul(ZD_dyn,Y2(:,ij)))
|
||||
|
||||
Om2_dyn(kl) = dot_product(X2(:,ij),matmul(KC_dyn - KC_sta,X2(:,ij))) &
|
||||
- dot_product(Y2(:,ij),matmul(KD_dyn - KD_sta,Y2(:,ij))) &
|
||||
+ dot_product(X2(:,ij),matmul(KB_dyn - KB_sta,Y2(:,ij))) &
|
||||
- dot_product(Y2(:,ij),matmul(transpose(KB_dyn - KB_sta),X2(:,ij)))
|
||||
|
||||
end if
|
||||
|
||||
Z2_dyn(kl) = 1d0/(1d0 - Z2_dyn(kl))
|
||||
Om2_dyn(kl) = Z2_dyn(kl)*Om2_dyn(kl)
|
||||
|
||||
write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') &
|
||||
kl,Om2(ij)*HaToeV,(Om2(ij)+Om2_dyn(kl))*HaToeV,Om2_dyn(kl)*HaToeV,Z2_dyn(kl)
|
||||
|
||||
end do
|
||||
write(*,*) '---------------------------------------------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
end subroutine
|
||||
Loading…
Reference in New Issue
Block a user