LCPQ
/
quack
mirror of https://github.com/pfloos/quack
4
1
Fork 0
Code Releases Activity

OK with LR in RPA routines

This commit is contained in:
Pierre-Francois Loos 2023-07-18 11:53:38 +02:00
parent 266be65f2e
commit 092a3f5e6d
13 changed files with 159 additions and 238 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
input/methods
View File

@ -9,7 +9,7 @@
# CIS* CIS(D) CID CISD FCI
F F F F F
# phRPA* phRPAx* crRPA ppRPA
F F F T
F F T F
# G0F2* evGF2* qsGF2* G0F3 evGF3
F F F F F
# G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW

2
input/options
View File

@ -13,6 +13,6 @@
# GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg
256 0.00001 T 5 T 0.0 F F
# ACFDT: AC Kx XBS
T F F
T T F
# BSE: phBSE phBSE2 ppBSE dBSE dTDA evDyn
F F F T F F

2
src/GW/G0W0.f90
View File

@ -229,7 +229,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD
end if
call phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF,eGW,EcAC)
call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF,eGW,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

2
src/GW/SRG_qsGW.f90
View File

@ -342,7 +342,7 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE
end if
call phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC)
call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

2
src/GW/evGW.f90
View File

@ -290,7 +290,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
end if
call phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC)
call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

2
src/GW/qsGW.f90
View File

@ -331,7 +331,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
end if
call phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC)
call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

60
src/LR/phLR.f90
View File

@ -1,4 +1,4 @@
subroutine phLR(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Ec,Omega,XpY,XmY)
subroutine phLR(TDA,nS,A,B,EcRPA,Om,XpY,XmY)
! Compute linear response
@ -7,27 +7,16 @@ subroutine phLR(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Ec,Omega,XpY
! Input variables
logical,intent(in) :: dRPA
logical,intent(in) :: TDA
double precision,intent(in) :: eta
integer,intent(in) :: ispin
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(nS,nS)
double precision,intent(in) :: B(nS,nS)
! Local variables
integer :: i,j,k
double precision :: trace_matrix
double precision,allocatable :: A(:,:)
double precision,allocatable :: B(:,:)
double precision,allocatable :: ApB(:,:)
double precision,allocatable :: AmB(:,:)
double precision,allocatable :: AmBSq(:,:)
@ -37,33 +26,26 @@ subroutine phLR(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Ec,Omega,XpY
! Output variables
double precision,intent(out) :: Ec
double precision,intent(out) :: Omega(nS)
double precision,intent(out) :: EcRPA
double precision,intent(out) :: Om(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),tmp(nS,nS))
! Build A and B matrices
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A)
allocate(ApB(nS,nS),AmB(nS,nS),AmBSq(nS,nS),AmBIv(nS,nS),Z(nS,nS),tmp(nS,nS))
! Tamm-Dancoff approximation
if(TDA) then
B(:,:) = 0d0
XpY(:,:) = A(:,:)
call diagonalize_matrix(nS,XpY,Omega)
call diagonalize_matrix(nS,XpY,Om)
XpY(:,:) = transpose(XpY(:,:))
XmY(:,:) = XpY(:,:)
else
call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B)
! Build A + B and A - B matrices
ApB = A + B
@ -71,42 +53,34 @@ subroutine phLR(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Ec,Omega,XpY
! Diagonalize linear response matrix
call diagonalize_matrix(nS,AmB,Omega)
call diagonalize_matrix(nS,AmB,Om)
if(minval(Omega) < 0d0) &
if(minval(Om) < 0d0) &
call print_warning('You may have instabilities in linear response: A-B is not positive definite!!')
! do ia=1,nS
! if(Omega(ia) < 0d0) Omega(ia) = 0d0
! end do
call ADAt(nS,AmB,1d0*dsqrt(Omega),AmBSq)
call ADAt(nS,AmB,1d0/dsqrt(Omega),AmBIv)
call ADAt(nS,AmB,1d0*dsqrt(Om),AmBSq)
call ADAt(nS,AmB,1d0/dsqrt(Om),AmBIv)
call dgemm('N','N',nS,nS,nS,1d0,ApB,size(ApB,1),AmBSq,size(AmBSq,1),0d0,tmp,size(tmp,1))
call dgemm('N','N',nS,nS,nS,1d0,AmBSq,size(AmBSq,1),tmp,size(tmp,1),0d0,Z,size(Z,1))
call diagonalize_matrix(nS,Z,Omega)
call diagonalize_matrix(nS,Z,Om)
if(minval(Omega) < 0d0) &
if(minval(Om) < 0d0) &
call print_warning('You may have instabilities in linear response: negative excitations!!')
! do ia=1,nS
! if(Omega(ia) < 0d0) Omega(ia) = 0d0
! end do
Omega = dsqrt(Omega)
Om = sqrt(Om)
call dgemm('T','N',nS,nS,nS,1d0,Z,size(Z,1),AmBSq,size(AmBSq,1),0d0,XpY,size(XpY,1))
call DA(nS,1d0/dsqrt(Omega),XpY)
call DA(nS,1d0/dsqrt(Om),XpY)
call dgemm('T','N',nS,nS,nS,1d0,Z,size(Z,1),AmBIv,size(AmBIv,1),0d0,XmY,size(XmY,1))
call DA(nS,1d0*dsqrt(Omega),XmY)
call DA(nS,1d0*dsqrt(Om),XmY)
end if
! Compute the RPA correlation energy
Ec = 0.5d0*(sum(Omega) - trace_matrix(nS,A))
EcRPA = 0.5d0*(sum(Om) - trace_matrix(nS,A))
end subroutine

6
src/QuAcK/QuAcK.f90
View File

@ -581,7 +581,7 @@ program QuAcK
else
call phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF)
call phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF)
end if
call cpu_time(end_RPA)
@ -606,7 +606,7 @@ program QuAcK
else
call phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF)
call phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF)
end if
call cpu_time(end_RPA)
@ -624,7 +624,7 @@ program QuAcK
if(docrRPA) then
call cpu_time(start_RPA)
call crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF)
call crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF)
call cpu_time(end_RPA)
t_RPA = end_RPA - start_RPA

