PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-09-27 12:00:50 +02:00
Code Releases Activity

fix bug in ppBSE dynamic kernel and subroutine names

Browse Source
This commit is contained in:
Pierre-Francois Loos 2024-09-12 15:46:11 +02:00
parent 90cc34c81a
commit ff7cff0963
6 changed files with 30 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
src/GW/RGW_Im_SigC.f90
View File

@ -1,4 +1,4 @@
double precision function RGW_ImSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) double precision function RGW_Im_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
! Compute diagonal of the correlation part of the self-energy ! Compute diagonal of the correlation part of the self-energy
@ -27,7 +27,7 @@ double precision function RGW_ImSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
! Initialize ! Initialize
RGW_ImSigC = 0d0 RGW_Im_SigC = 0d0
! Occupied part of the correlation self-energy ! Occupied part of the correlation self-energy
@ -35,7 +35,7 @@ double precision function RGW_ImSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
do m=1,nS do m=1,nS
eps = w - e(i) + Om(m) eps = w - e(i) + Om(m)
num = 2d0*rho(p,i,m)**2 num = 2d0*rho(p,i,m)**2
RGW_ImSigC = RGW_ImSigC + num*eta/(eps**2 + eta**2) RGW_Im_SigC = RGW_Im_SigC + num*eta/(eps**2 + eta**2)
end do end do
end do end do
@ -45,7 +45,7 @@ double precision function RGW_ImSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
do m=1,nS do m=1,nS
eps = w - e(a) - Om(m) eps = w - e(a) - Om(m)
num = 2d0*rho(p,a,m)**2 num = 2d0*rho(p,a,m)**2
RGW_ImSigC = RGW_ImSigC - num*eta/(eps**2 + eta**2) RGW_Im_SigC = RGW_Im_SigC - num*eta/(eps**2 + eta**2)
end do end do
end do end do

8
src/GW/RGW_Im_dSigC.f90
View File

@ -1,4 +1,4 @@
double precision function RGW_ImdSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) double precision function RGW_Im_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
! Compute the derivative of the correlation part of the self-energy ! Compute the derivative of the correlation part of the self-energy
@ -27,7 +27,7 @@ double precision function RGW_ImdSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
! Initialize ! Initialize
RGW_ImdSigC = 0d0 RGW_Im_dSigC = 0d0
! Occupied part of the correlation self-energy ! Occupied part of the correlation self-energy
@ -35,7 +35,7 @@ double precision function RGW_ImdSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
do m=1,nS do m=1,nS
eps = w - e(i) + Om(m) eps = w - e(i) + Om(m)
num = 2d0*rho(p,i,m)**2 num = 2d0*rho(p,i,m)**2
RGW_ImdSigC = RGW_ImdSigC - 2d0*num*eps*eta/(eps**2 + eta**2)**2 RGW_Im_dSigC = RGW_Im_dSigC - 2d0*num*eps*eta/(eps**2 + eta**2)**2
end do end do
end do end do
@ -45,7 +45,7 @@ double precision function RGW_ImdSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
do m=1,nS do m=1,nS
eps = w - e(a) - Om(m) eps = w - e(a) - Om(m)
num = 2d0*rho(p,a,m)**2 num = 2d0*rho(p,a,m)**2
RGW_ImdSigC = RGW_ImdSigC + 2d0*num*eps*eta/(eps**2 + eta**2)**2 RGW_Im_dSigC = RGW_Im_dSigC + 2d0*num*eps*eta/(eps**2 + eta**2)**2
end do end do
end do end do

8
src/GW/RGW_plot_self_energy.f90
View File

