diff --git a/src/GW/UGW_phBSE.f90 b/src/GW/UGW_phBSE.f90 index eeba00d..82ddc30 100644 --- a/src/GW/UGW_phBSE.f90 +++ b/src/GW/UGW_phBSE.f90 @@ -109,9 +109,9 @@ subroutine UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta, & !------------------------------------------------- if(dBSE) & - call unrestricted_Bethe_Salpeter_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,nS_sc, & - eW,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, & - OmRPA,rho_RPA,OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) + call UGW_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,nS_sc, & + eW,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, & + OmRPA,rho_RPA,OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) deallocate(OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) @@ -145,9 +145,9 @@ subroutine UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta, & !------------------------------------------------- if(dBSE) & - call unrestricted_Bethe_Salpeter_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,nS_sc, & - eW,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, & - OmRPA,rho_RPA,OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf) + call UGW_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,nS_sc, & + eW,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, & + OmRPA,rho_RPA,OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf) deallocate(OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf) diff --git a/src/GW/unrestricted_Bethe_Salpeter_A_matrix_dynamic.f90 b/src/GW/UGW_phBSE_dynamic_kernel_A.f90 similarity index 71% rename from src/GW/unrestricted_Bethe_Salpeter_A_matrix_dynamic.f90 rename to src/GW/UGW_phBSE_dynamic_kernel_A.f90 index 3418fe7..06bc3de 100644 --- a/src/GW/unrestricted_Bethe_Salpeter_A_matrix_dynamic.f90 +++ b/src/GW/UGW_phBSE_dynamic_kernel_A.f90 @@ -1,5 +1,5 @@ -subroutine unrestricted_Bethe_Salpeter_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,eGW, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,OmBSE,A_dyn) +subroutine UGW_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,eGW, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,OmBSE,A_dyn,ZA_dyn) ! Compute the extra term for dynamical Bethe-Salpeter equation for linear response in the unrestricted formalism @@ -37,6 +37,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV, ! Output variables double precision,intent(out) :: A_dyn(nSt,nSt) + double precision,intent(out) :: ZA_dyn(nSt,nSt) !--------------------------------------------------! ! Build BSE matrix for spin-conserving transitions ! @@ -76,6 +77,19 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV, enddo A_dyn(ia,jb) = A_dyn(ia,jb) - lambda*chi + + chi = 0d0 + do kc=1,nS_sc + + eps = + OmBSE - OmRPA(kc) - (eGW(a,1) - eGW(j,1)) + chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,1)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + eps = + OmBSE - OmRPA(kc) - (eGW(b,1) - eGW(i,1)) + chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,1)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + enddo + + ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + lambda*chi enddo enddo @@ -112,6 +126,19 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV, enddo A_dyn(nSa+ia,nSa+jb) = A_dyn(nSa+ia,nSa+jb) - lambda*chi + + chi = 0d0 + do kc=1,nS_sc + + eps = + OmBSE - OmRPA(kc) - (eGW(a,2) - eGW(j,2)) + chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,2)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + eps = + OmBSE - OmRPA(kc) - (eGW(b,2) - eGW(i,2)) + chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,2)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + enddo + + ZA_dyn(nSa+ia,nSa+jb) = ZA_dyn(nSa+ia,nSa+jb) + lambda*chi enddo enddo @@ -157,6 +184,19 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV, A_dyn(ia,jb) = A_dyn(ia,jb) - lambda*chi + chi = 0d0 + do kc=1,nS_sc + + eps = + OmBSE - OmRPA(kc) - (eGW(a,2) - eGW(j,1)) + chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,2)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + eps = + OmBSE - OmRPA(kc) - (eGW(b,2) - eGW(i,1)) + chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,2)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + enddo + + ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + lambda*chi + end do end do end do @@ -193,6 +233,19 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV, A_dyn(nSa+ia,nSa+jb) = A_dyn(nSa+ia,nSa+jb) - lambda*chi + chi = 0d0 + do kc=1,nS_sc + + eps = + OmBSE - OmRPA(kc) - (eGW(a,1) - eGW(j,2)) + chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,1)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + eps = + OmBSE - OmRPA(kc) - (eGW(b,1) - eGW(i,2)) + chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,1)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + enddo + + ZA_dyn(nSa+ia,nSa+jb) = ZA_dyn(nSa+ia,nSa+jb) + lambda*chi + end do end do end do @@ -200,4 +253,4 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV, end if -end subroutine unrestricted_Bethe_Salpeter_A_matrix_dynamic +end subroutine diff --git a/src/GW/unrestricted_Bethe_Salpeter_dynamic_perturbation.f90 b/src/GW/UGW_phBSE_dynamic_perturbation.f90 similarity index 77% rename from src/GW/unrestricted_Bethe_Salpeter_dynamic_perturbation.f90 rename to src/GW/UGW_phBSE_dynamic_perturbation.f90 index 56a8ed9..75655c9 100644 --- a/src/GW/unrestricted_Bethe_Salpeter_dynamic_perturbation.f90 +++ b/src/GW/UGW_phBSE_dynamic_perturbation.f90 @@ -1,6 +1,6 @@ -subroutine unrestricted_Bethe_Salpeter_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,nS_sc,eW,eGW, & - ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, & - OmRPA,rho_RPA,OmBSE,XpY_BSE,XmY_BSE) +subroutine UGW_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,nS_sc,eW,eGW, & + ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, & + OmRPA,rho_RPA,OmBSE,XpY_BSE,XmY_BSE) ! Compute dynamical effects via perturbation theory for BSE @@ -82,15 +82,10 @@ subroutine unrestricted_Bethe_Salpeter_dynamic_perturbation(ispin,dTDA,eta,nBas, if(dTDA) then - ! Resonant part of the BSE correction for dynamical TDA + ! Resonant part of the BSE correction for dynamical TDA and its renormalization factor - call unrestricted_Bethe_Salpeter_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,1d0,eGW, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,OmBSE(ia),A_dyn) - - ! Renormalization factor of the resonant parts for dynamical TDA - - call unrestricted_Bethe_Salpeter_ZA_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,1d0,eGW, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,OmBSE(ia),ZA_dyn) + call UGW_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,1d0,eGW, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,OmBSE(ia),A_dyn,ZA_dyn) ZDyn(ia) = dot_product(X,matmul(ZA_dyn,X)) OmDyn(ia) = dot_product(X,matmul( A_dyn,X)) @@ -107,4 +102,4 @@ subroutine unrestricted_Bethe_Salpeter_dynamic_perturbation(ispin,dTDA,eta,nBas, write(*,*) '---------------------------------------------------------------------------------------------------' write(*,*) -end subroutine unrestricted_Bethe_Salpeter_dynamic_perturbation +end subroutine diff --git a/src/GW/unrestricted_Bethe_Salpeter_ZA_matrix_dynamic.f90 b/src/GW/unrestricted_Bethe_Salpeter_ZA_matrix_dynamic.f90 deleted file mode 100644 index 6452cc1..0000000 --- a/src/GW/unrestricted_Bethe_Salpeter_ZA_matrix_dynamic.f90 +++ /dev/null @@ -1,175 +0,0 @@ -subroutine unrestricted_Bethe_Salpeter_ZA_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,eGW, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,OmBSE,ZA_dyn) - -! Compute the extra term for dynamical Bethe-Salpeter equation for linear response in the unrestricted formalism - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: ispin - integer,intent(in) :: nBas - integer,intent(in) :: nC(nspin) - integer,intent(in) :: nO(nspin) - integer,intent(in) :: nV(nspin) - integer,intent(in) :: nR(nspin) - integer,intent(in) :: nSa - integer,intent(in) :: nSb - integer,intent(in) :: nSt - integer,intent(in) :: nS_sc - double precision,intent(in) :: eta - double precision,intent(in) :: lambda - double precision,intent(in) :: eGW(nBas,nspin) - double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) - double precision,intent(in) :: OmRPA(nS_sc) - double precision,intent(in) :: rho_RPA(nBas,nBas,nS_sc,nspin) - double precision,intent(in) :: OmBSE - -! Local variables - - double precision :: chi - double precision :: eps - integer :: i,j,a,b,ia,jb,kc - -! Output variables - - double precision,intent(out) :: ZA_dyn(nSt,nSt) - -!--------------------------------------------------! -! Build BSE matrix for spin-conserving transitions ! -!--------------------------------------------------! - - ZA_dyn(:,:) = 0d0 - - if(ispin == 1) then - - ! aaaa block - - ia = 0 - do i=nC(1)+1,nO(1) - do a=nO(1)+1,nBas-nR(1) - ia = ia + 1 - jb = 0 - do j=nC(1)+1,nO(1) - do b=nO(1)+1,nBas-nR(1) - jb = jb + 1 - - chi = 0d0 - do kc=1,nS_sc - - eps = + OmBSE - OmRPA(kc) - (eGW(a,1) - eGW(j,1)) - chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,1)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - - eps = + OmBSE - OmRPA(kc) - (eGW(b,1) - eGW(i,1)) - chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,1)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - - enddo - - ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + lambda*chi - - enddo - enddo - enddo - enddo - - ! bbbb block - - ia = 0 - do i=nC(2)+1,nO(2) - do a=nO(2)+1,nBas-nR(2) - ia = ia + 1 - jb = 0 - do j=nC(2)+1,nO(2) - do b=nO(2)+1,nBas-nR(2) - jb = jb + 1 - - chi = 0d0 - do kc=1,nS_sc - - eps = + OmBSE - OmRPA(kc) - (eGW(a,2) - eGW(j,2)) - chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,2)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - - eps = + OmBSE - OmRPA(kc) - (eGW(b,2) - eGW(i,2)) - chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,2)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - - enddo - - ZA_dyn(nSa+ia,nSa+jb) = ZA_dyn(nSa+ia,nSa+jb) + lambda*chi - - enddo - enddo - enddo - enddo - - end if - -!--------------------------------------------! -! Build BSE matrix for spin-flip transitions ! -!--------------------------------------------! - - if(ispin == 2) then - - ! abab block - - ia = 0 - do i=nC(1)+1,nO(1) - do a=nO(2)+1,nBas-nR(2) - ia = ia + 1 - jb = 0 - do j=nC(1)+1,nO(1) - do b=nO(2)+1,nBas-nR(2) - jb = jb + 1 - - chi = 0d0 - do kc=1,nS_sc - - eps = + OmBSE - OmRPA(kc) - (eGW(a,2) - eGW(j,1)) - chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,2)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - - eps = + OmBSE - OmRPA(kc) - (eGW(b,2) - eGW(i,1)) - chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,2)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - - enddo - - ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + lambda*chi - - end do - end do - end do - end do - - ! baba block - - ia = 0 - do i=nC(2)+1,nO(2) - do a=nO(1)+1,nBas-nR(1) - ia = ia + 1 - jb = 0 - do j=nC(2)+1,nO(2) - do b=nO(1)+1,nBas-nR(1) - jb = jb + 1 - - chi = 0d0 - do kc=1,nS_sc - - eps = + OmBSE - OmRPA(kc) - (eGW(a,1) - eGW(j,2)) - chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,1)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - - eps = + OmBSE - OmRPA(kc) - (eGW(b,1) - eGW(i,2)) - chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,1)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - - enddo - - ZA_dyn(nSa+ia,nSa+jb) = ZA_dyn(nSa+ia,nSa+jb) + lambda*chi - - end do - end do - end do - end do - - end if - -end subroutine unrestricted_Bethe_Salpeter_ZA_matrix_dynamic