LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:35:36 +01:00
Code Releases Activity

Merge branch 'master' of github.com:pfloos/QuAcK

Browse Source
This commit is contained in:
Antoine Marie 2024-10-23 18:07:04 +02:00
commit 7c972806a0
11 changed files with 166 additions and 140 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
mol/H2.xyz
View File

@ -1,4 +1,4 @@
2 2
H 0.00000000 0.0 0.0 H 0.00000000 0.0 0.0
H 0.00000000 0.0 0.74031892 H 0.00000000 0.0 1.235

2
src/GW/RGW.f90
View File

@ -162,7 +162,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_dii
! Perform CC-based G0W0 calculation ! Perform CC-based G0W0 calculation
!------------------------------------------------------------------------ !------------------------------------------------------------------------
doccG0W0 = .true. doccG0W0 = .false.
if(doccG0W0) then if(doccG0W0) then

4
src/GW/RGW_phBSE.f90
View File

@ -65,9 +65,7 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
Aph(nS,nS),Bph(nS,nS),KA_sta(nS,nS),KB_sta(nS,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)) OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS))
!-----! ! TDA
! TDA !
!-----!
if(TDA) then if(TDA) then
write(*,*) 'Tamm-Dancoff approximation activated in phBSE!' write(*,*) 'Tamm-Dancoff approximation activated in phBSE!'

18
src/GW/UG0W0.f90
View File

@ -46,9 +46,9 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
! Local variables ! Local variables
logical :: print_W = .true. logical :: print_W = .true.
logical :: dRPA logical :: dRPA_W
integer :: is integer :: is
integer :: ispin integer :: isp_W
double precision :: flow double precision :: flow
double precision :: EcRPA double precision :: EcRPA
double precision :: EcGM(nspin) double precision :: EcGM(nspin)
@ -78,7 +78,7 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
! Initialization ! Initialization
EcRPA = 0d0 EcRPA = 0d0
dRPA = .true. dRPA_W = .true.
! TDA for W ! TDA for W
@ -113,10 +113,10 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
! Spin-conserving transitions ! Spin-conserving transitions
ispin = 1 isp_W = 1
call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) if(.not.TDA_W) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
call phULR(TDA_W,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY) call phULR(TDA_W,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY)
@ -172,8 +172,8 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
! Compute RPA correlation energy ! Compute RPA correlation energy
call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) if(.not.TDA_W) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
call phULR(TDA_W,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY) call phULR(TDA_W,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY)
@ -205,7 +205,7 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
if(doACFDT) then if(doACFDT) then
call UGW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,dophBSE,spin_conserved,spin_flip,eta, & call UGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,spin_conserved,spin_flip,eta, &
nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGW,EcBSE) nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGW,EcBSE)
write(*,*) write(*,*)

85
src/GW/UGW_phACFDT.f90
View File

