diff --git a/src/GT/RG0T0pp.f90 b/src/GT/RG0T0pp.f90 index 1bd6db4..8569885 100644 --- a/src/GT/RG0T0pp.f90 +++ b/src/GT/RG0T0pp.f90 @@ -256,10 +256,10 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp correlation energy (singlet) =',EcBSE(1),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp correlation energy (triplet) =',EcBSE(2),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp correlation energy =',sum(EcBSE),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp total energy =',ENuc + ERHF + sum(EcBSE),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp@RHF correlation energy (singlet) = ',EcBSE(1),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp@RHF correlation energy (triplet) = ',EcBSE(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp@RHF correlation energy = ',sum(EcBSE),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/GT/RGTpp_phBSE.f90 b/src/GT/RGTpp_phBSE.f90 index b225406..21a5126 100644 --- a/src/GT/RGTpp_phBSE.f90 +++ b/src/GT/RGTpp_phBSE.f90 @@ -1,6 +1,6 @@ -subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & - Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,rho1ab,rho2ab,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, & - ERI,dipole_int,eT,eGT,EcBSE) +subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & + Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,rho1ab,rho2ab,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, & + ERI,dipole_int,eT,eGT,EcBSE) ! Compute the Bethe-Salpeter excitation energies with the T-matrix kernel @@ -18,7 +18,7 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n logical,intent(in) :: triplet double precision,intent(in) :: eta - integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -30,10 +30,10 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n integer,intent(in) :: nVVab integer,intent(in) :: nVVaa - double precision,intent(in) :: eT(nBas) - double precision,intent(in) :: eGT(nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: eT(nOrb) + double precision,intent(in) :: eGT(nOrb) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart) double precision,intent(in) :: Om1ab(nVVab) double precision,intent(in) :: X1ab(nVVab,nVVab) @@ -41,16 +41,16 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n double precision,intent(in) :: Om2ab(nOOab) double precision,intent(in) :: X2ab(nVVab,nOOab) double precision,intent(in) :: Y2ab(nOOab,nOOab) - double precision,intent(in) :: rho1ab(nBas,nBas,nVVab) - double precision,intent(in) :: rho2ab(nBas,nBas,nOOab) + double precision,intent(in) :: rho1ab(nOrb,nOrb,nVVab) + double precision,intent(in) :: rho2ab(nOrb,nOrb,nOOab) double precision,intent(in) :: Om1aa(nVVaa) double precision,intent(in) :: X1aa(nVVaa,nVVaa) double precision,intent(in) :: Y1aa(nOOaa,nVVaa) double precision,intent(in) :: Om2aa(nOOaa) double precision,intent(in) :: X2aa(nVVaa,nOOaa) double precision,intent(in) :: Y2aa(nOOaa,nOOaa) - double precision,intent(in) :: rho1aa(nBas,nBas,nVVaa) - double precision,intent(in) :: rho2aa(nBas,nBas,nOOaa) + double precision,intent(in) :: rho1aa(nOrb,nOrb,nVVaa) + double precision,intent(in) :: rho2aa(nOrb,nOrb,nOOaa) ! Local variables @@ -100,16 +100,16 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n allocate(Bpp(nVVab,nOOab),Cpp(nVVab,nVVab),Dpp(nOOab,nOOab)) - if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,ERI,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVab,1d0,eT,ERI,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOab,1d0,eT,ERI,Dpp) + if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOab,nVVab,1d0,ERI,Bpp) + call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVab,1d0,eT,ERI,Cpp) + call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOab,1d0,eT,ERI,Dpp) call ppLR(TDA_T,nOOab,nVVab,Bpp,Cpp,Dpp,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin)) deallocate(Bpp,Cpp,Dpp) - call RGTpp_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TAab) - if(.not.TDA) call RGTpp_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TBab) + call RGTpp_phBSE_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,nOOab,nVVab,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TAab) + if(.not.TDA) call RGTpp_phBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,nOOab,nVVab,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TBab) !----------------------------------------! ! Compute T-matrix for alpha-alpha block ! @@ -120,16 +120,16 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n allocate(Bpp(nVVaa,nOOaa),Cpp(nVVaa,nVVaa),Dpp(nOOaa,nOOaa)) - if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,ERI,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVaa,1d0,eT,ERI,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOaa,1d0,eT,ERI,Dpp) + if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOaa,nVVaa,1d0,ERI,Bpp) + call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVaa,1d0,eT,ERI,Cpp) + call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOaa,1d0,eT,ERI,Dpp) call ppLR(TDA_T,nOOaa,nVVaa,Bpp,Cpp,Dpp,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin)) deallocate(Bpp,Cpp,Dpp) - call RGTpp_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,nOOaa,nVVaa,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TAaa) - if(.not.TDA) call RGTpp_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,nOOaa,nVVaa,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TBaa) + call RGTpp_phBSE_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,nOOaa,nVVaa,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TAaa) + if(.not.TDA) call RGTpp_phBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,nOOaa,nVVaa,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TBaa) !------------------! ! Singlet manifold ! @@ -141,8 +141,8 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n ! Compute BSE singlet excitation energies - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph) - if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph) + if(.not.TDA) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) Aph(:,:) = Aph(:,:) + TAaa(:,:) + TAab(:,:) if(.not.TDA) Bph(:,:) = Bph(:,:) + TBaa(:,:) + TBab(:,:) @@ -150,12 +150,12 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE) call print_excitation_energies('phBSE@GTpp','singlet',nS,OmBSE) - call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE) + call phLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE) ! Compute dynamic correction for BSE via renormalized perturbation theory if(dBSE) & - call RGTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & + call RGTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & Om1ab,Om2ab,Om1aa,Om2aa,rho1ab,rho2ab,rho1aa,rho2aa,eT,eGT, & dipole_int,OmBSE,XpY_BSE,XmY_BSE,TAab,TAaa) @@ -171,8 +171,8 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n ! Compute BSE triplet excitation energies - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph) - if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph) + if(.not.TDA) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) Aph(:,:) = Aph(:,:) + TAaa(:,:) - TAab(:,:) if(.not.TDA) Bph(:,:) = Bph(:,:) + TBaa(:,:) - TBab(:,:) @@ -180,12 +180,12 @@ subroutine RGTpp_phBSE(exchange_kernel,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,n call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE) call print_excitation_energies('phBSE@GTpp','triplet',nS,OmBSE) - call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE) + call phLR_transition_vectors(.false.,nOrb,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE) ! Compute dynamic correction for BSE via renormalized perturbation theory if(dBSE) & - call RGTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & + call RGTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & Om1ab,Om2ab,Om1aa,Om2aa,rho1ab,rho2ab,rho1aa,rho2aa,eT,eGT, & dipole_int,OmBSE,XpY_BSE,XmY_BSE,TAab,TAaa) end if diff --git a/src/GW/ufBSE.f90 b/src/GW/RGW_phBSE_upfolded.f90 similarity index 88% rename from src/GW/ufBSE.f90 rename to src/GW/RGW_phBSE_upfolded.f90 index 592a089..5e1f61a 100644 --- a/src/GW/ufBSE.f90 +++ b/src/GW/RGW_phBSE_upfolded.f90 @@ -1,6 +1,6 @@ -subroutine ufBSE(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) +subroutine RGW_phBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) -! Unfold BSE@GW equations +! Upfolded phBSE@GW (TDA singlets only!) implicit none include 'parameters.h' @@ -43,9 +43,9 @@ subroutine ufBSE(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) ! Hello world write(*,*) - write(*,*)'**********************************************' - write(*,*)'| Unfolded BSE@GW calculation |' - write(*,*)'**********************************************' + write(*,*)'*********************************' + write(*,*)'* Upfolded phBSE@GW Calculation *' + write(*,*)'*********************************' write(*,*) ! TDA for W @@ -71,11 +71,11 @@ subroutine ufBSE(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) ! Compute BSE supermatrix ! !---------------------------! ! ! -! | A -Ve -Vh | ! +! | A Jph Kph | ! ! | | ! -! H = | +Vh C2h2p 0 | ! +! H = | Kph C2h2p 0 | ! ! | | ! -! | +Ve 0 C2p2h | ! +! | Jph 0 C2p2h | ! ! ! !---------------------------! @@ -102,21 +102,19 @@ subroutine ufBSE(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) end do end do - !----------------! - ! Blocks Vp & Ve ! - !----------------! + !------------------! + ! Blocks Jph & Kph ! + !------------------! - iajb=0 + ia = 0 do i=nC+1,nO do a=nO+1,nOrb-nR - do j=nC+1,nO - do b=nO+1,nOrb-nR - iajb = iajb + 1 + ia = ia + 1 - kc = 0 - do k=nC+1,nO - do c=nO+1,nOrb-nR - kc = kc + 1 + do k=nC+1,nO + do c=nO+1,nOrn-nR + do m=1,nS + kcm = kcm + 1 tmp = sqrt(2d0)*Kronecker_delta(k,j)*ERI(b,a,c,i) H(n1h1p+iajb,kc ) = +tmp diff --git a/src/GW/RGW_ppBSE_dynamic_kernel_D.f90 b/src/GW/RGW_ppBSE_dynamic_kernel_D.f90 index d9d61c1..e56fa80 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(i) + dem = - OmBSE + eGW(k) - 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(k) - Om(m) + eGW(j) + dem = - OmBSE + eGW(k) - 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) 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(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(i) + dem = - OmBSE + eGW(l) - 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) @@ -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(i) + dem = - OmBSE + eGW(k) - 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(k) - Om(m) + eGW(j) + dem = - OmBSE + eGW(k) - 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) 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(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(i) + dem = - OmBSE + eGW(l) - 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) diff --git a/src/GW/RGW_ppBSE_upfolded.f90 b/src/GW/RGW_ppBSE_upfolded.f90 new file mode 100644 index 0000000..3a00333 --- /dev/null +++ b/src/GW/RGW_ppBSE_upfolded.f90 @@ -0,0 +1,227 @@ +subroutine RGW_ppBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) + +! Upfolded ppBSE@GW (TDA triplets only!) + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nOrb + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: eGW(nOrb) + +! Local variables + + integer :: s + integer :: i,j,k,l + integer :: a,b,c,d + integer :: ia,jb,kc,iajb,kcld + integer,parameter :: maxH = 20 + double precision :: tmp + + integer :: n1h1p,n2h2p,nH + double precision,external :: Kronecker_delta + + integer,allocatable :: order(:) + double precision,allocatable :: H(:,:) + double precision,allocatable :: X(:,:) + double precision,allocatable :: Om(:) + double precision,allocatable :: Z(:) + +! Output variables + +! Hello world + + write(*,*) + write(*,*)'**********************************************' + write(*,*)'| Unfolded ppBSE@GW calculation |' + write(*,*)'**********************************************' + write(*,*) + +! TDA for W + + write(*,*) 'Tamm-Dancoff approximation by default!' + write(*,*) + +! Dimension of the supermatrix + + n1h1p = nO*nV + n2h2p = nO*nO*nV*nV + nH = n1h1p + n2h2p + n2h2p + +! Memory allocation + + allocate(order(nH),H(nH,nH),X(nH,nH),Om(nH),Z(nH)) + +! Initialization + + H(:,:) = 0d0 + +!----------------------------------------! +! Compute BSE supermatrix ! +!----------------------------------------! +! ! +! | D -M1 -M1 -M2 -M2 | ! +! | | ! +! | +M1 E1 0 0 0 | ! +! | | ! +! H = | +M1 0 E2 0 0 | ! +! | | ! +! | +M2 0 0 E3 0 | ! +! | | ! +! | +M2 0 0 0 E4 | ! +! ! +!----------------------------------------! + + !---------! + ! Block D ! + !---------! + + ij = 0 + do i=nC+1,nO + do j=i+1,nO + ij = ij + 1 + + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + + H(ij,kl) = - (eGW(i) + eGW(j))*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + + (ERI(i,j,k,l) - ERI(i,j,l,k)) + + end do + end do + + end do + end do + + !----------------! + ! Blocks M1 ! + !----------------! + + ijm = 0 + do i=nC+1,nO + do j=i+1,nO + do m=1,nS + ijm = ijm + 1 + + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + + tmp = sqrt(2d0)*Kronecker_delta(k,j)*M(i,k,m) + H(n2h+ijm,kl ) = +tmp + H(kl ,n2h+ijm) = -tmp + + tmp = sqrt(2d0)*Kronecker_delta(k,j)*M(i,k,m) + H(n2h+1*n3h1p+ijm,kl ) = +tmp + H(kl ,n2h+n3h1p+ijm) = -tmp + + tmp = sqrt(2d0)*Kronecker_delta(b,c)*M(j,k,m) + H(n2h+2*n3h1p+iajb,kc ) = +tmp + H(kc ,n1h1p+iajb) = -tmp + + tmp = sqrt(2d0)*Kronecker_delta(b,c)*M(j,k,m) + H(n2h+3*n2h2p+iajb,kc ) = +tmp + H(kc ,n1h1p+iajb) = -tmp + + end do + end do + + end do + end do + end do + end do + + !-------------! + ! Block 2h2p ! + !-------------! + + iajb = 0 + do i=nC+1,nO + do a=nO+1,nOrb-nR + do j=nC+1,nO + do b=nO+1,nOrb-nR + iajb = iajb + 1 + + kcld = 0 + do k=nC+1,nO + do c=nO+1,nOrb-nR + do l=nC+1,nO + do d=nO+1,nOrb-nR + kcld = kcld + 1 + + tmp = ((eHF(a) + eGW(b) - eHF(i) - eGW(j))*Kronecker_delta(i,k)*Kronecker_delta(a,c) & + + 2d0*ERI(a,k,i,c))*Kronecker_delta(j,l)*Kronecker_delta(b,d) + H(n1h1p +iajb,n1h1p +kcld) = tmp + H(n1h1p+n2h2p+iajb,n1h1p+n2h2p+kcld) = tmp + + end do + end do + end do + end do + + end do + end do + end do + end do + +!-------------------------! +! Diagonalize supermatrix ! +!-------------------------! + +! call matout(nH,nH,H) + + call diagonalize_general_matrix(nH,H,Om,X) + + do s=1,nH + order(s) = s + end do + + call quick_sort(Om,order,nH) + call set_order(X,order,nH,nH) + +!-----------------! +! Compute weights ! +!-----------------! + + Z(:) = 0d0 + do s=1,nH + do ia=1,n1h1p + Z(s) = Z(s) + X(ia,s)**2 + end do + end do + +!--------------! +! Dump results ! +!--------------! + + write(*,*)'-------------------------------------------' + write(*,*)' unfolded BSE excitation energies (eV) ' + write(*,*)'-------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') & + '|','#','|','Omega (eV)','|','Z','|' + write(*,*)'-------------------------------------------' + + do s=1,min(nH,maxH) + if(Z(s) > 1d-7) & + write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & + '|',s,'|',Om(s)*HaToeV,'|',Z(s),'|' + end do + + write(*,*)'-------------------------------------------' + write(*,*) + +end subroutine