4
1
mirror of https://github.com/pfloos/quack synced 2024-11-07 06:33:55 +01:00

crossed ring ACFDT

This commit is contained in:
Pierre-Francois Loos 2021-11-10 22:53:50 +01:00
parent 23d7d83091
commit 6f886a61a6
4 changed files with 177 additions and 4 deletions

View File

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

View File

@ -13,6 +13,6 @@
# ACFDT: AC Kx XBS # ACFDT: AC Kx XBS
T T T T T T
# BSE: BSE dBSE dTDA evDyn # BSE: BSE dBSE dTDA evDyn
T T T F T F T F
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T 1000000 100000 10 0.3 10000 1234 T

173
src/RPA/ACFDT_cr.f90 Normal file
View File

@ -0,0 +1,173 @@
subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,e,EcAC)
! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem
! for the crossed-ring contribution
implicit none
include 'parameters.h'
include 'quadrature.h'
! 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,nC,nO,nV,nR,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 :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:)
double precision,allocatable :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:)
! Output variables
double precision,intent(out) :: EcAC(nspin)
! Memory allocation
allocate(Ec(nAC,nspin))
allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS))
allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
! Antisymmetrized kernel version
if(exchange_kernel) then
write(*,*)
write(*,*) '*** Exchange kernel version ***'
write(*,*)
end if
EcAC(:) = 0d0
Ec(:,:) = 0d0
! Compute (singlet) RPA screening
isp_W = 1
EcRPA = 0d0
call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,OmRPA, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
! Singlet manifold
if(singlet) then
ispin = 1
write(*,*) '--------------'
write(*,*) 'Singlet states'
write(*,*) '--------------'
write(*,*)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)'
write(*,*) '-----------------------------------------------------------------------------------'
do iAC=1,nAC
lambda = -rAC(iAC)
if(doXBS) then
call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
! call print_excitation('W^lambda: ',isp_W,nS,OmRPA)
end if
call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, &
rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, &
ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
end do
EcAC(ispin) = -0.5d0*dot_product(wAC,Ec(:,ispin))
if(exchange_kernel) EcAC(ispin) = 0.5d0*EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,*)
end if
! Triplet manifold
if(triplet) then
ispin = 2
write(*,*) '--------------'
write(*,*) 'Triplet states'
write(*,*) '--------------'
write(*,*)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)'
write(*,*) '-----------------------------------------------------------------------------------'
do iAC=1,nAC
lambda = -rAC(iAC)
if(doXBS) then
call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
end if
call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, &
rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
end do
EcAC(ispin) = -0.5d0*dot_product(wAC,Ec(:,ispin))
if(exchange_kernel) EcAC(ispin) = 1.5d0*EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,*)
end if
end subroutine ACFDT_cr

View File

@ -113,8 +113,8 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
write(*,*) '-------------------------------------------------------' write(*,*) '-------------------------------------------------------'
write(*,*) write(*,*)
call ACFDT(exchange_kernel,.false.,.false.,.false.,TDA,.false.,singlet,triplet,eta, & call ACFDT_cr(exchange_kernel,.false.,.false.,.false.,TDA,.false.,singlet,triplet,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC) nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC)
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'