@ -1,5 +1,5 @@
subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin_flip,eta, & subroutine UGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,spin_conserved,spin_flip,eta, &
nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eW,e,EcAC) nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eW,eGW,EcAC)
! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem ! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem
@ -10,10 +10,8 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
! Input variables ! Input variables
logical,intent(in) :: doXBS logical,intent(in) :: doXBS
logical,intent(in) :: dRPA
logical,intent(in) :: TDA_W logical,intent(in) :: TDA_W
logical,intent(in) :: TDA logical,intent(in) :: TDA
logical,intent(in) :: BSE
logical,intent(in) :: exchange_kernel logical,intent(in) :: exchange_kernel
logical,intent(in) :: spin_conserved logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip logical,intent(in) :: spin_flip
@ -26,19 +24,27 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
integer,intent(in) :: nR(nspin) integer,intent(in) :: nR(nspin)
integer,intent(in) :: nS(nspin) integer,intent(in) :: nS(nspin)
double precision,intent(in) :: eW(nBas,nspin) double precision,intent(in) :: eW(nBas,nspin)
double precision,intent(in) :: e(nBas,nspin) double precision,intent(in) :: eGW(nBas,nspin)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
! Local variables ! Local variables
logical :: dRPA = .false.
logical :: dRPA_W = .true.
integer :: ispin integer :: ispin
integer :: isp_W integer :: isp_W
integer :: iAC integer :: iAC
double precision :: lambda double precision :: lambda
double precision,allocatable :: Ec(:,:) double precision,allocatable :: Ec(:,:)
double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:)
double precision,allocatable :: KA(:,:)
double precision,allocatable :: KB(:,:)
double precision :: EcRPA double precision :: EcRPA
double precision,allocatable :: OmRPA(:) double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:) double precision,allocatable :: XpY_RPA(:,:)
@ -46,14 +52,10 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
double precision,allocatable :: rho_RPA(:,:,:,:) double precision,allocatable :: rho_RPA(:,:,:,:)
integer :: nS_aa,nS_bb,nS_sc integer :: nS_aa,nS_bb,nS_sc
double precision,allocatable :: Om_sc(:)
double precision,allocatable :: XpY_sc(:,:)
double precision,allocatable :: XmY_sc(:,:)
integer :: nS_ab,nS_ba,nS_sf integer :: nS_ab,nS_ba,nS_sf
double precision,allocatable :: Om_sf(:) double precision,allocatable :: Om(:)
double precision,allocatable :: XpY_sf(:,:) double precision,allocatable :: XpY(:,:)
double precision,allocatable :: XmY_sf(:,:) double precision,allocatable :: XmY(:,:)
! Output variables ! Output variables
@ -63,7 +65,6 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
allocate(Ec(nAC,nspin)) allocate(Ec(nAC,nspin))
! Hello World ! Hello World
write(*,*) '-----------------------------------------------------------' write(*,*) '-----------------------------------------------------------'
@ -109,18 +110,24 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
nS_sf = nS_ab + nS_ba nS_sf = nS_ab + nS_ba
allocate(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin)) allocate(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin))
allocate(Aph(nS_sc,nS_sc),Bph(nS_sc,nS_sc),KA(nS_sc,nS_sc),KB(nS_sc,nS_sc))
call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,eW, & call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) if(.not.TDA) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,OmRPA,rho_RPA,KA)
call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,OmRPA,rho_RPA,KB)
! Spin-conserved manifold ! Spin-conserved manifold
if(spin_conserved) then if(spin_conserved) then
ispin = 1 ispin = 1
allocate(Om_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) allocate(Om(nS_sc),XpY(nS_sc,nS_sc),XmY(nS_sc,nS_sc))
write(*,*) '------------------------' write(*,*) '------------------------'
write(*,*) 'Spin-conserved manifold ' write(*,*) 'Spin-conserved manifold '
@ -137,17 +144,27 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
if(doXBS) then if(doXBS) then
call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) if(.not.TDA) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KA)
call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KB)
end if end if
call phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,e, & call phULR_A(isp_W,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),Om_sc,XpY_sc,XmY_sc) if(.not.TDA) call phULR_B(isp_W,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
Aph(:,:) = Aph(:,:) + KA(:,:)
if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:)
call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcAC(ispin),OmRPA,XpY_RPA,XmY_RPA)
call phUACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc, & call phUACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc, &
ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sc,XmY_sc,Ec(iAC,ispin)) ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,XmY,Ec(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
@ -162,7 +179,7 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
write(*,*) '-----------------------------------------------------------------------------------' write(*,*) '-----------------------------------------------------------------------------------'
write(*,*) write(*,*)
deallocate(Om_sc,XpY_sc,XmY_sc) deallocate(Om,XpY,XmY)
end if end if
@ -174,7 +191,7 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
! Memory allocation ! Memory allocation
allocate(Om_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) allocate(Om(nS_sf),XpY(nS_sf,nS_sf),XmY(nS_sf,nS_sf))
write(*,*) '--------------------' write(*,*) '--------------------'
write(*,*) ' Spin-flip manifold ' write(*,*) ' Spin-flip manifold '
@ -191,17 +208,27 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
if(doXBS) then if(doXBS) then
call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) if(.not.TDA) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KA)
call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KB)
end if end if
call phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,lambda,e, & call phULR_A(isp_W,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),Om_sf,XpY_sf,XmY_sf) if(.not.TDA) call phULR_B(isp_W,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
Aph(:,:) = Aph(:,:) + KA(:,:)
if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:)
call phULR(TDA_W,nS_aa,nS_bb,nS_sf,Aph,Bph,EcAC(ispin),OmRPA,XpY_RPA,XmY_RPA)
call phUACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf, & call phUACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf, &
ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sf,XmY_sf,Ec(iAC,ispin)) ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,XmY,Ec(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
@ -216,7 +243,7 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
write(*,*) '-----------------------------------------------------------------------------------' write(*,*) '-----------------------------------------------------------------------------------'
write(*,*) write(*,*)
deallocate(Om_sf,XpY_sf,XmY_sf) deallocate(Om,XpY,XmY)
end if end if

96
src/GW/UGW_phBSE.f90
View File

@ -37,22 +37,25 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli
integer :: ispin integer :: ispin
integer :: isp_W integer :: isp_W
logical :: dRPA = .false.
logical :: dRPA_W = .true.
double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:)
double precision,allocatable :: KA(:,:)
double precision,allocatable :: KB(:,:)
integer :: nS_aa,nS_bb,nS_sc
double precision :: EcRPA double precision :: EcRPA
double precision,allocatable :: OmRPA(:) double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:) double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:) double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:,:) double precision,allocatable :: rho_RPA(:,:,:,:)
double precision,allocatable :: OmBSE_sc(:) integer :: nS_aa,nS_bb,nS_sc
double precision,allocatable :: XpY_BSE_sc(:,:)
double precision,allocatable :: XmY_BSE_sc(:,:)
integer :: nS_ab,nS_ba,nS_sf integer :: nS_ab,nS_ba,nS_sf
double precision,allocatable :: OmBSE_sf(:) double precision,allocatable :: OmBSE(:)
double precision,allocatable :: XpY_BSE_sf(:,:) double precision,allocatable :: XpY_BSE(:,:)
double precision,allocatable :: XmY_BSE_sf(:,:) double precision,allocatable :: XmY_BSE(:,:)
! Output variables ! Output variables
@ -68,31 +71,35 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli
nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1)) nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1))
nS_sf = nS_ab + nS_ba nS_sf = nS_ab + nS_ba
allocate(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin)) ! TDA
!-----!
! TDA !
!-----!
if(TDA) then if(TDA) then
write(*,*) 'Tamm-Dancoff approximation activated in phBSE!' write(*,*) 'Tamm-Dancoff approximation activated in phBSE!'
write(*,*) write(*,*)
end if end if
! Initialization
EcBSE(:) = 0d0
!--------------------------! !--------------------------!
! Spin-conserved screening ! ! Spin-conserved screening !
!--------------------------! !--------------------------!
isp_W = 1 isp_W = 1
EcRPA = 0d0
! Compute spin-conserved RPA screening ! Compute spin-conserved RPA screening
call phULR(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & allocate(Aph(nS_sc,nS_sc),Bph(nS_sc,nS_sc))
eW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) allocate(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin))
call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
if(.not.TDA_W) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
deallocate(Aph,Bph)
!----------------------------! !----------------------------!
! Spin-conserved excitations ! ! Spin-conserved excitations !
@ -101,17 +108,26 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli
if(spin_conserved) then if(spin_conserved) then
ispin = 1 ispin = 1
EcBSE(ispin) = 0d0
allocate(OmBSE_sc(nS_sc),XpY_BSE_sc(nS_sc,nS_sc),XmY_BSE_sc(nS_sc,nS_sc)) allocate(Aph(nS_sc,nS_sc),Bph(nS_sc,nS_sc),KA(nS_sc,nS_sc),KB(nS_sc,nS_sc))
allocate(OmBSE(nS_sc),XpY_BSE(nS_sc,nS_sc),XmY_BSE(nS_sc,nS_sc))
! Compute spin-conserved BSE excitation energies ! Compute spin-conserved BSE excitation energies
call phULR(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcBSE(ispin),OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
call print_excitation_energies('phBSE@GW@UHF','spin-conserved',nS_sc,OmBSE_sc)
call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,OmRPA,rho_RPA,KA)
if(.not.TDA) call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,OmRPA,rho_RPA,KB)
Aph(:,:) = Aph(:,:) + KA(:,:)
! if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:)
call phULR(TDA,nS_aa,nS_bb,nS_sc,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
call print_excitation_energies('phBSE@GW@UHF','spin-conserved',nS_sc,OmBSE)
call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
cW,S,OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) cW,S,OmBSE,XpY_BSE,XmY_BSE)
!------------------------------------------------- !-------------------------------------------------
! Compute the dynamical screening at the BSE level ! Compute the dynamical screening at the BSE level
@ -120,9 +136,10 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli
if(dBSE) & if(dBSE) &
call UGW_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,nS_sc, & call UGW_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,nS_sc, &
eW,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, & eW,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, &
OmRPA,rho_RPA,OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) OmRPA,rho_RPA,OmBSE,XpY_BSE,XmY_BSE)
deallocate(OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) deallocate(Aph,Bph,KA,KB)
deallocate(OmBSE,XpY_BSE,XmY_BSE)
end if end if
@ -133,21 +150,25 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli
if(spin_flip) then if(spin_flip) then
ispin = 2 ispin = 2
EcBSE(ispin) = 0d0
! Memory allocation allocate(Aph(nS_sf,nS_sf),Bph(nS_sf,nS_sf),KA(nS_sf,nS_sf),KB(nS_sf,nS_sf))
allocate(OmBSE(nS_sf),XpY_BSE(nS_sf,nS_sf),XmY_BSE(nS_sf,nS_sf))
allocate(OmBSE_sf(nS_sf),XpY_BSE_sf(nS_sf,nS_sf),XmY_BSE_sf(nS_sf,nS_sf)) ! Compute spin-conserved BSE excitation energies
call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,1d0,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
! Compute spin-flip BSE excitation energies call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sf,1d0,OmRPA,rho_RPA,KA)
if(.not.TDA) call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sf,1d0,OmRPA,rho_RPA,KB)
call phULR(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,1d0, & Aph(:,:) = Aph(:,:) + KA(:,:)
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcBSE(ispin), & if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:)
OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf)
call print_excitation_energies('phBSE@GW@UHF','spin-flip',nS_sf,OmBSE_sf) call phULR(TDA,nS_aa,nS_bb,nS_sc,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
call print_excitation_energies('phBSE@GW@UHF','spin-flip',nS_sf,OmBSE)
call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, &
cW,S,OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf) cW,S,OmBSE,XpY_BSE,XmY_BSE)
!------------------------------------------------- !-------------------------------------------------
! Compute the dynamical screening at the BSE level ! Compute the dynamical screening at the BSE level
@ -156,19 +177,18 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli
if(dBSE) & if(dBSE) &
call UGW_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,nS_sc, & call UGW_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,nS_sc, &
eW,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, & eW,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, &
OmRPA,rho_RPA,OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf) OmRPA,rho_RPA,OmBSE,XpY_BSE,XmY_BSE)
deallocate(OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf) deallocate(Aph,Bph,KA,KB)
deallocate(OmBSE,XpY_BSE,XmY_BSE)
end if end if
! Scale properly correlation energy if exchange is included in interaction kernel ! Scale properly correlation energy if exchange is included in interaction kernel
if(exchange_kernel) then if(exchange_kernel) then
EcBSE(1) = 0.5d0*EcBSE(1) EcBSE(:) = 0.5d0*EcBSE(:)
EcBSE(2) = 0.5d0*EcBSE(2)
else else