57
src/RPA/crACFDT.f90
View File

@ -1,4 +1,4 @@
subroutine crACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,e,EcAC)
subroutine crACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem
! for the crossed-ring contribution
@ -9,41 +9,31 @@ subroutine crACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,
! Input variables
logical,intent(in) :: doXBS
logical,intent(in) :: exchange_kernel
logical,intent(in) :: dRPA
logical,intent(in) :: TDA_W
logical,intent(in) :: TDA
logical,intent(in) :: BSE
logical,intent(in) :: singlet
logical,intent(in) :: triplet
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) :: eW(nBas)
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: ispin
integer :: isp_W
integer :: iAC
double precision :: lambda
double precision,allocatable :: Ec(:,:)
double precision :: EcRPA
double precision,allocatable :: KA(:,:)
double precision,allocatable :: KB(:,:)
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:)
double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:)
double precision,allocatable :: Om(:)
double precision,allocatable :: XpY(:,:)
@ -56,8 +46,7 @@ subroutine crACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,
! Memory allocation
allocate(Ec(nAC,nspin))
allocate(KA(nS,nS),KB(nS,nS),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS))
allocate(Om(nS),XpY(nS,nS),XmY(nS,nS))
allocate(Aph(nS,nS),Bph(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS))
! Antisymmetrized kernel version
@ -72,17 +61,6 @@ subroutine crACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,
EcAC(:) = 0d0
Ec(:,:) = 0d0
! Compute (singlet) RPA screening
isp_W = 1
EcRPA = 0d0
call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA)
call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB)
! Singlet manifold
if(singlet) then
@ -102,18 +80,10 @@ subroutine crACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,
lambda = -rAC(iAC)
if(doXBS) then
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
! call print_excitation('W^lambda: ',isp_W,nS,OmRPA)
call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA)
call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB)
end if
call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,KA,KB,EcAC(ispin),Om,XpY,XmY)
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))
@ -151,17 +121,10 @@ subroutine crACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,
lambda = -rAC(iAC)
if(doXBS) then
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA)
call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB)
end if
call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,KA,KB,EcAC(ispin),Om,XpY,XmY)
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))

