diff --git a/src/GT/G0T0pp.f90 b/src/GT/G0T0pp.f90 index 5c4eac9..385014a 100644 --- a/src/GT/G0T0pp.f90 +++ b/src/GT/G0T0pp.f90 @@ -144,9 +144,9 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp ! Compute T-matrix version of the self-energy !---------------------------------------------- - EcGM = 0d0 + EcGM = 0d0 Sig(:) = 0d0 - Z(:) = 0d0 + Z(:) = 0d0 iblock = 3 @@ -285,9 +285,9 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp end if - call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,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_MO,eHF,eGT,EcAC) + call GTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,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_MO,eHF,eGT,EcAC) if(exchange_kernel) then diff --git a/src/GT/ACFDT_Tmatrix.f90 b/src/GT/GTpp_phACFTD.f90 similarity index 66% rename from src/GT/ACFDT_Tmatrix.f90 rename to src/GT/GTpp_phACFTD.f90 index 5953979..e15d00c 100644 --- a/src/GT/ACFDT_Tmatrix.f90 +++ b/src/GT/GTpp_phACFTD.f90 @@ -1,6 +1,6 @@ -subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & - nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & - Om2t,X2t,Y2t,rho1t,rho2t,ERI,eT,eGT,EcAC) +subroutine GTpp_phACFDT(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & + nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & + Om2t,X2t,Y2t,rho1t,rho2t,ERI,eT,eGT,EcAC) ! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem for the T-matrix @@ -63,13 +63,18 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple double precision,allocatable :: Ec(:,:) double precision :: EcRPA(nspin) + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) + double precision,allocatable :: Bpp(:,:) + double precision,allocatable :: Cpp(:,:) + double precision,allocatable :: Dpp(:,:) double precision,allocatable :: TAs(:,:) double precision,allocatable :: TBs(:,:) double precision,allocatable :: TAt(:,:) double precision,allocatable :: TBt(:,:) - double precision,allocatable :: Om(:,:) - double precision,allocatable :: XpY(:,:,:) - double precision,allocatable :: XmY(:,:,:) + double precision,allocatable :: Om(:) + double precision,allocatable :: XpY(:,:) + double precision,allocatable :: XmY(:,:) ! Output variables @@ -77,8 +82,7 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple ! Memory allocation - allocate(TAs(nS,nS),TBs(nS,nS),TAt(nS,nS),TBt(nS,nS), & - Om(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) + allocate(Aph(nS,nS),Bph(nS,nS),TAs(nS,nS),TBs(nS,nS),TAt(nS,nS),TBt(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS)) allocate(Ec(nAC,nspin)) ! Antisymmetrized kernel version @@ -118,7 +122,15 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple isp_T = 1 iblock = 3 - call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,eT,ERI,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T)) + allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) + + if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,ERI,Bpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,lambda,eT,ERI,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,lambda,eT,ERI,Dpp) + + call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T)) + + deallocate(Bpp,Cpp,Dpp) call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) @@ -128,7 +140,15 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple isp_T = 2 iblock = 4 - call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,eT,ERI,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T)) + allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) + + if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,ERI,Bpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,lambda,eT,ERI,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,lambda,eT,ERI,Dpp) + + call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T)) + + deallocate(Bpp,Cpp,Dpp) call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) @@ -137,10 +157,15 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple end if - call linear_response_BSE(ispin,.false.,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TAt+TAs,TBt+TBs, & - EcAC(ispin),Om(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,Aph) + if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph) + + Aph(:,:) = Aph(:,:) + TAt(:,:) + TAs(:,:) + if(.not.TDA) Bph(:,:) = Bph(:,:) + TBt(:,:) + TBs(:,:) + + call phLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY) - call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) + call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin)) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) @@ -181,7 +206,15 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple isp_T = 1 iblock = 3 - call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,eT,ERI,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T)) + allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) + + if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,ERI,Bpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,lambda,eT,ERI,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,lambda,eT,ERI,Dpp) + + call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T)) + + deallocate(Bpp,Cpp,Dpp) call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) @@ -191,19 +224,32 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple isp_T = 2 iblock = 4 - call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,eT,ERI,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T)) + allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) + + if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,ERI,Bpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,lambda,eT,ERI,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,lambda,eT,ERI,Dpp) + + call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T)) + + deallocate(Bpp,Cpp,Dpp) call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Om1t,rho1t,Om2t,rho2t,TAt) if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Om1t,rho1t,Om2t,rho2t,TBt) - end if + end if - call linear_response_BSE(ispin,.false.,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TAt-TAs,TBt-TBs, & - EcAC(ispin),Om(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,Aph) + if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph) - call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) + Aph(:,:) = Aph(:,:) + TAt(:,:) - TAs(:,:) + if(.not.TDA) Bph(:,:) = Bph(:,:) + TBt(:,:) - TBs(:,:) + + call phLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY) + + call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin)) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) diff --git a/src/GT/GTpp_phBSE.f90 b/src/GT/GTpp_phBSE.f90 index 6c81ae7..e0f5137 100644 --- a/src/GT/GTpp_phBSE.f90 +++ b/src/GT/GTpp_phBSE.f90 @@ -54,15 +54,24 @@ subroutine GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n ! Local variables + logical :: dRPA + integer :: ispin integer :: iblock + double precision,allocatable :: Bpp(:,:) + double precision,allocatable :: Cpp(:,:) + double precision,allocatable :: Dpp(:,:) + + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) + double precision :: EcRPA(nspin) double precision,allocatable :: TAab(:,:),TBab(:,:) double precision,allocatable :: TAaa(:,:),TBaa(:,:) - double precision,allocatable :: OmBSE(:,:) - double precision,allocatable :: XpY_BSE(:,:,:) - double precision,allocatable :: XmY_BSE(:,:,:) + double precision,allocatable :: OmBSE(:) + double precision,allocatable :: XpY_BSE(:,:) + double precision,allocatable :: XmY_BSE(:,:) ! Output variables @@ -70,8 +79,8 @@ subroutine GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n ! Memory allocation - allocate(TAab(nS,nS),TBab(nS,nS),TAaa(nS,nS),TBaa(nS,nS), & - OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin)) + allocate(Aph(nS,nS),Bph(nS,nS),TAab(nS,nS),TBab(nS,nS),TAaa(nS,nS),TBaa(nS,nS), & + OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS)) !---------------------------------------! ! Compute T-matrix for alpha-beta block ! @@ -80,7 +89,15 @@ subroutine GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n ispin = 1 iblock = 3 - call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,eT,ERI,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin)) + 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) + + call ppLR(TDA_T,nOOab,nVVab,Bpp,Cpp,Dpp,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin)) + + deallocate(Bpp,Cpp,Dpp) call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TAab) if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,1d0,Om1ab,rho1ab,Om2ab,rho2ab,TBab) @@ -92,7 +109,15 @@ subroutine GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n ispin = 2 iblock = 4 - call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,eT,ERI,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin)) + 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) + + call ppLR(TDA_T,nOOaa,nVVaa,Bpp,Cpp,Dpp,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin)) + + deallocate(Bpp,Cpp,Dpp) call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOaa,nVVaa,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TAaa) if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOaa,nVVaa,1d0,Om1aa,rho1aa,Om2aa,rho2aa,TBaa) @@ -104,15 +129,19 @@ subroutine GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n if(singlet) then ispin = 1 - EcBSE(ispin) = 0d0 ! Compute BSE singlet excitation energies - call linear_response_BSE(ispin,.false.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TAab+TAaa,TBab+TBaa, & - EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + 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 print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin)) - call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + Aph(:,:) = Aph(:,:) + TAaa(:,:) + TAab(:,:) + if(.not.TDA) Bph(:,:) = Bph(:,:) + TBaa(:,:) + TBab(:,:) + + call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE) + + call print_excitation('phBSE@GTpp ',ispin,nS,OmBSE) + call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE) if(dBSE) then @@ -121,13 +150,12 @@ subroutine GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n if(evDyn) then print*, ' Iterative dynamical correction for BSE@GT NYI' -! call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, & -! OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + else call Bethe_Salpeter_Tmatrix_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & Om1ab,Om2ab,Om1aa,Om2aa,rho1ab,rho2ab,rho1aa,rho2aa,eT,eGT, & - dipole_int,OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),TAab,TAaa) + dipole_int,OmBSE,XpY_BSE,XmY_BSE,TAab,TAaa) end if end if @@ -141,14 +169,19 @@ subroutine GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n if(triplet) then ispin = 2 - EcBSE(ispin) = 0d0 ! Compute BSE triplet excitation energies - call linear_response_BSE(ispin,.false.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TAaa-TAab,TBaa-TBab, & - EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) - call print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin)) - call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + 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) + + Aph(:,:) = Aph(:,:) + TAaa(:,:) - TAab(:,:) + if(.not.TDA) Bph(:,:) = Bph(:,:) + TBaa(:,:) - TBab(:,:) + + call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE) + + call print_excitation('phBSE@GTpp ',ispin,nS,OmBSE) + call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE) if(dBSE) then @@ -157,13 +190,12 @@ subroutine GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n if(evDyn) then print*, ' Iterative dynamical correction for BSE@GT NYI' -! call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, & -! OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + else call Bethe_Salpeter_Tmatrix_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & Om1ab,Om2ab,Om1aa,Om2aa,rho1ab,rho2ab,rho1aa,rho2aa,eT,eGT, & - dipole_int,OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),TAab,TAaa) + dipole_int,OmBSE,XpY_BSE,XmY_BSE,TAab,TAaa) end if end if diff --git a/src/GT/evGTpp.f90 b/src/GT/evGTpp.f90 index 3cf172c..1412d78 100644 --- a/src/GT/evGTpp.f90 +++ b/src/GT/evGTpp.f90 @@ -269,9 +269,9 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T end if - call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & - nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & - Om2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcAC) + call GTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & + nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & + Om2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcAC) if(exchange_kernel) then diff --git a/src/GT/qsGTpp.f90 b/src/GT/qsGTpp.f90 index 33ddb83..0ed0487 100644 --- a/src/GT/qsGTpp.f90 +++ b/src/GT/qsGTpp.f90 @@ -391,9 +391,9 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T end if - call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & - nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & - Om2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcAC) + call GTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & + nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & + Om2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GW/GW_ppBSE.f90 b/src/GW/GW_ppBSE.f90 index 152023e..83a9dd2 100644 --- a/src/GW/GW_ppBSE.f90 +++ b/src/GW/GW_ppBSE.f90 @@ -29,20 +29,29 @@ subroutine GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole integer :: ispin integer :: isp_W + logical :: dRPA = .false. + logical :: dRPA_W = .true. + integer :: nOO integer :: nVV + double precision,allocatable :: A_sta(:,:) + double precision,allocatable :: B_sta(:,:) + double precision,allocatable :: KA_sta(:,:) + double precision,allocatable :: KB_sta(:,:) + + double precision :: EcRPA double precision,allocatable :: OmRPA(:) double precision,allocatable :: XpY_RPA(:,:) double precision,allocatable :: XmY_RPA(:,:) double precision,allocatable :: rho_RPA(:,:,:) - double precision,allocatable :: Omega1(:) + double precision,allocatable :: Om1(:) double precision,allocatable :: X1(:,:) double precision,allocatable :: Y1(:,:) - double precision,allocatable :: Omega2(:) + double precision,allocatable :: Om2(:) double precision,allocatable :: X2(:,:) double precision,allocatable :: Y2(:,:) @@ -61,13 +70,15 @@ subroutine GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole isp_W = 1 EcRPA = 0d0 - ! Memory allocation + call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,A_sta) + if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta) - allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) - - call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call phLR(TDA_W,nS,A_sta,B_sta,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + call GW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta) + call GW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta) + !------------------- ! Singlet manifold !------------------- @@ -85,8 +96,8 @@ subroutine GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole nOO = nO*(nO+1)/2 nVV = nV*(nV+1)/2 - allocate(Omega1(nVV),X1(nVV,nVV),Y1(nOO,nVV), & - Omega2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & + allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV), & + Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & WB(nVV,nOO),WC(nVV,nVV),WD(nOO,nOO)) if(.not.TDA) call static_screening_WB_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WB) @@ -95,15 +106,15 @@ subroutine GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole ! Compute BSE excitation energies - call linear_response_pp_BSE(ispin,TDA,.true.,nBas,nC,nO,nV,nR,nOO,nVV,1d0,eGW,ERI,WB,WC,WD, & - Omega1,X1,Y1,Omega2,X2,Y2,EcBSE(ispin)) + call linear_response_pp_BSE(ispin,TDA,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,eGW,ERI,WB,WC,WD, & + Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin)) -! call print_excitation('pp-BSE (N+2)',ispin,nVV,Omega1) -! call print_excitation('pp-BSE (N-2)',ispin,nOO,Omega2) +! call print_excitation('pp-BSE (N+2)',ispin,nVV,Om1) +! call print_excitation('pp-BSE (N-2)',ispin,nOO,Om2) - call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2) + call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) - deallocate(Omega1,X1,Y1,Omega2,X2,Y2,WB,WC,WD) + deallocate(Om1,X1,Y1,Om2,X2,Y2,WB,WC,WD) end if @@ -124,8 +135,8 @@ subroutine GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole nOO = nO*(nO-1)/2 nVV = nV*(nV-1)/2 - allocate(Omega1(nVV),X1(nVV,nVV),Y1(nOO,nVV), & - Omega2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & + allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV), & + Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & WB(nVV,nOO),WC(nVV,nVV),WD(nOO,nOO)) if(.not.TDA) call static_screening_WB_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WB) @@ -134,15 +145,15 @@ subroutine GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole ! Compute BSE excitation energies - call linear_response_pp_BSE(ispin,TDA,.true.,nBas,nC,nO,nV,nR,nOO,nVV,1d0,eGW,ERI,WB,WC,WD, & - Omega1,X1,Y1,Omega2,X2,Y2,EcBSE(ispin)) + call linear_response_pp_BSE(ispin,TDA,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,eGW,ERI,WB,WC,WD, & + Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin)) -! call print_excitation('pp-BSE (N+2)',ispin,nVV,Omega1) -! call print_excitation('pp-BSE (N-2)',ispin,nOO,Omega2) +! call print_excitation('pp-BSE (N+2)',ispin,nVV,Om1) +! call print_excitation('pp-BSE (N-2)',ispin,nOO,Om2) - call print_transition_vectors_pp(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2) + call print_transition_vectors_pp(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) - deallocate(Omega1,X1,Y1,Omega2,X2,Y2,WB,WC,WD) + deallocate(Om1,X1,Y1,Om2,X2,Y2,WB,WC,WD) end if diff --git a/src/LR/linear_response_BSE.f90 b/src/LR/linear_response_BSE.f90 deleted file mode 100644 index 5aa3835..0000000 --- a/src/LR/linear_response_BSE.f90 +++ /dev/null @@ -1,107 +0,0 @@ -subroutine linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A_BSE,B_BSE,Ec,Omega,XpY,XmY) - -! Compute linear response with BSE additional terms - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: ispin - logical,intent(in) :: dRPA - logical,intent(in) :: TDA - logical,intent(in) :: BSE - double precision,intent(in) :: eta - integer,intent(in) :: nBas - integer,intent(in) :: nC - integer,intent(in) :: nO - integer,intent(in) :: nV - integer,intent(in) :: nR - integer,intent(in) :: nS - double precision,intent(in) :: lambda - double precision,intent(in) :: e(nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: A_BSE(nS,nS) - double precision,intent(in) :: B_BSE(nS,nS) - -! Local variables - - double precision :: trace_matrix - double precision,allocatable :: A(:,:) - double precision,allocatable :: B(:,:) - double precision,allocatable :: ApB(:,:) - double precision,allocatable :: AmB(:,:) - double precision,allocatable :: AmBSq(:,:) - double precision,allocatable :: AmBIv(:,:) - double precision,allocatable :: Z(:,:) - -! Output variables - - double precision,intent(out) :: Ec - double precision,intent(out) :: Omega(nS) - double precision,intent(out) :: XpY(nS,nS) - double precision,intent(out) :: XmY(nS,nS) - -! Memory allocation - - allocate(A(nS,nS),B(nS,nS),ApB(nS,nS),AmB(nS,nS),AmBSq(nS,nS),AmBIv(nS,nS),Z(nS,nS)) - -! Build A and B matrices - - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A) - - if(BSE) A(:,:) = A(:,:) - A_BSE(:,:) - -! Tamm-Dancoff approximation - - if(TDA) then - - B(:,:) = 0d0 - XpY(:,:) = A(:,:) - call diagonalize_matrix(nS,XpY,Omega) - XpY(:,:) = transpose(XpY(:,:)) - XmY(:,:) = XpY(:,:) - - else - - call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B) - - if(BSE) B(:,:) = B(:,:) - B_BSE(:,:) - - ! Build A + B and A - B matrices - - ApB = A + B - AmB = A - B - - ! Diagonalize linear response matrix - - call diagonalize_matrix(nS,AmB,Omega) - - if(minval(Omega) < 0d0) & - call print_warning('You may have instabilities in linear response: A-B is not positive definite!!') - - call ADAt(nS,AmB,1d0*sqrt(Omega),AmBSq) - call ADAt(nS,AmB,1d0/sqrt(Omega),AmBIv) - - Z = matmul(AmBSq,matmul(ApB,AmBSq)) - - call diagonalize_matrix(nS,Z,Omega) - - if(minval(Omega) < 0d0) & - call print_warning('You may have instabilities in linear response: negative excitations!!') - - Omega = sqrt(Omega) - - XpY = matmul(transpose(Z),AmBSq) - call DA(nS,1d0/sqrt(Omega),XpY) - - XmY = matmul(transpose(Z),AmBIv) - call DA(nS,1d0*sqrt(Omega),XmY) - - end if - - ! Compute the RPA correlation energy - - Ec = 0.5d0*(sum(Omega) - trace_matrix(nS,A)) - -end subroutine diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index ef566ed..a0ee3dd 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -176,8 +176,6 @@ program QuAcK !------------------------------------------------------------------------ call read_basis_pyscf (nBas,nO,nV) -! call read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, & -! max_ang_mom,min_exponent,max_exponent) nS(:) = (nO(:) - nC(:))*(nV(:) - nR(:)) !------------------------------------------------------------------------ @@ -926,8 +924,6 @@ program QuAcK write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0W0 = ',t_GW,' seconds' write(*,*) -! if(dophBSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,epsHF,eG0W0) - end if !------------------------------------------------------------------------ @@ -938,15 +934,12 @@ program QuAcK call cpu_time(start_GW) call ufGW(nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,epsHF) -! call CCGW(maxSCF_CC,thresh_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF) call cpu_time(end_GW) t_GW = end_GW - start_GW write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufGW = ',t_GW,' seconds' write(*,*) -! if(dophBSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,epsHF,eG0W0) - end if !------------------------------------------------------------------------ @@ -959,7 +952,6 @@ program QuAcK if(unrestricted) then - !print*,'!!! G0T0 NYI at the unrestricted level !!!' call UG0T0(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,evDyn, & spin_conserved,spin_flip,linGT,eta_GT,regGT,nBas,nC,nO,nV, & nR,nS,ENuc,EHF,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & @@ -967,7 +959,6 @@ program QuAcK else -! call soG0T0(eta_GT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI_MO,epsHF) call G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,evDyn,doppBSE,singlet,triplet, & linGT,eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,epsHF,Vxc) @@ -1148,4 +1139,4 @@ program QuAcK write(*,'(A65,1X,F9.3,A8)') 'Total wall time for QuAcK = ',t_QuAcK,' seconds' write(*,*) -end program QuAcK +end program