35
src/GW/UGW_phBSE_static_kernel_A.f90
View File

@ -1,5 +1,4 @@
subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,eGW, & subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,Om,rho,KA)
ERI_aaaa,ERI_aabb,ERI_bbbb,Omega,rho,A_lr)
! Compute the extra term for Bethe-Salpeter equation for linear response in the unrestricted formalism ! Compute the extra term for Bethe-Salpeter equation for linear response in the unrestricted formalism
@ -20,11 +19,7 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
integer,intent(in) :: nS_sc integer,intent(in) :: nS_sc
double precision,intent(in) :: eta double precision,intent(in) :: eta
double precision,intent(in) :: lambda double precision,intent(in) :: lambda
double precision,intent(in) :: eGW(nBas,nspin) double precision,intent(in) :: Om(nS_sc)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Omega(nS_sc)
double precision,intent(in) :: rho(nBas,nBas,nS_sc,nspin) double precision,intent(in) :: rho(nBas,nBas,nS_sc,nspin)
! Local variables ! Local variables
@ -35,7 +30,7 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
! Output variables ! Output variables
double precision,intent(out) :: A_lr(nSt,nSt) double precision,intent(out) :: KA(nSt,nSt)
!--------------------------------------------------! !--------------------------------------------------!
! Build BSE matrix for spin-conserving transitions ! ! Build BSE matrix for spin-conserving transitions !
@ -56,11 +51,11 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
chi = 0d0 chi = 0d0
do kc=1,nS_sc do kc=1,nS_sc
eps = Omega(kc)**2 + eta**2 eps = Om(kc)**2 + eta**2
chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Om(kc)/eps
end do end do
A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_aaaa(i,b,j,a) + 2d0*lambda*chi KA(ia,jb) = KA(ia,jb) + 2d0*lambda*chi
end do end do
end do end do
@ -80,11 +75,11 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
chi = 0d0 chi = 0d0
do kc=1,nS_sc do kc=1,nS_sc
eps = Omega(kc)**2 + eta**2 eps = Om(kc)**2 + eta**2
chi = chi + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps chi = chi + rho(i,j,kc,2)*rho(a,b,kc,2)*Om(kc)/eps
end do end do
A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,b,j,a) + 2d0*lambda*chi KA(nSa+ia,nSa+jb) = KA(nSa+ia,nSa+jb) + 2d0*lambda*chi
end do end do
end do end do
@ -112,11 +107,11 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
chi = 0d0 chi = 0d0
do kc=1,nS_sc do kc=1,nS_sc
eps = Omega(kc)**2 + eta**2 eps = Om(kc)**2 + eta**2
chi = chi + rho(i,j,kc,1)*rho(a,b,kc,2)*Omega(kc)/eps chi = chi + rho(i,j,kc,1)*rho(a,b,kc,2)*Om(kc)/eps
end do end do
A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_aabb(i,b,j,a) + 2d0*lambda*chi KA(ia,jb) = KA(ia,jb) + 2d0*lambda*chi
end do end do
end do end do
@ -136,11 +131,11 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
chi = 0d0 chi = 0d0
do kc=1,nS_sc do kc=1,nS_sc
eps = Omega(kc)**2 + eta**2 eps = Om(kc)**2 + eta**2
chi = chi + rho(i,j,kc,2)*rho(a,b,kc,1)*Omega(kc)/eps chi = chi + rho(i,j,kc,2)*rho(a,b,kc,1)*Om(kc)/eps
end do end do
A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_aabb(b,i,a,j) + 2d0*lambda*chi KA(nSa+ia,nSa+jb) = KA(nSa+ia,nSa+jb) + 2d0*lambda*chi
end do end do
end do end do

