mirror of
https://github.com/pfloos/quack
synced 2025-01-08 20:33:19 +01:00
ok with linear response
This commit is contained in:
parent
ba40c9bdb8
commit
11794c6625
@ -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
|
||||
|
||||
|
||||
@ -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)
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
@ -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(*,*)'-------------------------------------------------------------------------------'
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
@ -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
|
||||
@ -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
|
||||
|
||||
Loading…
Reference in New Issue
Block a user