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

crGRPA

This commit is contained in:
Pierre-Francois Loos 2023-11-14 13:50:53 +01:00
parent 3318daf241
commit 0bb137ed57
14 changed files with 168 additions and 97 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
input/methods
View File

@ -1,7 +1,7 @@
# RHF UHF GHF ROHF
F F T F
# MP2 MP3
T T
F F
# CCD pCCD DCD CCSD CCSD(T)
F F F F F
# drCCD rCCD crCCD lCCD
@ -9,7 +9,7 @@
# CIS CIS(D) CID CISD FCI
F F F F F
# phRPA phRPAx crRPA ppRPA
F F F F
T T T F
# G0F2 evGF2 qsGF2 G0F3 evGF3
F F F F F
# G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW

2
input/options
View File

@ -1,5 +1,5 @@
# HF: maxSCF thresh DIIS guess mix shift stab search
10000 0.0000001 5 2 0.0 0.0 F T
10000 0.0000001 5 1 0.0 0.0 F F
# MP: reg
F
# CC: maxSCF thresh DIIS

6
src/GW/qsRGW.f90
View File

@ -349,9 +349,9 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
if(dotest) then
call dump_test_value('R','qsRGW correlation energy',EcRPA)
call dump_test_value('R','qsRGW HOMO energy',eGW(nO))
call dump_test_value('R','qsRGW LUMO energy',eGW(nO+1))
call dump_test_value('R','qsGW correlation energy',EcRPA)
call dump_test_value('R','qsGW HOMO energy',eGW(nO))
call dump_test_value('R','qsGW LUMO energy',eGW(nO+1))
end if

9
src/QuAcK/GQuAcK.f90
View File

@ -1,5 +1,5 @@
subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,doppRPA, &
dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, &
nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, &
maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, &
@ -20,7 +20,7 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
logical,intent(in) :: doMP3
logical,intent(in) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT
logical,intent(in) :: dodrCCD,dorCCD,docrCCD,dolCCD
logical,intent(in) :: dophRPA,dophRPAx,doppRPA
logical,intent(in) :: dophRPA,dophRPAx,docrRPA,doppRPA
logical,intent(in) :: doG0F2,doevGF2,doqsGF2
logical,intent(in) :: doG0W0,doevGW,doqsGW
@ -236,13 +236,12 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
! Random-phase approximation module !
!-----------------------------------!
doRPA = dophRPA .or. dophRPAx .or. doppRPA
doRPA = dophRPA .or. dophRPAx .or. docrRPA .or. doppRPA
if(doRPA) then
call wall_time(start_RPA)
call GRPA(dotest,dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas2,nC,nO,nV,nR,nS,ENuc,EGHF, &
ERI_MO,dipole_int_MO,eHF,cHF,S)
call GRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,nBas2,nC,nO,nV,nR,nS,ENuc,EGHF,ERI_MO,dipole_int_MO,eHF)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA

2
src/QuAcK/QuAcK.f90
View File

@ -227,7 +227,7 @@ program QuAcK
if(doGQuAcK) &
call GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,doppRPA, &
dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, &
nNuc,nBas,sum(nC),sum(nO),sum(nV),sum(nR),ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, &
maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, &

31
src/RPA/GRPA.f90
View File