34
src/GW/UGW_phBSE_static_kernel_B.f90
View File

@ -1,5 +1,4 @@
subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,Om,rho,KB)
ERI_aaaa,ERI_aabb,ERI_bbbb,Omega,rho,B_lr)
! Compute the extra term for Bethe-Salpeter equation for linear response ! Compute the extra term for Bethe-Salpeter equation for linear response
@ -20,10 +19,7 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
integer,intent(in) :: nS_sc integer,intent(in) :: nS_sc
double precision,intent(in) :: eta double precision,intent(in) :: eta
double precision,intent(in) :: lambda double precision,intent(in) :: lambda
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(in) :: Om(nS_sc)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Omega(nS_sc)
double precision,intent(in) :: rho(nBas,nBas,nS_sc,nspin) double precision,intent(in) :: rho(nBas,nBas,nS_sc,nspin)
! Local variables ! Local variables
@ -34,7 +30,7 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
! Output variables ! Output variables
double precision,intent(out) :: B_lr(nSt,nSt) double precision,intent(out) :: KB(nSt,nSt)
!--------------------------------------------------! !--------------------------------------------------!
! Build BSE matrix for spin-conserving transitions ! ! Build BSE matrix for spin-conserving transitions !
@ -55,11 +51,11 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
chi = 0d0 chi = 0d0
do kc=1,nS_sc do kc=1,nS_sc
eps = Omega(kc)**2 + eta**2 eps = Om(kc)**2 + eta**2
chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Om(kc)/eps
end do end do
B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI_aaaa(i,j,b,a) + 2d0*lambda*chi KB(ia,jb) = KB(ia,jb) + 2d0*lambda*chi
end do end do
end do end do
@ -80,11 +76,11 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
chi = 0d0 chi = 0d0
do kc=1,nS_sc do kc=1,nS_sc
eps = Omega(kc)**2 + eta**2 eps = Om(kc)**2 + eta**2
chi = chi + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps chi = chi + rho(i,b,kc,2)*rho(a,j,kc,2)*Om(kc)/eps
end do end do
B_lr(nSa+ia,nSa+jb) = B_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,j,b,a) + 2d0*lambda*chi KB(nSa+ia,nSa+jb) = KB(nSa+ia,nSa+jb) + 2d0*lambda*chi
end do end do
end do end do
@ -113,11 +109,11 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
chi = 0d0 chi = 0d0
do kc=1,nS_sc do kc=1,nS_sc
eps = Omega(kc)**2 + eta**2 eps = Om(kc)**2 + eta**2
chi = chi + rho(i,b,kc,1)*rho(a,j,kc,2)*Omega(kc)/eps chi = chi + rho(i,b,kc,1)*rho(a,j,kc,2)*Om(kc)/eps
end do end do
B_lr(ia,nSa+jb) = B_lr(ia,nSa+jb) - lambda*ERI_aabb(i,j,b,a) + 2d0*lambda*chi KB(ia,nSa+jb) = KB(ia,nSa+jb) + 2d0*lambda*chi
end do end do
end do end do
@ -137,11 +133,11 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s
chi = 0d0 chi = 0d0
do kc=1,nS_sc do kc=1,nS_sc
eps = Omega(kc)**2 + eta**2 eps = Om(kc)**2 + eta**2
chi = chi + rho(i,b,kc,2)*rho(a,j,kc,1)*Omega(kc)/eps chi = chi + rho(i,b,kc,2)*rho(a,j,kc,1)*Om(kc)/eps
end do end do
B_lr(nSa+ia,jb) = B_lr(nSa+ia,jb) - lambda*ERI_aabb(j,i,a,b) + 2d0*lambda*chi KB(nSa+ia,jb) = KB(nSa+ia,jb) + 2d0*lambda*chi
end do end do
end do end do