@ -24,7 +24,7 @@ subroutine RGW_plot_self_energy(nBas,eta,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
integer :: p,g integer :: p,g
integer :: nGrid integer :: nGrid
double precision :: wmin,wmax,dw double precision :: wmin,wmax,dw
double precision,external :: RGW_ReSigC,RGW_ImSigC,RGW_RedSigC double precision,external :: RGW_Re_SigC,RGW_Im_SigC,RGW_Re_dSigC
double precision,allocatable :: w(:) double precision,allocatable :: w(:)
double precision,allocatable :: ReSigC(:,:),ImSigC(:,:) double precision,allocatable :: ReSigC(:,:),ImSigC(:,:)
double precision,allocatable :: Z(:,:) double precision,allocatable :: Z(:,:)
@ -56,9 +56,9 @@ subroutine RGW_plot_self_energy(nBas,eta,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
do g=1,nGrid do g=1,nGrid
do p=nC+1,nBas-nR do p=nC+1,nBas-nR
ReSigC(p,g) = RGW_ReSigC(p,w(g),eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho) ReSigC(p,g) = RGW_Re_SigC(p,w(g),eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho)
ImSigC(p,g) = RGW_ImSigC(p,w(g),eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho) ImSigC(p,g) = RGW_Im_SigC(p,w(g),eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho)
Z(p,g) = RGW_RedSigC(p,w(g),eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho) Z(p,g) = RGW_Re_dSigC(p,w(g),eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho)
end do end do
end do end do

14
src/GW/RGW_ppBSE.f90
View File

@ -176,7 +176,7 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
!call wall_time(tt2) !call wall_time(tt2)
!write(*,'(A,1X,F10.3)'), 'wall time for ppLR (sec)', tt2-tt1 !write(*,'(A,1X,F10.3)'), 'wall time for ppLR (sec)', tt2-tt1
deallocate(Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta) deallocate(Bpp,Cpp,Dpp)
! !
! print*, 'LAPACK:' ! print*, 'LAPACK:'
! print*, Om2 ! print*, Om2
@ -241,10 +241,11 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
if(dBSE) & if(dBSE) &
call RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, & call RGW_ppBSE_dynamic_perturbation(ispin,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) Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta)
deallocate(KB_sta,KC_sta,KD_sta)
deallocate(Om1,X1,Y1,Om2,X2,Y2) deallocate(Om1,X1,Y1,Om2,X2,Y2)
end if end if
!------------------- !-------------------
@ -289,7 +290,7 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
Dpp(:,:) = Dpp(:,:) + KD_sta(:,:) Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin)) call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
deallocate(Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta) deallocate(Bpp,Cpp,Dpp)
!print*, 'LAPACK:' !print*, 'LAPACK:'
!print*, Om2 !print*, Om2
@ -351,9 +352,10 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
!----------------------------------------------------! !----------------------------------------------------!
if(dBSE) & if(dBSE) &
call RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, & call RGW_ppBSE_dynamic_perturbation(ispin,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) Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta)
deallocate(KB_sta,KC_sta,KD_sta)
deallocate(Om1,X1,Y1,Om2,X2,Y2) deallocate(Om1,X1,Y1,Om2,X2,Y2)
end if end if

16
src/GW/RGW_ppBSE_dynamic_kernel_D.f90
View File

@ -56,25 +56,25 @@ subroutine RGW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,e
do m=1,nS do m=1,nS
dem = - OmBSE + eGW(k) - Om(m) + eGW(j) dem = - OmBSE + eGW(k) - Om(m) + eGW(i)
num = rho(i,k,m)*rho(j,l,m) num = rho(i,k,m)*rho(j,l,m)
KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2) 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 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) dem = - OmBSE + eGW(k) - Om(m) + eGW(j)
num = rho(j,k,m)*rho(i,l,m) num = rho(j,k,m)*rho(i,l,m)
KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2) 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 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) dem = - OmBSE + eGW(l) - Om(m) + eGW(j)
num = rho(i,k,m)*rho(j,l,m) num = rho(i,k,m)*rho(j,l,m)
KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2) 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 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) dem = - OmBSE + eGW(l) - Om(m) + eGW(i)
num = rho(j,k,m)*rho(i,l,m) num = rho(j,k,m)*rho(i,l,m)
KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2) KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2)
@ -106,25 +106,25 @@ subroutine RGW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,e
kl = kl + 1 kl = kl + 1
do m=1,nS do m=1,nS
dem = - OmBSE + eGW(k) - Om(m) + eGW(j) dem = - OmBSE + eGW(k) - Om(m) + eGW(i)
num = rho(i,k,m)*rho(j,l,m) num = rho(i,k,m)*rho(j,l,m)
KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2) 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 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) dem = - OmBSE + eGW(k) - Om(m) + eGW(j)
num = rho(j,k,m)*rho(i,l,m) num = rho(j,k,m)*rho(i,l,m)
KD_dyn(ij,kl) = KD_dyn(ij,kl) - num*dem/(dem**2 + eta**2) 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 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) dem = - OmBSE + eGW(l) - Om(m) + eGW(j)
num = rho(i,k,m)*rho(j,l,m) num = rho(i,k,m)*rho(j,l,m)
KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2) 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 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) dem = - OmBSE + eGW(l) - Om(m) + eGW(i)
num = rho(j,k,m)*rho(i,l,m) num = rho(j,k,m)*rho(i,l,m)
KD_dyn(ij,kl) = KD_dyn(ij,kl) - num*dem/(dem**2 + eta**2) KD_dyn(ij,kl) = KD_dyn(ij,kl) - num*dem/(dem**2 + eta**2)

4
src/GW/RGW_ppBSE_dynamic_perturbation.f90
View File

@ -1,5 +1,5 @@
subroutine RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int, & subroutine RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,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) OmRPA,rho_RPA,Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta)
! Compute dynamical effects via perturbation theory for BSE ! Compute dynamical effects via perturbation theory for BSE
@ -121,7 +121,7 @@ subroutine RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO
call RGW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,eGW,OmRPA,rho_RPA,Om2(ij),KD_dyn,ZD_dyn) call RGW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,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))) 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))) Om2_dyn(kl) = - dot_product(Y2(:,ij),matmul(KD_dyn - KD_sta,Y2(:,ij)))
else else