@ -1,4 +1,4 @@
subroutine GRPA(dotest,dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
subroutine GRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
! Random-phase approximation module
@ -11,11 +11,10 @@ subroutine GRPA(dotest,dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas
logical,intent(in) :: dophRPA
logical,intent(in) :: dophRPAx
logical,intent(in) :: docrRPA
logical,intent(in) :: doppRPA
logical,intent(in) :: TDA
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
@ -23,8 +22,8 @@ subroutine GRPA(dotest,dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: EHF
double precision,intent(in) :: epsHF(nBas)
double precision,intent(in) :: EGHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
@ -39,7 +38,7 @@ subroutine GRPA(dotest,dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas
if(dophRPA) then
call wall_time(start_RPA)
call phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
@ -55,7 +54,7 @@ subroutine GRPA(dotest,dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas
if(dophRPAx) then
call wall_time(start_RPA)
call phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
call wall_time(end_RPA)
@ -65,6 +64,22 @@ subroutine GRPA(dotest,dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas
end if
!------------------------------------------------------------------------
! Compute crRPA excitations
!------------------------------------------------------------------------
if(docrRPA) then
call wall_time(start_RPA)
call crGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for cr-RPA = ',t_RPA,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute ppRPA excitations
!------------------------------------------------------------------------
@ -72,7 +87,7 @@ subroutine GRPA(dotest,dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas
if(doppRPA) then
call wall_time(start_RPA)
call ppGRPA(dotest,TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,epsHF)
call ppGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA

85
src/RPA/crGRPA.f90 Normal file
View File

@ -0,0 +1,85 @@
subroutine crGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
! Crossed-ring channel of the random phase approximation
implicit none
include 'parameters.h'
include 'quadrature.h'
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: TDA
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) :: ENuc
double precision,intent(in) :: EGHF
double precision,intent(in) :: eHF(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
! Hello world
write(*,*)
write(*,*)'**********************************'
write(*,*)'* Generalized cr-RPA Calculation *'
write(*,*)'**********************************'
write(*,*)
! TDA
if(TDA) then
write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*)
end if
! Initialization
dRPA = .false.
EcRPA = 0d0
! Memory allocation
allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS))
if(.not.TDA) allocate(Bph(nS,nS))
ispin = 3
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,eHF,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,EcRPA,Om,XpY,XmY)
call print_excitation_energies('crRPA@GHF',ispin,nS,Om)
call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'Tr@crRPA correlation energy = ',EcRPA,' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@crRPA total energy = ',ENuc + EGHF + EcRPA,' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
if(dotest) then
call dump_test_value('G','crRPA correlation energy',EcRPA)
end if
end subroutine

30
src/RPA/crRRPA.f90
View File

@ -1,4 +1,4 @@
subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Crossed-ring channel of the random phase approximation
@ -22,8 +22,8 @@ subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: EHF
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
@ -42,9 +42,9 @@ subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
! Hello world
write(*,*)
write(*,*)'***********************************************************'
write(*,*)'| Random phase approximation calculation: cr channel |'
write(*,*)'***********************************************************'
write(*,*)'*********************************'
write(*,*)'* Restricted cr-RPA Calculation *'
write(*,*)'*********************************'
write(*,*)
! TDA
@ -56,10 +56,8 @@ subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
! Initialization
dRPA = .false.
dRPA = .false.
EcRPA(:) = 0d0
EcRPA(:) = 0d0
! Memory allocation
@ -72,11 +70,11 @@ subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
ispin = 1
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,e,ERI,Aph)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,eHF,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,EcRPA(ispin),Om,XpY,XmY)
call print_excitation_energies('crRPA@HF',ispin,nS,Om)
call print_excitation_energies('crRPA@RHF',ispin,nS,Om)
call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
endif
@ -87,7 +85,7 @@ subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
ispin = 2
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,e,ERI,Aph)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,eHF,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,EcRPA(ispin),Om,XpY,XmY)
@ -108,7 +106,7 @@ subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
write(*,'(2X,A50,F20.10,A3)') 'Tr@crRRPA correlation energy (singlet) =',EcRPA(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@crRRPA correlation energy (triplet) =',EcRPA(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@crRRPA correlation energy =',sum(EcRPA),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@crRRPA total energy =',ENuc + EHF + sum(EcRPA),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@crRRPA total energy =',ENuc + ERHF + sum(EcRPA),' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -121,14 +119,14 @@ subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
write(*,*) '-------------------------------------------------------'
write(*,*)
call crACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcRPA)
call crACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,eHF,EcRPA)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'AC@crRRPA correlation energy (singlet) =',EcRPA(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@crRRPA correlation energy (triplet) =',EcRPA(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@crRRPA correlation energy =',sum(EcRPA),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@crRRPA total energy =',ENuc + EHF + sum(EcRPA),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@crRRPA total energy =',ENuc + ERHF + sum(EcRPA),' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -136,7 +134,7 @@ subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
if(dotest) then
call dump_test_value('R','crRRPA correlation energy',sum(EcRPA))
call dump_test_value('R','crRPA correlation energy',sum(EcRPA))
end if

17
src/RPA/phGRPA.f90
View File

@ -1,4 +1,4 @@
subroutine phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
subroutine phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
! Perform a direct random phase approximation calculation
@ -18,8 +18,8 @@ subroutine phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: EHF
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: EGHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
@ -52,8 +52,7 @@ subroutine phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
! Initialization
dRPA = .true.
dRPA = .true.
EcRPA = 0d0
! Memory allocation
@ -63,7 +62,7 @@ subroutine phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
ispin = 3
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,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,EcRPA,Om,XpY,XmY)
@ -72,14 +71,14 @@ subroutine phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@phGRPA correlation energy = ',EcRPA
write(*,'(2X,A50,F20.10)') 'Tr@phGRPA total energy = ',ENuc + EHF + EcRPA
write(*,'(2X,A50,F20.10,A3)') 'Tr@phGRPA correlation energy = ',EcRPA,' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@phGRPA total energy = ',ENuc + EGHF + EcRPA,' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
if(dotest) then
call dump_test_value('G','phGRPA corrlation energy',EcRPA)
call dump_test_value('G','phRPA corrlation energy',EcRPA)
end if

16
src/RPA/phGRPAx.f90
View File

@ -1,4 +1,4 @@
subroutine phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
subroutine phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
! Perform random phase approximation calculation with exchange (aka TDHF)
@ -18,8 +18,8 @@ subroutine phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: EHF
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: EGHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
@ -47,14 +47,12 @@ subroutine phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
if(TDA) then
write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*) ' => RPAx + TDA = CIS '
write(*,*)
end if
! Initialization
dRPA = .false.
dRPA = .false.
EcRPA = 0d0
! Memory allocation
@ -64,7 +62,7 @@ subroutine phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
ispin = 3
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,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,EcRPA,Om,XpY,XmY)
@ -74,13 +72,13 @@ subroutine phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'Tr@phRPAx correlation energy = ',EcRPA,' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@phRPAx total energy = ',ENuc + EHF + EcRPA,' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@phRPAx total energy = ',ENuc + EGHF + EcRPA,' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
if(dotest) then
call dump_test_value('G','phGRPAx correlation energy',EcRPA)
call dump_test_value('G','phRPAx correlation energy',EcRPA)
end if

