PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-03 20:53:53 +01:00
Code Releases Activity

ppBSE dynamic kernel B and D

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-07-24 21:46:03 +02:00
parent c0da089370
commit b53f18fa79
1 changed files with 140 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

140
src/GW/GW_ppBSE_dynamic_kernel_D.f90 Normal file
View File

@ -0,0 +1,140 @@
subroutine GW_ppBSE_dynamic_kernel_D(ispin,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) :: ispin
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 :: 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 A matrix
if(ispin == 1) then
ij = 0
do i=nC+1,nO
do j=i,nO
ij = ij + 1
kl = 0
do k=nC+1,nO
do l=k,nO
kl = kl + 1
do m=1,nS
dem = - OmBSE + eGW(k) - 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(k) - 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
dem = - OmBSE + eGW(l) - 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(l) - 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
end do
end do
end do
end do
end do
end if
if(ispin == 2) then
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(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(k) - 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
dem = - OmBSE + eGW(l) - 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(l) - 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
end do
end do
end do
end do
end do
end if
end subroutine