LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-06-21 20:52:21 +02:00
Code Releases Activity

problem minus sign

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-06-01 17:26:52 +02:00
parent 96207d9c58
commit 9a0c7bfac7
4 changed files with 73 additions and 112 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

60
src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90
View File

@ -1,4 +1,4 @@
subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,Ap,Am,Bp,Bm) subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,A_dyn,B_dyn)
! Compute the dynamic part of the Bethe-Salpeter equation matrices ! Compute the dynamic part of the Bethe-Salpeter equation matrices
@ -18,24 +18,18 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O
! Local variables ! Local variables
integer :: maxS integer :: maxS
double precision :: chi_A,chi_B,eps double precision :: chi_A,chi_B,eps,eps_A,eps_B
double precision :: chi_Ap,chi_Am,chi_Bp,chi_Bm
double precision :: eps_Ap,eps_Am,eps_Bp,eps_Bm
integer :: i,j,a,b,ia,jb,kc integer :: i,j,a,b,ia,jb,kc
! Output variables ! Output variables
double precision,intent(out) :: Ap(nS,nS) double precision,intent(out) :: A_dyn(nS,nS)
double precision,intent(out) :: Am(nS,nS) double precision,intent(out) :: B_dyn(nS,nS)
double precision,intent(out) :: Bp(nS,nS)
double precision,intent(out) :: Bm(nS,nS)
! Initialization ! Initialization
Ap(:,:) = 0d0 A_dyn(:,:) = 0d0
Am(:,:) = 0d0 B_dyn(:,:) = 0d0
Bp(:,:) = 0d0
Bm(:,:) = 0d0
! Number of poles taken into account ! Number of poles taken into account
@ -63,45 +57,31 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O
enddo enddo
Ap(ia,jb) = Ap(ia,jb) - 4d0*lambda*chi_A A_dyn(ia,jb) = A_dyn(ia,jb) - 4d0*lambda*chi_A
Am(ia,jb) = Am(ia,jb) - 4d0*lambda*chi_A B_dyn(ia,jb) = B_dyn(ia,jb) - 4d0*lambda*chi_B
Bp(ia,jb) = Bp(ia,jb) - 4d0*lambda*chi_B
Bm(ia,jb) = Bm(ia,jb) - 4d0*lambda*chi_B
chi_Ap = 0d0 chi_A = 0d0
chi_Am = 0d0 chi_B = 0d0
chi_Bp = 0d0
chi_Bm = 0d0
do kc=1,maxS do kc=1,maxS
eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 eps_A = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2
eps_Am = (- OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 chi_A = chi_A + rho(i,j,kc)*rho(a,b,kc)*(+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_A
chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Ap
chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(- OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Am
eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 eps_A = (+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2
eps_Am = (- OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 chi_A = chi_A + rho(i,j,kc)*rho(a,b,kc)*(+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_A
chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Ap
chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(- OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Am
eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 eps_B = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2
eps_Bm = (- OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 chi_B = chi_B + rho(i,b,kc)*rho(a,j,kc)*(+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_B
chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_Bp
chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(- OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_Bm
eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 eps_B = (+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2
eps_Bm = (- OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 chi_B = chi_B + rho(i,b,kc)*rho(a,j,kc)*(+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_B
chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bp
chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(- OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bm
enddo enddo
Ap(ia,jb) = Ap(ia,jb) - 2d0*lambda*chi_Ap A_dyn(ia,jb) = A_dyn(ia,jb) - 2d0*lambda*chi_A
Am(ia,jb) = Am(ia,jb) - 2d0*lambda*chi_Am
Bp(ia,jb) = Bp(ia,jb) - 2d0*lambda*chi_Bp B_dyn(ia,jb) = B_dyn(ia,jb) - 2d0*lambda*chi_B
Bm(ia,jb) = Bm(ia,jb) - 2d0*lambda*chi_Bm
enddo enddo
enddo enddo

54
src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90
View File

@ -1,4 +1,4 @@
subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,ZAp,ZAm,ZBp,ZBm) subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,ZA,ZB)
! Compute the dynamic part of the Bethe-Salpeter equation matrices ! Compute the dynamic part of the Bethe-Salpeter equation matrices
@ -18,23 +18,19 @@ subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,
! Local variables ! Local variables
integer :: maxS integer :: maxS
double precision :: chi_Ap,chi_Am,chi_Bp,chi_Bm double precision :: chi_A,chi_B
double precision :: eps_Ap,eps_Am,eps_Bp,eps_Bm double precision :: eps_A,eps_B
integer :: i,j,a,b,ia,jb,kc integer :: i,j,a,b,ia,jb,kc
! Output variables ! Output variables
double precision,intent(out) :: ZAp(nS,nS) double precision,intent(out) :: ZA(nS,nS)
double precision,intent(out) :: ZAm(nS,nS) double precision,intent(out) :: ZB(nS,nS)
double precision,intent(out) :: ZBp(nS,nS)
double precision,intent(out) :: ZBm(nS,nS)
! Initialization ! Initialization
ZAp(:,:) = 0d0 ZA(:,:) = 0d0
ZAm(:,:) = 0d0 ZB(:,:) = 0d0
ZBp(:,:) = 0d0
ZBm(:,:) = 0d0
! Number of poles taken into account ! Number of poles taken into account
@ -51,40 +47,28 @@ subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,
do b=nO+1,nBas-nR do b=nO+1,nBas-nR
jb = jb + 1 jb = jb + 1
chi_Ap = 0d0 chi_A = 0d0
chi_Am = 0d0 chi_B = 0d0
chi_Bp = 0d0
chi_Bm = 0d0
do kc=1,maxS do kc=1,maxS
eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 eps_A = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2
eps_Am = (- OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 chi_A = chi_A + rho(i,j,kc)*rho(a,b,kc)*((+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_A)**2
chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*((+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Ap)**2
chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*((- OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Am)**2
eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 eps_A = (+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2
eps_Am = (- OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 chi_A = chi_A + rho(i,j,kc)*rho(a,b,kc)*((+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_A)**2
chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*((+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Ap)**2
chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*((- OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Am)**2
eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 eps_B = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2
eps_Bm = (- OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 chi_B = chi_B + rho(i,b,kc)*rho(a,j,kc)*((+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_B)**2
chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*((+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_Bp)**2
chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*((- OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_Bm)**2
eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 eps_B = (+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2
eps_Bm = (- OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 chi_B = chi_B + rho(i,b,kc)*rho(a,j,kc)*((+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_B)**2
chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*((+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bp)**2
chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*((- OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bm)**2
enddo enddo
ZAp(ia,jb) = ZAp(ia,jb) + 2d0*lambda*chi_Ap ZA(ia,jb) = ZA(ia,jb) + 2d0*lambda*chi_A
ZAm(ia,jb) = ZAm(ia,jb) - 2d0*lambda*chi_Am
ZBp(ia,jb) = ZBp(ia,jb) + 2d0*lambda*chi_Bp ZB(ia,jb) = ZB(ia,jb) + 2d0*lambda*chi_B
ZBm(ia,jb) = ZBm(ia,jb) - 2d0*lambda*chi_Bm
enddo enddo
enddo enddo

46
src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90
View File

@ -35,21 +35,17 @@ subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW,O
double precision,allocatable :: X(:) double precision,allocatable :: X(:)
double precision,allocatable :: Y(:) double precision,allocatable :: Y(:)
double precision,allocatable :: Ap_dyn(:,:) double precision,allocatable :: A_dyn(:,:)
double precision,allocatable :: Am_dyn(:,:) double precision,allocatable :: ZA_dyn(:,:)
double precision,allocatable :: ZAp_dyn(:,:)
double precision,allocatable :: ZAm_dyn(:,:)
double precision,allocatable :: Bp_dyn(:,:) double precision,allocatable :: B_dyn(:,:)
double precision,allocatable :: Bm_dyn(:,:) double precision,allocatable :: ZB_dyn(:,:)
double precision,allocatable :: ZBp_dyn(:,:)
double precision,allocatable :: ZBm_dyn(:,:)
! Memory allocation ! Memory allocation
allocate(OmDyn(nS),ZDyn(nS),X(nS),Y(nS),Ap_dyn(nS,nS),ZAp_dyn(nS,nS)) allocate(OmDyn(nS),ZDyn(nS),X(nS),Y(nS),A_dyn(nS,nS),ZA_dyn(nS,nS))
if(.not.dTDA) allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),Bp_dyn(nS,nS),Bm_dyn(nS,nS),ZBp_dyn(nS,nS),ZBm_dyn(nS,nS)) if(.not.dTDA) allocate(B_dyn(nS,nS),ZB_dyn(nS,nS))
gapGW = eGW(nO+1) - eGW(nO) gapGW = eGW(nO+1) - eGW(nO)
@ -70,36 +66,38 @@ subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW,O
! Resonant part of the BSE correction for dynamical TDA ! Resonant part of the BSE correction for dynamical TDA
call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:),Ap_dyn(:,:)) call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:), &
A_dyn(:,:))
! Renormalization factor of the resonant parts for dynamical TDA ! Renormalization factor of the resonant parts for dynamical TDA
call Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:),ZAp_dyn(:,:)) call Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:), &
ZA_dyn(:,:))
ZDyn(ia) = dot_product(X(:),matmul(ZAp_dyn(:,:),X(:))) ZDyn(ia) = dot_product(X(:),matmul(ZA_dyn(:,:),X(:)))
OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) OmDyn(ia) = dot_product(X(:),matmul(A_dyn(:,:),X(:)))
else else
! Resonant and anti-resonant part of the BSE correction ! Resonant and anti-resonant part of the BSE correction
call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:), & call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:), &
Ap_dyn(:,:),Am_dyn(:,:),Bp_dyn(:,:),Bm_dyn(:,:)) A_dyn(:,:),B_dyn(:,:))
! Renormalization factor of the resonant and anti-resonant parts ! Renormalization factor of the resonant and anti-resonant parts
call Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:), & call Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:), &
ZAp_dyn(:,:),ZAm_dyn(:,:),ZBp_dyn(:,:),ZBm_dyn(:,:)) ZA_dyn(:,:),ZB_dyn(:,:))
ZDyn(ia) = dot_product(X(:),matmul(ZAp_dyn(:,:),X(:))) & ZDyn(ia) = dot_product(X(:),matmul(ZA_dyn(:,:),X(:))) &
- dot_product(Y(:),matmul(ZAm_dyn(:,:),Y(:))) & - dot_product(Y(:),matmul(ZA_dyn(:,:),Y(:))) &
+ dot_product(X(:),matmul(ZBp_dyn(:,:),Y(:))) & + dot_product(X(:),matmul(ZB_dyn(:,:),Y(:))) &
- dot_product(Y(:),matmul(ZBm_dyn(:,:),X(:))) - dot_product(Y(:),matmul(ZB_dyn(:,:),X(:)))
OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) & OmDyn(ia) = dot_product(X(:),matmul(A_dyn(:,:),X(:))) &
- dot_product(Y(:),matmul(Am_dyn(:,:),Y(:))) & - dot_product(Y(:),matmul(A_dyn(:,:),Y(:))) &
+ dot_product(X(:),matmul(Bp_dyn(:,:),Y(:))) & + dot_product(X(:),matmul(B_dyn(:,:),Y(:))) &
- dot_product(Y(:),matmul(Bm_dyn(:,:),X(:))) - dot_product(Y(:),matmul(B_dyn(:,:),X(:)))
end if end if

25
src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90
View File

@ -40,16 +40,14 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(TDA,eta,nBas,nC,nO,nV,n
double precision,allocatable :: OmOld(:) double precision,allocatable :: OmOld(:)
double precision,allocatable :: X(:) double precision,allocatable :: X(:)
double precision,allocatable :: Y(:) double precision,allocatable :: Y(:)
double precision,allocatable :: Ap_dyn(:,:) double precision,allocatable :: A_dyn(:,:)
double precision,allocatable :: Am_dyn(:,:) double precision,allocatable :: B_dyn(:,:)
double precision,allocatable :: Bp_dyn(:,:)
double precision,allocatable :: Bm_dyn(:,:)
! Memory allocation ! Memory allocation
allocate(OmDyn(nS),OmOld(nS),X(nS),Y(nS),Ap_dyn(nS,nS)) allocate(OmDyn(nS),OmOld(nS),X(nS),Y(nS),A_dyn(nS,nS))
if(.not.dTDA) allocate(Am_dyn(nS,nS),Bp_dyn(nS,nS),Bm_dyn(nS,nS)) if(.not.dTDA) allocate(B_dyn(nS,nS))
gapGW = eGW(nO+1) - eGW(nO) gapGW = eGW(nO+1) - eGW(nO)
@ -84,21 +82,22 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(TDA,eta,nBas,nC,nO,nV,n
! Resonant part of the BSE correction ! Resonant part of the BSE correction
call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmOld(ia),rho(:,:,:),Ap_dyn(:,:)) call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmOld(ia),rho(:,:,:), &
A_dyn(:,:))
OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) OmDyn(ia) = dot_product(X(:),matmul(A_dyn(:,:),X(:)))
else else
! Anti-resonant part of the BSE correction ! Anti-resonant part of the BSE correction
call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmOld(ia),rho(:,:,:), & call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmOld(ia),rho(:,:,:), &
Ap_dyn(:,:),Am_dyn(:,:),Bp_dyn(:,:),Bm_dyn(:,:)) A_dyn(:,:),B_dyn(:,:))
OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) & OmDyn(ia) = dot_product(X(:),matmul(A_dyn(:,:),X(:))) &
- dot_product(Y(:),matmul(Am_dyn(:,:),Y(:))) & - dot_product(Y(:),matmul(A_dyn(:,:),Y(:))) &
+ dot_product(X(:),matmul(Bp_dyn(:,:),Y(:))) & + dot_product(X(:),matmul(B_dyn(:,:),Y(:))) &
- dot_product(Y(:),matmul(Bm_dyn(:,:),X(:))) - dot_product(Y(:),matmul(B_dyn(:,:),X(:)))
end if end if