19
src/RPA/phRRPA.f90
View File

@ -1,4 +1,4 @@
subroutine phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
subroutine phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform a direct random phase approximation calculation
@ -22,8 +22,8 @@ subroutine phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: EHF
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
@ -57,7 +57,6 @@ subroutine phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
! Initialization
dRPA = .true.
EcRPA(:) = 0d0
! Memory allocation
@ -71,7 +70,7 @@ subroutine phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
ispin = 1
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,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,EcRPA(ispin),Om,XpY,XmY)
@ -86,7 +85,7 @@ subroutine phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
ispin = 2
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,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,EcRPA(ispin),Om,XpY,XmY)
@ -107,7 +106,7 @@ subroutine phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
write(*,'(2X,A50,F20.10,A3)') 'Tr@phRRPA correlation energy (singlet) = ',EcRPA(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@phRRPA correlation energy (triplet) = ',EcRPA(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@phRRPA correlation energy = ',sum(EcRPA),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@phRRPA total energy = ',ENuc + EHF + sum(EcRPA),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@phRRPA total energy = ',ENuc + ERHF + sum(EcRPA),' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -122,14 +121,14 @@ subroutine phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
write(*,*) '--------------------------------------------------------'
write(*,*)
call phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcRPA)
call phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,eHF,EcRPA)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'AC@phRRPA correlation energy (singlet) = ',EcRPA(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@phRRPA correlation energy (triplet) = ',EcRPA(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@phRRPA correlation energy = ',sum(EcRPA),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@phRRPA total energy = ',ENuc + EHF + sum(EcRPA),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@phRRPA total energy = ',ENuc + ERHF + sum(EcRPA),' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -137,7 +136,7 @@ subroutine phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,
if(dotest) then
call dump_test_value('R','phRRPA correlation energy',sum(EcRPA))
call dump_test_value('R','phRPA correlation energy',sum(EcRPA))
end if

7
src/RPA/phRRPAx.f90
View File

@ -42,16 +42,15 @@ subroutine phRRPAx(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO
! Hello world
write(*,*)
write(*,*)'***********************************************************'
write(*,*)'| Random phase approximation calculation with exchange |'
write(*,*)'***********************************************************'
write(*,*)'**********************************'
write(*,*)'* Restricted ph-RPAx Calculation *'
write(*,*)'**********************************'
write(*,*)
! TDA
if(TDA) then
write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*) ' => RPAx + TDA = CIS '
write(*,*)
end if

2
src/RPA/phURPA.f90
View File

@ -167,7 +167,7 @@ subroutine phURPA(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nB
if(dotest) then
call dump_test_value('U','phURPA correlation energy',sum(EcRPA))
call dump_test_value('U','phRPA correlation energy',sum(EcRPA))
end if

35
src/RPA/ppGRPA.f90
View File

@ -1,4 +1,4 @@
subroutine ppGRPA(dotest,TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
subroutine ppGRPA(dotest,TDA,nBas,nC,nO,nV,nR,ENuc,EGHF,ERI,dipole_int,eHF)
! Perform ppGRPA calculation
@ -10,15 +10,14 @@ subroutine ppGRPA(dotest,TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
logical,intent(in) :: dotest
logical,intent(in) :: TDA
logical,intent(in) :: doACFDT
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
double precision,intent(in) :: ENuc
double precision,intent(in) :: EHF
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: EGHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
@ -60,8 +59,8 @@ subroutine ppGRPA(dotest,TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO))
if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,e,ERI,Cpp)
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,e,ERI,Dpp)
call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eHF,ERI,Cpp)
call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eHF,ERI,Dpp)
call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
@ -73,33 +72,13 @@ subroutine ppGRPA(dotest,TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppGRPA correlation energy = ',EcRPA,' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppGRPA total energy = ',ENuc + EHF + EcRPA,' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppGRPA total energy = ',ENuc + EGHF + EcRPA,' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Compute the correlation energy via the adiabatic connection
! if(doACFDT) then
! write(*,*) '--------------------------------------------------------'
! write(*,*) 'Adiabatic connection version of ppRPA correlation energy'
! write(*,*) '--------------------------------------------------------'
! write(*,*)
! call ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ERI,e,EcRPA)
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy =',EcRPA,' au'
! write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA total energy =',ENuc + EHF + EcRPA(1) + EcRPA(2),' au'
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
! end if
if(dotest) then
call dump_test_value('G','ppGRPA correlation energy',EcRPA)
call dump_test_value('G','ppRPA correlation energy',EcRPA)
end if