From 9a0c7bfac7033c9c481fbb8f5e84b1e81d2280ea Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 1 Jun 2020 17:26:52 +0200 Subject: [PATCH] problem minus sign --- .../Bethe_Salpeter_AB_matrix_dynamic.f90 | 60 +++++++------------ .../Bethe_Salpeter_ZAB_matrix_dynamic.f90 | 54 ++++++----------- .../Bethe_Salpeter_dynamic_perturbation.f90 | 46 +++++++------- ...alpeter_dynamic_perturbation_iterative.f90 | 25 ++++---- 4 files changed, 73 insertions(+), 112 deletions(-) diff --git a/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 index 20c1bb4..a3bcc02 100644 --- a/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 +++ b/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 @@ -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 @@ -18,24 +18,18 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O ! Local variables integer :: maxS - double precision :: chi_A,chi_B,eps - double precision :: chi_Ap,chi_Am,chi_Bp,chi_Bm - double precision :: eps_Ap,eps_Am,eps_Bp,eps_Bm + double precision :: chi_A,chi_B,eps,eps_A,eps_B integer :: i,j,a,b,ia,jb,kc ! Output variables - double precision,intent(out) :: Ap(nS,nS) - double precision,intent(out) :: Am(nS,nS) - double precision,intent(out) :: Bp(nS,nS) - double precision,intent(out) :: Bm(nS,nS) + double precision,intent(out) :: A_dyn(nS,nS) + double precision,intent(out) :: B_dyn(nS,nS) ! Initialization - Ap(:,:) = 0d0 - Am(:,:) = 0d0 - Bp(:,:) = 0d0 - Bm(:,:) = 0d0 + A_dyn(:,:) = 0d0 + B_dyn(:,:) = 0d0 ! 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 - Ap(ia,jb) = Ap(ia,jb) - 4d0*lambda*chi_A - Am(ia,jb) = Am(ia,jb) - 4d0*lambda*chi_A - Bp(ia,jb) = Bp(ia,jb) - 4d0*lambda*chi_B - Bm(ia,jb) = Bm(ia,jb) - 4d0*lambda*chi_B + A_dyn(ia,jb) = A_dyn(ia,jb) - 4d0*lambda*chi_A + B_dyn(ia,jb) = B_dyn(ia,jb) - 4d0*lambda*chi_B - chi_Ap = 0d0 - chi_Am = 0d0 - chi_Bp = 0d0 - chi_Bm = 0d0 + chi_A = 0d0 + chi_B = 0d0 do kc=1,maxS - eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 - eps_Am = (- OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 - 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_A = (+ 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 - eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 - eps_Am = (- OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 - 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_A = (+ 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 - eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 - eps_Bm = (- OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 - 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_B = (+ 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 - eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 - eps_Bm = (- OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 - 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 + eps_B = (+ 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 enddo - Ap(ia,jb) = Ap(ia,jb) - 2d0*lambda*chi_Ap - Am(ia,jb) = Am(ia,jb) - 2d0*lambda*chi_Am + A_dyn(ia,jb) = A_dyn(ia,jb) - 2d0*lambda*chi_A - Bp(ia,jb) = Bp(ia,jb) - 2d0*lambda*chi_Bp - Bm(ia,jb) = Bm(ia,jb) - 2d0*lambda*chi_Bm + B_dyn(ia,jb) = B_dyn(ia,jb) - 2d0*lambda*chi_B enddo enddo diff --git a/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 index f89999d..4bacbcc 100644 --- a/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 +++ b/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 @@ -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 @@ -18,23 +18,19 @@ subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW, ! Local variables integer :: maxS - double precision :: chi_Ap,chi_Am,chi_Bp,chi_Bm - double precision :: eps_Ap,eps_Am,eps_Bp,eps_Bm + double precision :: chi_A,chi_B + double precision :: eps_A,eps_B integer :: i,j,a,b,ia,jb,kc ! Output variables - double precision,intent(out) :: ZAp(nS,nS) - double precision,intent(out) :: ZAm(nS,nS) - double precision,intent(out) :: ZBp(nS,nS) - double precision,intent(out) :: ZBm(nS,nS) + double precision,intent(out) :: ZA(nS,nS) + double precision,intent(out) :: ZB(nS,nS) ! Initialization - ZAp(:,:) = 0d0 - ZAm(:,:) = 0d0 - ZBp(:,:) = 0d0 - ZBm(:,:) = 0d0 + ZA(:,:) = 0d0 + ZB(:,:) = 0d0 ! 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 jb = jb + 1 - chi_Ap = 0d0 - chi_Am = 0d0 - chi_Bp = 0d0 - chi_Bm = 0d0 + chi_A = 0d0 + chi_B = 0d0 do kc=1,maxS - eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 - eps_Am = (- OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**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_A = (+ 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 - eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 - eps_Am = (- OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**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_A = (+ 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 - eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 - eps_Bm = (- OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**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_B = (+ 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 - eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 - eps_Bm = (- OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**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 + eps_B = (+ 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 enddo - ZAp(ia,jb) = ZAp(ia,jb) + 2d0*lambda*chi_Ap - ZAm(ia,jb) = ZAm(ia,jb) - 2d0*lambda*chi_Am + ZA(ia,jb) = ZA(ia,jb) + 2d0*lambda*chi_A - ZBp(ia,jb) = ZBp(ia,jb) + 2d0*lambda*chi_Bp - ZBm(ia,jb) = ZBm(ia,jb) - 2d0*lambda*chi_Bm + ZB(ia,jb) = ZB(ia,jb) + 2d0*lambda*chi_B enddo enddo diff --git a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 index 3e41fae..cf480a9 100644 --- a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 +++ b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 @@ -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 :: Y(:) - double precision,allocatable :: Ap_dyn(:,:) - double precision,allocatable :: Am_dyn(:,:) - double precision,allocatable :: ZAp_dyn(:,:) - double precision,allocatable :: ZAm_dyn(:,:) + double precision,allocatable :: A_dyn(:,:) + double precision,allocatable :: ZA_dyn(:,:) - double precision,allocatable :: Bp_dyn(:,:) - double precision,allocatable :: Bm_dyn(:,:) - double precision,allocatable :: ZBp_dyn(:,:) - double precision,allocatable :: ZBm_dyn(:,:) + double precision,allocatable :: B_dyn(:,:) + double precision,allocatable :: ZB_dyn(:,:) ! 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) @@ -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 - 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 - 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(:))) - OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) + ZDyn(ia) = dot_product(X(:),matmul(ZA_dyn(:,:),X(:))) + OmDyn(ia) = dot_product(X(:),matmul(A_dyn(:,:),X(:))) else ! 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(:,:,:), & - Ap_dyn(:,:),Am_dyn(:,:),Bp_dyn(:,:),Bm_dyn(:,:)) + A_dyn(:,:),B_dyn(:,:)) ! 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(:,:,:), & - ZAp_dyn(:,:),ZAm_dyn(:,:),ZBp_dyn(:,:),ZBm_dyn(:,:)) + ZA_dyn(:,:),ZB_dyn(:,:)) - ZDyn(ia) = dot_product(X(:),matmul(ZAp_dyn(:,:),X(:))) & - - dot_product(Y(:),matmul(ZAm_dyn(:,:),Y(:))) & - + dot_product(X(:),matmul(ZBp_dyn(:,:),Y(:))) & - - dot_product(Y(:),matmul(ZBm_dyn(:,:),X(:))) + ZDyn(ia) = dot_product(X(:),matmul(ZA_dyn(:,:),X(:))) & + - dot_product(Y(:),matmul(ZA_dyn(:,:),Y(:))) & + + dot_product(X(:),matmul(ZB_dyn(:,:),Y(:))) & + - dot_product(Y(:),matmul(ZB_dyn(:,:),X(:))) - OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) & - - dot_product(Y(:),matmul(Am_dyn(:,:),Y(:))) & - + dot_product(X(:),matmul(Bp_dyn(:,:),Y(:))) & - - dot_product(Y(:),matmul(Bm_dyn(:,:),X(:))) + OmDyn(ia) = dot_product(X(:),matmul(A_dyn(:,:),X(:))) & + - dot_product(Y(:),matmul(A_dyn(:,:),Y(:))) & + + dot_product(X(:),matmul(B_dyn(:,:),Y(:))) & + - dot_product(Y(:),matmul(B_dyn(:,:),X(:))) end if diff --git a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 index 2cd4eea..09926dd 100644 --- a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 +++ b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 @@ -40,16 +40,14 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(TDA,eta,nBas,nC,nO,nV,n double precision,allocatable :: OmOld(:) double precision,allocatable :: X(:) double precision,allocatable :: Y(:) - double precision,allocatable :: Ap_dyn(:,:) - double precision,allocatable :: Am_dyn(:,:) - double precision,allocatable :: Bp_dyn(:,:) - double precision,allocatable :: Bm_dyn(:,:) + double precision,allocatable :: A_dyn(:,:) + double precision,allocatable :: B_dyn(:,:) ! 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) @@ -84,21 +82,22 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(TDA,eta,nBas,nC,nO,nV,n ! 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 ! 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(:,:,:), & - Ap_dyn(:,:),Am_dyn(:,:),Bp_dyn(:,:),Bm_dyn(:,:)) + A_dyn(:,:),B_dyn(:,:)) - OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) & - - dot_product(Y(:),matmul(Am_dyn(:,:),Y(:))) & - + dot_product(X(:),matmul(Bp_dyn(:,:),Y(:))) & - - dot_product(Y(:),matmul(Bm_dyn(:,:),X(:))) + OmDyn(ia) = dot_product(X(:),matmul(A_dyn(:,:),X(:))) & + - dot_product(Y(:),matmul(A_dyn(:,:),Y(:))) & + + dot_product(X(:),matmul(B_dyn(:,:),Y(:))) & + - dot_product(Y(:),matmul(B_dyn(:,:),X(:))) end if