diff --git a/src/GT/GG0T0pp.f90 b/src/GT/GG0T0pp.f90 index b5453a8..7be6092 100644 --- a/src/GT/GG0T0pp.f90 +++ b/src/GT/GG0T0pp.f90 @@ -176,13 +176,12 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,T if(doppBSE) then - call GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV, & - ERI,dipole_int,eHF,eGT,EcBSE) + call GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV,ERI,dipole_int,eHF,eGT,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp@RHF correlation energy = ',EcBSE,' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp@RHF total energy = ',ENuc + EGHF + EcBSE,' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp@GHF correlation energy = ',EcBSE,' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp@GHF total energy = ',ENuc + EGHF + EcBSE,' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/GT/GGT_Tmatrix.f90 b/src/GT/GGT_Tmatrix.f90 index 909f2ab..3ef9984 100644 --- a/src/GT/GGT_Tmatrix.f90 +++ b/src/GT/GGT_Tmatrix.f90 @@ -46,21 +46,26 @@ subroutine GGT_Tmatrix(nOrb,nC,nO,nV,nR,nOO,nVV,lambda,ERI,eGF,Om1,rho1,Om2,rho2 do r=nC+1,nOrb-nR do q=nC+1,nOrb-nR do p=nC+1,nOrb-nR + T(p,q,r,s) = ERI(p,q,r,s) - ERI(p,q,s,r) cd = 0 - do c = nO+1, nOrb-nR - do d = c+1, nOrb-nR + do c=nO+1,nOrb-nR + do d=c+1,nOrb-nR cd = cd + 1 - T(p,q,r,s) = T(p,q,r,s) - rho1(p,q,cd) * rho1(r,s,cd) / Om1(cd) + + T(p,q,r,s) = T(p,q,r,s) - rho1(p,q,cd)*rho1(r,s,cd)/Om1(cd) + end do ! d end do ! c kl = 0 - do k = nC+1, nO - do l = k+1, nO + do k=nC+1,nO + do l=k+1,nO kl = kl + 1 - T(p,q,r,s) = T(p,q,r,s) + rho2(p,q,kl) * rho2(r,s,kl) / Om2(kl) + + T(p,q,r,s) = T(p,q,r,s) + rho2(p,q,kl)*rho2(r,s,kl)/Om2(kl) + enddo enddo diff --git a/src/GT/GGTpp_ppBSE.f90 b/src/GT/GGTpp_ppBSE.f90 index 5f8a861..dbc7e15 100644 --- a/src/GT/GGTpp_ppBSE.f90 +++ b/src/GT/GGTpp_ppBSE.f90 @@ -1,5 +1,4 @@ -subroutine GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV, & - ERI,dipole_int,eT,eGT,EcBSE) +subroutine GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV,ERI,dipole_int,eT,eGT,EcBSE) ! Compute the Bethe-Salpeter excitation energies with the T-matrix kernel @@ -30,7 +29,7 @@ subroutine GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV, & ! Local variables - double precision :: EcRPA(nspin) + double precision :: EcRPA double precision,allocatable :: Bpp(:,:),Cpp(:,:),Dpp(:,:) double precision,allocatable :: Om1(:), Om2(:) double precision,allocatable :: X1(:,:), X2(:,:) @@ -46,14 +45,13 @@ subroutine GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV, & !---------------------------------------------- ! Compute linear response !---------------------------------------------- + allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO)) allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)) - Bpp(:,:) = 0d0 - Cpp(:,:) = 0d0 - Dpp(:,:) = 0d0 - call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eT,ERI,Cpp) - call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eT,ERI,Dpp) + if(.not.TDA_T) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) + call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eT,ERI,Cpp) + call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eT,ERI,Dpp) call ppGLR(TDA_T,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) @@ -62,7 +60,9 @@ subroutine GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV, & !---------------------------------------------- ! Compute excitation densities !---------------------------------------------- + allocate(rho1(nOrb,nOrb,nVV),rho2(nOrb,nOrb,nOO)) + call GGTpp_excitation_density(nOrb,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2) deallocate(X1,Y1,X2,Y2) @@ -73,9 +73,9 @@ subroutine GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV, & allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)) - call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eGT,ERI,Cpp) - call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eGT,ERI,Dpp) if(.not.TDA_T) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) + call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eGT,ERI,Cpp) + call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eGT,ERI,Dpp) !---------------------------------------------- ! Compute T matrix tensor @@ -91,12 +91,10 @@ subroutine GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV, & allocate(KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO)) - call GGTpp_ppBSE_static_kernel_C(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT, & - Om1,rho1,Om2,rho2,T,KC_sta) - call GGTpp_ppBSE_static_kernel_D(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT, & - Om1,rho1,Om2,rho2,T,KD_sta) - if(.not.TDA_T) call GGTpp_ppBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT, & - Om1,rho1,Om2,rho2,T,KB_sta) + call GGTpp_ppBSE_static_kernel_C(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT,Om1,rho1,Om2,rho2,T,KC_sta) + call GGTpp_ppBSE_static_kernel_D(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT,Om1,rho1,Om2,rho2,T,KD_sta) + if(.not.TDA_T) & + call GGTpp_ppBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT,Om1,rho1,Om2,rho2,T,KB_sta) deallocate(Om1,Om2,rho1,rho2) ! Deallocate the 4-tensor T diff --git a/src/GT/GGTpp_self_energy_diag.f90 b/src/GT/GGTpp_self_energy_diag.f90 index c34bb73..1f34638 100644 --- a/src/GT/GGTpp_self_energy_diag.f90 +++ b/src/GT/GGTpp_self_energy_diag.f90 @@ -81,7 +81,7 @@ subroutine GGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rh do cd=1,nVV eps = e(i) + e(j) - Om1(cd) - num = rho1(i,j,cd)**2 + num = 0.5d0*rho1(i,j,cd)**2 EcGM = EcGM + num*eps/(eps**2 + eta**2) end do @@ -93,7 +93,7 @@ subroutine GGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rh do kl=1,nOO eps = e(a) + e(b) - Om2(kl) - num = rho2(a,b,kl)**2 + num = 0.5d0*rho2(a,b,kl)**2 EcGM = EcGM - num*eps/(eps**2 + eta**2) end do diff --git a/src/GT/RGTpp_self_energy_diag.f90 b/src/GT/RGTpp_self_energy_diag.f90 index 5a1fcc8..6138996 100644 --- a/src/GT/RGTpp_self_energy_diag.f90 +++ b/src/GT/RGTpp_self_energy_diag.f90 @@ -48,14 +48,14 @@ subroutine RGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1 do cd=1,nVVs eps = e(p) + e(i) - Om1s(cd) - num = (1d0/2d0)*rho1s(p,i,cd)**2 + num = 0.5d0*rho1s(p,i,cd)**2 Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do do cd=1,nVVt eps = e(p) + e(i) - Om1t(cd) - num = (3d0/2d0)*rho1t(p,i,cd)**2 + num = 1.5d0*rho1t(p,i,cd)**2 Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do @@ -72,14 +72,14 @@ subroutine RGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1 do kl=1,nOOs eps = e(p) + e(a) - Om2s(kl) - num = (1d0/2d0)*rho2s(p,a,kl)**2 + num = 0.5d0*rho2s(p,a,kl)**2 Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do do kl=1,nOOt eps = e(p) + e(a) - Om2t(kl) - num = (3d0/2d0)*rho2t(p,a,kl)**2 + num = 1.5d0*rho2t(p,a,kl)**2 Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do @@ -96,13 +96,13 @@ subroutine RGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1 do cd=1,nVVs eps = e(i) + e(j) - Om1s(cd) - num = (1d0/2d0)*rho1s(i,j,cd)**2 + num = 0.5d0*rho1s(i,j,cd)**2 EcGM = EcGM + num*eps/(eps**2 + eta**2) end do do cd=1,nVVt eps = e(i) + e(j) - Om1t(cd) - num = (3d0/2d0)*rho1t(i,j,cd)**2 + num = 1.5d0*rho1t(i,j,cd)**2 EcGM = EcGM + num*eps/(eps**2 + eta**2) end do @@ -114,13 +114,13 @@ subroutine RGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1 do kl=1,nOOs eps = e(a) + e(b) - Om2s(kl) - num = (1d0/2d0)*rho2s(a,b,kl)**2 + num = 0.5d0*rho2s(a,b,kl)**2 EcGM = EcGM - num*eps/(eps**2 + eta**2) end do do kl=1,nOOt eps = e(a) + e(b) - Om2t(kl) - num = (3d0/2d0)*rho2t(a,b,kl)**2 + num = 1.5d0*rho2t(a,b,kl)**2 EcGM = EcGM - num*eps/(eps**2 + eta**2) end do diff --git a/src/GT/print_GG0T0pp.f90 b/src/GT/print_GG0T0pp.f90 index 468c4c3..02787e1 100644 --- a/src/GT/print_GG0T0pp.f90 +++ b/src/GT/print_GG0T0pp.f90 @@ -1,4 +1,4 @@ -subroutine print_GG0T0pp(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) +subroutine print_GG0T0pp(nBas,nO,eHF,ENuc,EGHF,SigT,Z,eGT,EcGM,EcRPA) ! Print one-electron energies and other stuff for G0T0 @@ -8,7 +8,7 @@ subroutine print_GG0T0pp(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) integer,intent(in) :: nBas integer,intent(in) :: nO double precision,intent(in) :: ENuc - double precision,intent(in) :: ERHF + double precision,intent(in) :: EGHF double precision,intent(in) :: EcGM double precision,intent(in) :: EcRPA double precision,intent(in) :: eHF(nBas) @@ -48,14 +48,14 @@ subroutine print_GG0T0pp(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) end do write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@RHF HOMO energy = ',eGT(HOMO)*HaToeV,' eV' - write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@RHF LUMO energy = ',eGT(LUMO)*HaToeV,' eV' - write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@RHF HOMO-LUMO gap = ',Gap*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@GHF HOMO energy = ',eGT(HOMO)*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@GHF LUMO energy = ',eGT(LUMO)*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@GHF HOMO-LUMO gap = ',Gap*HaToeV,' eV' write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp@RHF correlation energy = ',EcRPA,' au' - write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp@RHF total energy = ',ENuc + ERHF + EcRPA,' au' - write(*,'(2X,A60,F15.6,A3)') ' GM@G0T0pp@RHF correlation energy = ',EcGM,' au' - write(*,'(2X,A60,F15.6,A3)') ' GM@G0T0pp@RHF total energy = ',ENuc + ERHF + EcGM,' au' + write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp@GHF correlation energy = ',EcRPA,' au' + write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp@GHF total energy = ',ENuc + EGHF + EcRPA,' au' + write(*,'(2X,A60,F15.6,A3)') ' GM@G0T0pp@GHF correlation energy = ',EcGM,' au' + write(*,'(2X,A60,F15.6,A3)') ' GM@G0T0pp@GHF total energy = ',ENuc + EGHF + EcGM,' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/GW/UG0W0.f90 b/src/GW/UG0W0.f90 index 44773a2..493f464 100644 --- a/src/GW/UG0W0.f90 +++ b/src/GW/UG0W0.f90 @@ -194,10 +194,10 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy (spin-conserved) = ',EcBSE(1),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy (spin-flip) = ',EcBSE(2),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy = ',sum(EcBSE),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@UHF total energy = ',ENuc + EUHF + sum(EcBSE),' au' + write(*,'(2X,A60,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy (spin-conserved) = ',EcBSE(1),' au' + write(*,'(2X,A60,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy (spin-flip) = ',EcBSE(2),' au' + write(*,'(2X,A60,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy = ',sum(EcBSE),' au' + write(*,'(2X,A60,F20.10,A3)') 'Tr@BSE@G0W0@UHF total energy = ',ENuc + EUHF + sum(EcBSE),' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) @@ -210,10 +210,10 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0@UHF correlation energy (spin-conserved) = ',EcBSE(1),' au' - write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0@UHF correlation energy (spin-flip) = ',EcBSE(2),' au' - write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0@UHF correlation energy = ',sum(EcBSE),' au' - write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0@UHF total energy = ',ENuc + EUHF + sum(EcBSE),' au' + write(*,'(2X,A60,F20.10,A3)') 'AC@BSE@G0W0@UHF correlation energy (spin-conserved) = ',EcBSE(1),' au' + write(*,'(2X,A60,F20.10,A3)') 'AC@BSE@G0W0@UHF correlation energy (spin-flip) = ',EcBSE(2),' au' + write(*,'(2X,A60,F20.10,A3)') 'AC@BSE@G0W0@UHF correlation energy = ',sum(EcBSE),' au' + write(*,'(2X,A60,F20.10,A3)') 'AC@BSE@G0W0@UHF total energy = ',ENuc + EUHF + sum(EcBSE),' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*)