65
src/RPA/crRPA.f90
View File

@ -1,5 +1,4 @@
subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, &
ERI,dipole_int,eHF)
subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
! Crossed-ring channel of the random phase approximation
@ -13,7 +12,6 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
logical,intent(in) :: singlet
double precision,intent(in) :: eta
logical,intent(in) :: triplet
integer,intent(in) :: nBas
integer,intent(in) :: nC
@ -22,19 +20,22 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: EHF
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
integer :: ispin
double precision,allocatable :: Om(:,:)
double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:)
logical :: dRPA
double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:)
double precision,allocatable :: Om(:)
double precision,allocatable :: XpY(:,:)
double precision,allocatable :: XmY(:,:)
double precision :: EcRPAx(nspin)
double precision :: EcTr(nspin)
double precision :: EcAC(nspin)
! Hello world
@ -54,12 +55,15 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
! Initialization
EcRPAx(:) = 0d0
dRPA = .false.
EcTr(:) = 0d0
EcAC(:) = 0d0
! Memory allocation
allocate(Om(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS))
if(.not.TDA) allocate(Bph(nS,nS))
! Singlet manifold
@ -67,9 +71,12 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
ispin = 1
call phLR(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,EcRPAx(ispin),Om(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('crRPA@HF ',ispin,nS,Om(:,ispin))
call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,e,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,ERI,Bph)
call phLR(TDA,nS,Aph,Bph,EcTr(ispin),Om,XpY,XmY)
call print_excitation('crRPA@HF ',ispin,nS,Om)
call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
endif
@ -79,25 +86,28 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
ispin = 2
call phLR(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,EcRPAx(ispin),Om(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('crRPA@HF ',ispin,nS,Om(:,ispin))
call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,e,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,ERI,Bph)
call phLR(TDA,nS,Aph,Bph,EcTr(ispin),Om,XpY,XmY)
call print_excitation('crRPA@HF ',ispin,nS,Om)
call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
endif
! if(exchange_kernel) then
if(exchange_kernel) then
EcRPAx(1) = 0.5d0*EcRPAx(1)
EcRPAx(2) = 1.5d0*EcRPAx(2)
EcTr(1) = 0.5d0*EcTr(1)
EcTr(2) = 1.5d0*EcTr(2)
! end if
end if
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (singlet) =',EcRPAx(1)
write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (triplet) =',EcRPAx(2)
write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy =',EcRPAx(1) + EcRPAx(2)
write(*,'(2X,A50,F20.10)') 'Tr@crRPA total energy =',ENuc + ERHF + EcRPAx(1) + EcRPAx(2)
write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (singlet) =',EcTr(1)
write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (triplet) =',EcTr(2)
write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy =',EcTr(1) + EcTr(2)
write(*,'(2X,A50,F20.10)') 'Tr@crRPA total energy =',ENuc + EHF + EcTr(1) + EcTr(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -110,15 +120,14 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
write(*,*) '-------------------------------------------------------'
write(*,*)
call crACFDT(exchange_kernel,.false.,.false.,.false.,TDA,.false.,singlet,triplet,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC)
call crACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@crRPA total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@crRPA total energy =',ENuc + EHF + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)

56
src/RPA/phACFDT.f90
View File

@ -1,4 +1,4 @@
subroutine phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,e,EcAC)
subroutine phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem
@ -8,41 +8,31 @@ subroutine phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,
! Input variables
logical,intent(in) :: doXBS
logical,intent(in) :: exchange_kernel
logical,intent(in) :: dRPA
logical,intent(in) :: TDA_W
logical,intent(in) :: TDA
logical,intent(in) :: BSE
logical,intent(in) :: singlet
logical,intent(in) :: triplet
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) :: eW(nBas)
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: ispin
integer :: isp_W
integer :: iAC
double precision :: lambda
double precision,allocatable :: Ec(:,:)
double precision :: EcRPA
double precision,allocatable :: KA(:,:)
double precision,allocatable :: KB(:,:)
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:)
double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:)
double precision,allocatable :: Om(:)
double precision,allocatable :: XpY(:,:)
@ -55,8 +45,7 @@ subroutine phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,
! Memory allocation
allocate(Ec(nAC,nspin))
allocate(KA(nS,nS),KB(nS,nS),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS))
allocate(Om(nS),XpY(nS,nS),XmY(nS,nS))
allocate(Aph(nS,nS),Bph(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS))
! Antisymmetrized kernel version
@ -71,17 +60,6 @@ subroutine phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,
EcAC(:) = 0d0
Ec(:,:) = 0d0
! Compute (singlet) RPA screening
isp_W = 1
EcRPA = 0d0
call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA)
call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB)
! Singlet manifold
if(singlet) then
@ -101,17 +79,10 @@ subroutine phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,
lambda = rAC(iAC)
if(doXBS) then
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA)
call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB)
end if
call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,KA,KB,EcAC(ispin),Om,XpY,XmY)
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))
@ -149,17 +120,10 @@ subroutine phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,
lambda = rAC(iAC)
if(doXBS) then
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA)
call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB)
end if
call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,KA,KB,EcAC(ispin),Om,XpY,XmY)
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))