10
src/LR/phULR.f90
View File

@ -34,12 +34,6 @@ subroutine phULR(TDA,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY)
allocate(ApB(nSt,nSt),AmB(nSt,nSt),AmBSq(nSt,nSt),AmBIv(nSt,nSt),Z(nSt,nSt)) allocate(ApB(nSt,nSt),AmB(nSt,nSt),AmBSq(nSt,nSt),AmBIv(nSt,nSt),Z(nSt,nSt))
! Build A and B matrices
! if(BSE) &
! call UGW_phBSE_static_kernel_A(ispin,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, &
! ERI_aaaa,ERI_aabb,ERI_bbbb,Om,rho,Aph)
! Tamm-Dancoff approximation ! Tamm-Dancoff approximation
if(TDA) then if(TDA) then
@ -51,10 +45,6 @@ subroutine phULR(TDA,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY)
else else
! if(BSE) &
! call UGW_phBSE_static_kernel_B(ispin,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, &
! ERI_aaaa,ERI_aabb,ERI_bbbb,Om,rho,Bph)
ApB(:,:) = Aph(:,:) + Bph(:,:) ApB(:,:) = Aph(:,:) + Bph(:,:)
AmB(:,:) = Aph(:,:) - Bph(:,:) AmB(:,:) = Aph(:,:) - Bph(:,:)

