diff --git a/src/GW/RGW_phBSE.f90 b/src/GW/RGW_phBSE.f90 index edb4111..1438194 100644 --- a/src/GW/RGW_phBSE.f90 +++ b/src/GW/RGW_phBSE.f90 @@ -1,5 +1,5 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta, & - nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) + nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) ! Compute the Bethe-Salpeter excitation energies @@ -18,16 +18,16 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple 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 integer,intent(in) :: nR integer,intent(in) :: nS - double precision,intent(in) :: eW(nBas) - double precision,intent(in) :: eGW(nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: eW(nOrb) + double precision,intent(in) :: eGW(nOrb) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart) ! Local variables @@ -61,10 +61,19 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple ! Memory allocation - allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), & + allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nOrb,nOrb,nS), & Aph(nS,nS),Bph(nS,nS),KA_sta(nS,nS),KB_sta(nS,nS), & OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS)) +!-----! +! TDA ! +!-----! + + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated in phBSE!' + write(*,*) + end if + ! Initialization EcBSE(:) = 0d0 @@ -76,23 +85,14 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple isp_W = 1 EcRPA = 0d0 - call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph) - if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph) + if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - call RGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta) - call RGW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta) - -!-----! -! TDA ! -!-----! - - if(TDA) then - write(*,*) 'Tamm-Dancoff approximation activated in phBSE!' - write(*,*) - end if + call RGW_phBSE_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta) + call RGW_phBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta) !------------------- ! Singlet manifold @@ -104,23 +104,23 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple ! Compute BSE excitation energies - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,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,eGW,ERI,Aph) + if(.not.TDA) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) ! Second-order BSE static kernel if(dophBSE2) then - allocate(W(nBas,nBas,nBas,nBas)) + allocate(W(nOrb,nOrb,nOrb,nOrb)) write(*,*) write(*,*) '*** Second-order BSE static kernel activated! ***' write(*,*) - call RGW_phBSE_static_kernel(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,W) - call RGW_phBSE2_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,eW,W,KA_sta) + call RGW_phBSE_static_kernel(eta,nOrb,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,W) + call RGW_phBSE2_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,1d0,eW,W,KA_sta) - if(.not.TDA) call RGW_phBSE2_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,eW,W,KB_sta) + if(.not.TDA) call RGW_phBSE2_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,1d0,eW,W,KB_sta) deallocate(W) @@ -132,22 +132,22 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE) call print_excitation_energies('phBSE@GW@RHF','singlet',nS,OmBSE) - call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE) - - !--------------------! - ! Cumulant expansion ! - !--------------------! - - ! call RGWC(.false.,nBas,nC,nO,nR,nS,OmBSE,rho_RPA,eGW) + call phLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE) !----------------------------------------------------! ! Compute the dynamical screening at the phBSE level ! !----------------------------------------------------! if(dBSE) & - call RGW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, & + call RGW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nOrb,nC,nO,nV,nR,nS,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, & OmBSE,XpY_BSE,XmY_BSE,KA_sta,KB_sta) + !----------------! + ! Upfolded phBSE ! + !----------------! + + call RGW_phBSE_upfolded(ispin,nOrb,nOrb,nC,nO,nV,nR,nS,ERI,rho_RPA,OmRPA,eGW) + end if !------------------- @@ -160,8 +160,8 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple ! Compute BSE excitation energies - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,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,eGW,ERI,Aph) + if(.not.TDA) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) Aph(:,:) = Aph(:,:) + KA_sta(:,:) if(.not.TDA) Bph(:,:) = Bph(:,:) + KB_sta(:,:) @@ -169,16 +169,22 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE) call print_excitation_energies('phBSE@GW@RHF','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 the dynamical screening at the BSE level !------------------------------------------------- if(dBSE) & - call RGW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, & + call RGW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nOrb,nC,nO,nV,nR,nS,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, & OmBSE,XpY_BSE,XmY_BSE,KA_sta,KB_sta) + !----------------! + ! Upfolded phBSE ! + !----------------! + + call RGW_phBSE_upfolded(ispin,nOrb,nOrb,nC,nO,nV,nR,nS,ERI,rho_RPA,OmRPA,eGW) + end if ! Scale properly correlation energy if exchange is included in interaction kernel diff --git a/src/GW/RGW_phBSE_upfolded.f90 b/src/GW/RGW_phBSE_upfolded.f90 index 5e1f61a..6bb1342 100644 --- a/src/GW/RGW_phBSE_upfolded.f90 +++ b/src/GW/RGW_phBSE_upfolded.f90 @@ -1,12 +1,13 @@ -subroutine RGW_phBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) +subroutine RGW_phBSE_upfolded(ispin,nBas,nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW) -! Upfolded phBSE@GW (TDA singlets only!) +! Upfolded phBSE@GW (TDA only!) implicit none include 'parameters.h' ! Input variables + integer,intent(in) :: ispin integer,intent(in) :: nBas integer,intent(in) :: nOrb integer,intent(in) :: nC @@ -14,10 +15,9 @@ subroutine RGW_phBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) 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) :: rho(nOrb,nOrb,nS) + double precision,intent(in) :: Om(nS) double precision,intent(in) :: eGW(nOrb) ! Local variables @@ -25,17 +25,17 @@ subroutine RGW_phBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) integer :: s integer :: i,j,k,l integer :: a,b,c,d - integer :: ia,jb,kc,iajb,kcld + integer :: m,ia,jb,iam,kcm integer,parameter :: maxH = 20 double precision :: tmp - integer :: n1h1p,n2h2p,nH + integer :: n1h,n1p,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 :: OmBSE(:) double precision,allocatable :: Z(:) ! Output variables @@ -55,13 +55,16 @@ subroutine RGW_phBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) ! Dimension of the supermatrix - n1h1p = nO*nV - n2h2p = nO*nO*nV*nV + n1h = nO + n1p = nV + + n1h1p = n1h*n1p + n2h2p = n1h1p*n1h1p nH = n1h1p + n2h2p + n2h2p ! Memory allocation - allocate(order(nH),H(nH,nH),X(nH,nH),Om(nH),Z(nH)) + allocate(order(nH),H(nH,nH),X(nH,nH),OmBSE(nH),Z(nH)) ! Initialization @@ -79,28 +82,59 @@ subroutine RGW_phBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) ! ! !---------------------------! - !---------! - ! Block A ! - !---------! + !----------------------! + ! Block A for singlets ! + !----------------------! - ia = 0 - do i=nC+1,nO - do a=nO+1,nOrb-nR - ia = ia + 1 - - jb = 0 - do j=nC+1,nO - do b=nO+1,nOrb-nR - jb = jb + 1 - - H(ia,jb) = (eGW(a) - eGW(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & - + 2d0*ERI(i,b,a,j) - ERI(i,b,j,a) + if(ispin == 1) then + ia = 0 + do i=nC+1,nO + do a=nO+1,nOrb-nR + ia = ia + 1 + + jb = 0 + do j=nC+1,nO + do b=nO+1,nOrb-nR + jb = jb + 1 + + H(ia,jb) = (eGW(a) - eGW(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + + 2d0*ERI(i,b,a,j) - ERI(i,b,j,a) + + end do end do + end do - end do - end do + + end if + + !----------------------! + ! Block A for triplets ! + !----------------------! + + if(ispin == 2) then + + ia = 0 + do i=nC+1,nO + do a=nO+1,nOrb-nR + ia = ia + 1 + + jb = 0 + do j=nC+1,nO + do b=nO+1,nOrb-nR + jb = jb + 1 + + H(ia,jb) = (eGW(a) - eGW(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + - ERI(i,b,j,a) + + end do + end do + + end do + end do + + end if !------------------! ! Blocks Jph & Kph ! @@ -111,24 +145,28 @@ subroutine RGW_phBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) do a=nO+1,nOrb-nR ia = ia + 1 + kcm = 0 do k=nC+1,nO - do c=nO+1,nOrn-nR + do c=nO+1,nOrb-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 - H(kc ,n1h1p+n2h2p+iajb) = -tmp + ! Jph + + tmp = - sqrt(2d0)*Kronecker_delta(a,c)*rho(i,k,m) + H(ia ,n1h1p+kcm) = tmp + H(n1h1p+n2h2p+kcm,ia) = tmp + + ! Kph - tmp = sqrt(2d0)*Kronecker_delta(b,c)*ERI(a,k,i,j) - H(n1h1p+n2h2p+iajb,kc ) = +tmp - H(kc ,n1h1p+iajb) = -tmp + tmp = sqrt(2d0)*Kronecker_delta(i,k)*rho(a,c,m) + H(ia ,n1h1p+n2h2p+kcm) = tmp + H(n1h1p+kcm,ia) = tmp - end do end do - end do end do + end do end do @@ -136,31 +174,16 @@ subroutine RGW_phBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) ! Block 2h2p ! !-------------! - iajb = 0 + iam = 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 + do m=1,nS + iam = iam + 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 = Om(m) + eGW(a) - eGW(i) + H(n1h1p +iam,n1h1p +iam) = tmp + H(n1h1p+n2h2p+iam,n1h1p+n2h2p+iam) = tmp - 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 @@ -169,15 +192,13 @@ subroutine RGW_phBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) ! Diagonalize supermatrix ! !-------------------------! -! call matout(nH,nH,H) - - call diagonalize_general_matrix(nH,H,Om,X) + call diagonalize_general_matrix(nH,H,OmBSE,X) do s=1,nH order(s) = s end do - call quick_sort(Om,order,nH) + call quick_sort(OmBSE,order,nH) call set_order(X,order,nH,nH) !-----------------! @@ -196,16 +217,16 @@ subroutine RGW_phBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) !--------------! write(*,*)'-------------------------------------------' - write(*,*)' unfolded BSE excitation energies (eV) ' + write(*,*)' Upfolded phBSE 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','|' + '|','#','|','OmBSE (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),'|' + '|',s,'|',OmBSE(s)*HaToeV,'|',Z(s),'|' end do write(*,*)'-------------------------------------------' diff --git a/src/GW/RGW_ppBSE.f90 b/src/GW/RGW_ppBSE.f90 index 808caa9..b64fd70 100644 --- a/src/GW/RGW_ppBSE.f90 +++ b/src/GW/RGW_ppBSE.f90 @@ -76,12 +76,24 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS double precision,intent(out) :: EcBSE(nspin) +!-----! +! TDA ! +!-----! + + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated in ppBSE!' + write(*,*) + end if + +! Initialization + + EcBSE(:) = 0d0 + !--------------------------------- ! Compute (singlet) RPA screening !--------------------------------- isp_W = 1 - EcRPA = 0d0 allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nOrb,nOrb,nS), & Aph(nS,nS),Bph(nS,nS)) @@ -108,7 +120,6 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS write(*,*) ispin = 1 - EcBSE(ispin) = 0d0 nOO = nO*(nO+1)/2 nVV = nV*(nV+1)/2 @@ -259,7 +270,6 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS write(*,*) ispin = 2 - EcBSE(ispin) = 0d0 nOO = nO*(nO-1)/2 nVV = nV*(nV-1)/2 @@ -354,6 +364,13 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS 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) + + !----------------! + ! Upfolded ppBSE ! + !----------------! + + call RGW_ppBSE_upfolded(ispin,nOrb,nC,nO,nV,nR,nS,ERI,rho_RPA,OmRPA,eGW) + deallocate(KB_sta,KC_sta,KD_sta) deallocate(Om1,X1,Y1,Om2,X2,Y2) diff --git a/src/GW/RGW_ppBSE_upfolded.f90 b/src/GW/RGW_ppBSE_upfolded.f90 index 3a00333..5e96d9d 100644 --- a/src/GW/RGW_ppBSE_upfolded.f90 +++ b/src/GW/RGW_ppBSE_upfolded.f90 @@ -1,23 +1,22 @@ -subroutine RGW_ppBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) +subroutine RGW_ppBSE_upfolded(ispin,nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW) -! Upfolded ppBSE@GW (TDA triplets only!) +! Upfolded ppBSE@GW (TDA only!) implicit none include 'parameters.h' ! Input variables - integer,intent(in) :: nBas + integer,intent(in) :: ispin 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) :: rho(nOrb,nOrb,nS) + double precision,intent(in) :: Om(nS) double precision,intent(in) :: eGW(nOrb) ! Local variables @@ -25,17 +24,17 @@ subroutine RGW_ppBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) integer :: s integer :: i,j,k,l integer :: a,b,c,d - integer :: ia,jb,kc,iajb,kcld + integer :: m,ij,kl,ijm integer,parameter :: maxH = 20 - double precision :: tmp + double precision :: tmp,tmp1,tmp2,tmp3,tmp4 - integer :: n1h1p,n2h2p,nH + integer :: n1h,n1p,n2h,n2p,n1h1p,n3h1p,n3p1h,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 :: OmBSE(:) double precision,allocatable :: Z(:) ! Output variables @@ -55,13 +54,33 @@ subroutine RGW_ppBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) ! Dimension of the supermatrix - n1h1p = nO*nV - n2h2p = nO*nO*nV*nV - nH = n1h1p + n2h2p + n2h2p + n1h = nO + n1p = nV + + if(ispin == 1) then + + n2h = nO*(nO+1)/2 + n2p = nV*(nV+1)/2 + + end if + + if(ispin == 2) then + + n2h = nO*(nO-1)/2 + n2p = nV*(nV-1)/2 + + end if + + n1h1p = n1h*n1p + + n3h1p = n2h*n1h1p + n3p1h = n2p*n1h1p + + nH = n1h1p + 4*n3h1p ! Memory allocation - allocate(order(nH),H(nH,nH),X(nH,nH),Om(nH),Z(nH)) + allocate(order(nH),H(nH,nH),X(nH,nH),OmBSE(nH),Z(nH)) ! Initialization @@ -83,31 +102,62 @@ subroutine RGW_ppBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) ! ! !----------------------------------------! - !---------! - ! Block D ! - !---------! + !----------------------! + ! Block D for singlets ! + !----------------------! - 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)) + if(ispin == 1) then + ij = 0 + do i=nC+1,nO + do j=i,nO + ij = ij + 1 + + kl = 0 + do k=nC+1,nO + do l=k,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))/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) + + end do end do + end do - end do - end do + + end if + + !----------------------! + ! Block D for triplets ! + !----------------------! + + if(ispin == 2) then + + 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 + + end if !----------------! - ! Blocks M1 ! + ! Blocks M1 & M2 ! !----------------! ijm = 0 @@ -121,59 +171,47 @@ subroutine RGW_ppBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) 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 + tmp1 = sqrt(2d0)*Kronecker_delta(j,l)*rho(i,k,m) + tmp2 = sqrt(2d0)*Kronecker_delta(j,l)*rho(i,l,m) + tmp3 = sqrt(2d0)*Kronecker_delta(i,l)*rho(j,k,m) + tmp4 = sqrt(2d0)*Kronecker_delta(i,k)*rho(j,l,m) - tmp = sqrt(2d0)*Kronecker_delta(b,c)*M(j,k,m) - H(n2h+3*n2h2p+iajb,kc ) = +tmp - H(kc ,n1h1p+iajb) = -tmp + H(n2h+0*n3h1p+ijm,kl ) = -tmp1 + H(kl ,n2h+0*n3h1p+ijm) = +tmp3 + + H(n2h+1*n3h1p+ijm,kl ) = -tmp1 + H(kl ,n2h+1*n3h1p+ijm) = +tmp4 + + H(n2h+2*n3h1p+ijm,kl ) = -tmp2 + H(kl ,n2h+2*n3h1p+ijm) = +tmp3 + + H(n2h+3*n3h1p+ijm,kl ) = -tmp2 + H(kl ,n2h+3*n3h1p+ijm) = +tmp4 - end do end do - end do + end do end do end do - !-------------! - ! Block 2h2p ! - !-------------! + !------------! + ! Block 3h1p ! + !------------! - iajb = 0 + ijm = 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 + do j=i+1,nO + do m=1,nS + ijm = ijm + 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 = - eGW(i) - eGW(j) + Om(m) - 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 + H(n2h+0*n3h1p+ijm,n2h+0*n3h1p+ijm) = tmp + H(n2h+1*n3h1p+ijm,n2h+1*n3h1p+ijm) = tmp + H(n2h+2*n3h1p+ijm,n2h+2*n3h1p+ijm) = tmp + H(n2h+3*n3h1p+ijm,n2h+3*n3h1p+ijm) = tmp - end do - end do - end do - end do - - end do end do end do end do @@ -182,15 +220,13 @@ subroutine RGW_ppBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) ! Diagonalize supermatrix ! !-------------------------! -! call matout(nH,nH,H) - - call diagonalize_general_matrix(nH,H,Om,X) + call diagonalize_general_matrix(nH,H,OmBSE,X) do s=1,nH order(s) = s end do - call quick_sort(Om,order,nH) + call quick_sort(OmBSE,order,nH) call set_order(X,order,nH,nH) !-----------------! @@ -199,8 +235,8 @@ subroutine RGW_ppBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) Z(:) = 0d0 do s=1,nH - do ia=1,n1h1p - Z(s) = Z(s) + X(ia,s)**2 + do ij=1,n2h + Z(s) = Z(s) + X(ij,s)**2 end do end do @@ -209,16 +245,16 @@ subroutine RGW_ppBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) !--------------! write(*,*)'-------------------------------------------' - write(*,*)' unfolded BSE excitation energies (eV) ' + write(*,*)' Upfolded ppBSE 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','|' + '|','#','|','OmBSE (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),'|' + '|',s,'|',OmBSE(s)*HaToeV,'|',Z(s),'|' end do write(*,*)'-------------------------------------------'