71
src/RPA/phRPA.f90
View File

@ -1,4 +1,4 @@
subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
! Perform a direct random phase approximation calculation
@ -13,7 +13,6 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
logical,intent(in) :: exchange_kernel
logical,intent(in) :: singlet
logical,intent(in) :: triplet
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
@ -21,19 +20,22 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: EHF
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
integer :: ispin
logical :: dRPA
double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:)
double precision,allocatable :: Om(:)
double precision,allocatable :: XpY(:,:)
double precision,allocatable :: XmY(:,:)
double precision :: EcRPA(nspin)
double precision :: EcTr(nspin)
double precision :: EcAC(nspin)
! Hello world
@ -53,12 +55,15 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
! Initialization
EcRPA(:) = 0d0
dRPA = .true.
EcTr(:) = 0d0
EcAC(:) = 0d0
! Memory allocation
allocate(Om(nS),XpY(nS,nS),XmY(nS,nS))
allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS))
if(.not.TDA) allocate(Bph(nS,nS))
! Singlet manifold
@ -66,8 +71,11 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
ispin = 1
call phLR(ispin,.true.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,EcRPA(ispin),Om,XpY,XmY)
call print_excitation('RPA@HF ',ispin,nS,Om)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA,nS,Aph,Bph,EcTr(ispin),Om,XpY,XmY)
call print_excitation('phRPA@HF ',ispin,nS,Om)
call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
endif
@ -78,60 +86,53 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
ispin = 2
call phLR(ispin,.true.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,EcRPA(ispin),Om,XpY,XmY)
call print_excitation('RPA@HF ',ispin,nS,Om)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA,nS,Aph,Bph,EcTr(ispin),Om,XpY,XmY)
call print_excitation('phRPA@HF ',ispin,nS,Om)
call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
endif
! if(exchange_kernel) then
if(exchange_kernel) then
! EcRPA(1) = 0.5d0*EcRPA(1)
! EcRPA(2) = 1.5d0*EcRPA(2)
EcTr(1) = 0.5d0*EcTr(1)
EcTr(2) = 1.5d0*EcTr(2)
! end if
end if
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@phRPA total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy (singlet) =',EcTr(1)
write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy (triplet) =',EcTr(2)
write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy =',EcTr(1) + EcTr(2)
write(*,'(2X,A50,F20.10)') 'Tr@phRPA total energy =',ENuc + EHF + EcTr(1) + EcTr(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
deallocate(Om,XpY,XmY)
deallocate(Om,XpY,XmY,Aph,Bph)
! Compute the correlation energy via the adiabatic connection
! Switch off ACFDT for RPA as the trace formula is equivalent
if(doACFDT) then
write(*,*) '------------------------------------------------------'
write(*,*) 'Adiabatic connection version of RPA correlation energy'
write(*,*) '------------------------------------------------------'
write(*,*) '--------------------------------------------------------'
write(*,*) 'Adiabatic connection version of phRPA correlation energy'
write(*,*) '--------------------------------------------------------'
write(*,*)
call phACFDT(exchange_kernel,.false.,.true.,.false.,TDA,.false.,singlet,triplet,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC)
if(exchange_kernel) then
EcAC(1) = 0.5d0*EcAC(1)
EcAC(2) = 1.5d0*EcAC(2)
end if
call phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'AC@phRPA correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A50,F20.10)') 'AC@phRPA correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@phRPA correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@phRPA total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@phRPA total energy =',ENuc + EHF + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end if
end subroutine