14
src/LR/phULR_transition_vectors.f90
View File

@ -1,4 +1,4 @@
subroutine phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb,c,S,Omega,XpY,XmY) subroutine phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY)
! Print transition vectors for linear response calculation ! Print transition vectors for linear response calculation
@ -21,14 +21,14 @@ subroutine phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole
double precision :: dipole_int_bb(nBas,nBas,ncart) double precision :: dipole_int_bb(nBas,nBas,ncart)
double precision,intent(in) :: c(nBas,nBas,nspin) double precision,intent(in) :: c(nBas,nBas,nspin)
double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: Omega(nSt) double precision,intent(in) :: Om(nSt)
double precision,intent(in) :: XpY(nSt,nSt) double precision,intent(in) :: XpY(nSt,nSt)
double precision,intent(in) :: XmY(nSt,nSt) double precision,intent(in) :: XmY(nSt,nSt)
! Local variables ! Local variables
integer :: ia,jb,j,b integer :: ia,jb,j,b
integer :: maxS = 20 integer :: maxS = 10
double precision,parameter :: thres_vec = 0.1d0 double precision,parameter :: thres_vec = 0.1d0
double precision,allocatable :: X(:) double precision,allocatable :: X(:)
double precision,allocatable :: Y(:) double precision,allocatable :: Y(:)
@ -44,11 +44,11 @@ subroutine phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole
os(:) = 0d0 os(:) = 0d0
if(ispin == 1) call phULR_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS, & if(ispin == 1) call phULR_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS, &
dipole_int_aa,dipole_int_bb,Omega,XpY,XmY,os) dipole_int_aa,dipole_int_bb,Om,XpY,XmY,os)
! Compute <S**2> ! Compute <S**2>
call S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,Omega,XpY,XmY,S2) call S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,Om,XpY,XmY,S2)
! Print details about spin-conserved excitations ! Print details about spin-conserved excitations
@ -61,7 +61,7 @@ subroutine phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole
print*,'-------------------------------------------------------------' print*,'-------------------------------------------------------------'
write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') &
' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' f = ',os(ia),' <S**2> = ',S2(ia) ' Excitation n. ',ia,': ',Om(ia)*HaToeV,' eV',' f = ',os(ia),' <S**2> = ',S2(ia)
print*,'-------------------------------------------------------------' print*,'-------------------------------------------------------------'
! Spin-up transitions ! Spin-up transitions
@ -117,7 +117,7 @@ subroutine phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole
print*,'-------------------------------------------------------------' print*,'-------------------------------------------------------------'
write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') &
' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' f = ',os(ia),' <S**2> = ',S2(ia) ' Excitation n. ',ia,': ',Om(ia)*HaToeV,' eV',' f = ',os(ia),' <S**2> = ',S2(ia)
print*,'-------------------------------------------------------------' print*,'-------------------------------------------------------------'
! Spin-up transitions ! Spin-up transitions

2
src/LR/print_excitation_energies.f90
View File

@ -14,7 +14,7 @@ subroutine print_excitation_energies(method,manifold,nS,Om)
! Local variables ! Local variables
integer,parameter :: maxS = 20 integer,parameter :: maxS = 10
integer :: m integer :: m
write(*,*) write(*,*)