diff --git a/src/GW/RGW_Im_SigC.f90 b/src/GW/RGW_Im_SigC.f90 index 2763afa..97f6748 100644 --- a/src/GW/RGW_Im_SigC.f90 +++ b/src/GW/RGW_Im_SigC.f90 @@ -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 @@ -27,7 +27,7 @@ double precision function RGW_ImSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Initialize - RGW_ImSigC = 0d0 + RGW_Im_SigC = 0d0 ! 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 eps = w - e(i) + Om(m) 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 @@ -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 eps = w - e(a) - Om(m) 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 diff --git a/src/GW/RGW_Im_dSigC.f90 b/src/GW/RGW_Im_dSigC.f90 index e32d26c..22a050e 100644 --- a/src/GW/RGW_Im_dSigC.f90 +++ b/src/GW/RGW_Im_dSigC.f90 @@ -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 @@ -27,7 +27,7 @@ double precision function RGW_ImdSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho) ! Initialize - RGW_ImdSigC = 0d0 + RGW_Im_dSigC = 0d0 ! 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 eps = w - e(i) + Om(m) 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 @@ -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 eps = w - e(a) - Om(m) 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 diff --git a/src/GW/RGW_plot_self_energy.f90 b/src/GW/RGW_plot_self_energy.f90 index f7bdf06..5d34b80 100644 --- a/src/GW/RGW_plot_self_energy.f90 +++ b/src/GW/RGW_plot_self_energy.f90 @@ -24,7 +24,7 @@ subroutine RGW_plot_self_energy(nBas,eta,nC,nO,nV,nR,nS,eHF,eGW,Om,rho) integer :: p,g integer :: nGrid 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 :: ReSigC(:,:),ImSigC(:,:) 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 p=nC+1,nBas-nR - ReSigC(p,g) = RGW_ReSigC(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) - Z(p,g) = RGW_RedSigC(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_Im_SigC(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 diff --git a/src/GW/RGW_ppBSE.f90 b/src/GW/RGW_ppBSE.f90 index 48bc225..a823c2f 100644 --- a/src/GW/RGW_ppBSE.f90 +++ b/src/GW/RGW_ppBSE.f90 @@ -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) !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*, Om2 @@ -241,10 +241,11 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS 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, & - 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) + 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(:,:) 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*, Om2 @@ -351,9 +352,10 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS !----------------------------------------------------! 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, & - Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta) + 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) + deallocate(KB_sta,KC_sta,KD_sta) deallocate(Om1,X1,Y1,Om2,X2,Y2) end if diff --git a/src/GW/RGW_ppBSE_dynamic_kernel_D.f90 b/src/GW/RGW_ppBSE_dynamic_kernel_D.f90 index e56fa80..d9d61c1 100644 --- a/src/GW/RGW_ppBSE_dynamic_kernel_D.f90 +++ b/src/GW/RGW_ppBSE_dynamic_kernel_D.f90 @@ -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 - 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) 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) + 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(i) + 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(j) + 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) @@ -106,25 +106,25 @@ subroutine RGW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,e kl = kl + 1 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) 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) + 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(i) + 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(j) + 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) diff --git a/src/GW/RGW_ppBSE_dynamic_perturbation.f90 b/src/GW/RGW_ppBSE_dynamic_perturbation.f90 index 058a68a..d7422a5 100644 --- a/src/GW/RGW_ppBSE_dynamic_perturbation.f90 +++ b/src/GW/RGW_ppBSE_dynamic_perturbation.f90 @@ -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, & - 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 @@ -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) - 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))) else