70
src/RPA/phRPAx.f90
View File

@ -1,4 +1,4 @@
subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
! Perform random phase approximation calculation with exchange (aka TDHF)
@ -13,7 +13,6 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,
logical,intent(in) :: exchange_kernel
logical,intent(in) :: singlet
logical,intent(in) :: triplet
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
@ -21,19 +20,22 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: EHF
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
integer :: ispin
logical :: dRPA
double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:)
double precision,allocatable :: Om(:)
double precision,allocatable :: XpY(:,:)
double precision,allocatable :: XmY(:,:)
double precision :: EcRPAx(nspin)
double precision :: EcTr(nspin)
double precision :: EcAC(nspin)
! Hello world
@ -54,12 +56,15 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,
! Initialization
EcRPAx(:) = 0d0
EcAC(:) = 0d0
dRPA = .false.
EcTr(:) = 0d0
EcAC(:) = 0d0
! Memory allocation
allocate(Om(nS),XpY(nS,nS),XmY(nS,nS))
allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS))
if(.not.TDA) allocate(Bph(nS,nS))
! Singlet manifold
@ -67,8 +72,11 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,
ispin = 1
call phLR(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,EcRPAx(ispin),Om,XpY,XmY)
call print_excitation('RPAx@HF ',ispin,nS,Om)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA,nS,Aph,Bph,EcTr(ispin),Om,XpY,XmY)
call print_excitation('phRPAx@HF ',ispin,nS,Om)
call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
endif
@ -79,50 +87,52 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,
ispin = 2
call phLR(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,EcRPAx(ispin),Om,XpY,XmY)
call print_excitation('RPAx@HF ',ispin,nS,Om)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph)
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA,nS,Aph,Bph,EcTr(ispin),Om,XpY,XmY)
call print_excitation('phRPAx@HF ',ispin,nS,Om)
call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
endif
! if(exchange_kernel) then
if(exchange_kernel) then
EcRPAx(1) = 0.5d0*EcRPAx(1)
EcRPAx(2) = 1.5d0*EcRPAx(2)
EcTr(1) = 0.5d0*EcTr(1)
EcTr(2) = 1.5d0*EcTr(2)
! end if
end if
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@RPAx correlation energy (singlet) =',EcRPAx(1)
write(*,'(2X,A50,F20.10)') 'Tr@RPAx correlation energy (triplet) =',EcRPAx(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPAx correlation energy =',EcRPAx(1) + EcRPAx(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPAx total energy =',ENuc + ERHF + EcRPAx(1) + EcRPAx(2)
write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy (singlet) =',EcTr(1)
write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy (triplet) =',EcTr(2)
write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy =',EcTr(1) + EcTr(2)
write(*,'(2X,A50,F20.10)') 'Tr@phRPAx total energy =',ENuc + EHF + EcTr(1) + EcTr(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! deallocate memory
deallocate(Om,XpY,XmY)
deallocate(Om,XpY,XmY,Aph,Bph)
! Compute the correlation energy via the adiabatic connection
if(doACFDT) then
write(*,*) '-------------------------------------------------------'
write(*,*) 'Adiabatic connection version of RPAx correlation energy'
write(*,*) '-------------------------------------------------------'
write(*,*) '---------------------------------------------------------'
write(*,*) 'Adiabatic connection version of phRPAx correlation energy'
write(*,*) '---------------------------------------------------------'
write(*,*)
call phACFDT(exchange_kernel,.false.,.false.,.false.,TDA,.false.,singlet,triplet,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC)
call phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@RPAx total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@phRPAx correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A50,F20.10)') 'AC@phRPAx correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@phRPAx correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@phRPAx total energy =',ENuc + EHF + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)