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

clean up in RPA

This commit is contained in:
Pierre-Francois Loos 2023-07-31 16:01:12 +02:00
parent 29f3312cfc
commit 0dfa3b0071
14 changed files with 118 additions and 125 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
input/methods
View File

@ -15,5 +15,5 @@
# G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW
F F F F F F
# G0T0pp* evGTpp* qsGTpp* G0T0eh evGTeh qsGTeh
T F F F F F
F F F F F T
# * unrestricted version available

2
input/options
View File

@ -11,7 +11,7 @@
# GW: maxSCF thresh DIIS n_diis lin eta TDA_W reg
256 0.00001 T 5 T 0.0 F F
# GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg
256 0.00001 T 5 T 0.0 F F
256 0.00001 T 5 T 0.1 T F
# ACFDT: AC Kx XBS
F F T
# BSE: phBSE phBSE2 ppBSE dBSE dTDA

32
src/CC/CC.f90
View File

@ -42,9 +42,9 @@ subroutine CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,do_drCCD,do_rCCD,do_crCCD,do_lCC
if(doCCD) then
call cpu_time(start_CC)
call wall_time(start_CC)
call CCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
call cpu_time(end_CC)
call wall_time(end_CC)
t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CC,' seconds'
@ -58,10 +58,10 @@ subroutine CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,do_drCCD,do_rCCD,do_crCCD,do_lCC
if(doDCD) then
call cpu_time(start_CC)
call wall_time(start_CC)
call DCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR, &
ERI,ENuc,EHF,epsHF)
call cpu_time(end_CC)
call wall_time(end_CC)
t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for DCD = ',t_CC,' seconds'
@ -77,9 +77,9 @@ subroutine CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,do_drCCD,do_rCCD,do_crCCD,do_lCC
if(doCCSD) then
call cpu_time(start_CC)
call wall_time(start_CC)
call CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
call cpu_time(end_CC)
call wall_time(end_CC)
t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCSD or CCSD(T)= ',t_CC,' seconds'
@ -93,9 +93,9 @@ subroutine CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,do_drCCD,do_rCCD,do_crCCD,do_lCC
if(do_drCCD) then
call cpu_time(start_CC)
call wall_time(start_CC)
call drCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
call cpu_time(end_CC)
call wall_time(end_CC)
t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for direct ring CCD = ',t_CC,' seconds'
@ -109,9 +109,9 @@ subroutine CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,do_drCCD,do_rCCD,do_crCCD,do_lCC
if(do_rCCD) then
call cpu_time(start_CC)
call wall_time(start_CC)
call rCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF,epsHF)
call cpu_time(end_CC)
call wall_time(end_CC)
t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for rCCD = ',t_CC,' seconds'
@ -125,9 +125,9 @@ subroutine CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,do_drCCD,do_rCCD,do_crCCD,do_lCC
if(do_crCCD) then
call cpu_time(start_CC)
call wall_time(start_CC)
call crCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
call cpu_time(end_CC)
call wall_time(end_CC)
t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for crossed-ring CCD = ',t_CC,' seconds'
@ -141,9 +141,9 @@ subroutine CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,do_drCCD,do_rCCD,do_crCCD,do_lCC
if(do_lCCD) then
call cpu_time(start_CC)
call wall_time(start_CC)
call lCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
call cpu_time(end_CC)
call wall_time(end_CC)
t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ladder CCD = ',t_CC,' seconds'
@ -157,9 +157,9 @@ subroutine CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,do_drCCD,do_rCCD,do_crCCD,do_lCC
if(dopCCD) then
call cpu_time(start_CC)
call wall_time(start_CC)
call pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
call cpu_time(end_CC)
call wall_time(end_CC)
t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pair CCD = ',t_CC,' seconds'

4
src/CC/CCSD.f90
View File

@ -283,9 +283,9 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBasin,nCin,nOin,nVin,nRin,ERI,EN
!------------------------------------------------------------------------
if(doCCSDT) then
call cpu_time(start_CCSDT)
call wall_time(start_CCSDT)
call CCSDT(nC,nO,nV,nR,eO,eV,OOVV,VVVO,VOOO,t1,t2,EcCCT)
call cpu_time(end_CCSDT)
call wall_time(end_CCSDT)
t_CCSDT = end_CCSDT - start_CCSDT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for (T) = ',t_CCSDT,' seconds'

16
src/CI/CI.f90
View File

@ -49,14 +49,14 @@ subroutine CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin
if(doCIS) then
call cpu_time(start_CI)
call wall_time(start_CI)
if(unrestricted) then
call UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb, &
ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
else
call CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,epsHF)
end if
call cpu_time(end_CI)
call wall_time(end_CI)
t_CI = end_CI - start_CI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CIS = ',t_CI,' seconds'
@ -70,9 +70,9 @@ subroutine CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin
if(doCID) then
call cpu_time(start_CI)
call wall_time(start_CI)
call CID(singlet,triplet,nBas,nC,nO,nV,nR,ERI,F,EHF)
call cpu_time(end_CI)
call wall_time(end_CI)
t_CI = end_CI - start_CI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CID = ',t_CI,' seconds'
@ -86,9 +86,9 @@ subroutine CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin
if(doCISD) then
call cpu_time(start_CI)
call wall_time(start_CI)
call CISD(singlet,triplet,nBas,nC,nO,nV,nR,ERI,F,EHF)
call cpu_time(end_CI)
call wall_time(end_CI)
t_CI = end_CI - start_CI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CISD = ',t_CI,' seconds'
@ -102,10 +102,10 @@ subroutine CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin
if(doFCI) then
call cpu_time(start_CI)
call wall_time(start_CI)
write(*,*) ' FCI is not yet implemented! Sorry.'
! call FCI(nBas,nC,nO,nV,nR,ERI,epsHF)
call cpu_time(end_CI)
call wall_time(end_CI)
t_CI = end_CI - start_CI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for FCI = ',t_CI,' seconds'

20
src/GF/GF.f90
View File

@ -75,7 +75,7 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t
if(doG0F2) then
call cpu_time(start_GF)
call wall_time(start_GF)
if(unrestricted) then
call UG0F2(dophBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,linearize,eta,regularize, &
nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb, &
@ -84,7 +84,7 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t
call G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, &
nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if
call cpu_time(end_GF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds'
@ -98,7 +98,7 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t
if(doevGF2) then
call cpu_time(start_GF)
call wall_time(start_GF)
if(unrestricted) then
call evUGF2(maxSCF,thresh,max_diis,dophBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb, &
@ -108,7 +108,7 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF, &
ERI,dipole_int,epsHF)
end if
call cpu_time(end_GF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds'
@ -122,7 +122,7 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t
if(doqsGF2) then
call cpu_time(start_GF)
call wall_time(start_GF)
if(unrestricted) then
call qsUGF2(maxSCF,thresh,max_diis,dophBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,regularize, &
nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO, &
@ -131,7 +131,7 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t
call qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc, &
nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF)
end if
call cpu_time(end_GF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGF2 = ',t_GF,' seconds'
@ -145,13 +145,13 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t
if(doG0F3) then
call cpu_time(start_GF)
call wall_time(start_GF)
if(unrestricted) then
print*,'Unrestricted version of G0F3 not yet implemented! Sorry.'
else
call G0F3(renorm,nBas,nC,nO,nV,nR,ERI,epsHF)
end if
call cpu_time(end_GF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF,' seconds'
@ -165,13 +165,13 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t
if(doevGF3) then
call cpu_time(start_GF)
call wall_time(start_GF)
if(unrestricted) then
print*,'Unrestricted version of evGF3 not yet implemented! Sorry.'
else
call evGF3(maxSCF,thresh,max_diis,renorm,nBas,nC,nO,nV,nR,ERI,epsHF)
end if
call cpu_time(end_GF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF,' seconds'

4
src/GT/GT.f90
View File

@ -158,8 +158,8 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted
else
call G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call soG0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
! call soG0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
! linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if
call wall_time(end_GT)

8
src/MP/MP.f90
View File

@ -37,13 +37,13 @@ subroutine MP(doMP2,doMP3,unrestricted,regularize,nBas,nC,nO,nV,nR,ERI,ERI_aaaa,
if(doMP2) then
call cpu_time(start_MP)
call wall_time(start_MP)
if(unrestricted) then
call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF)
else
call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
end if
call cpu_time(end_MP)
call wall_time(end_MP)
t_MP = end_MP - start_MP
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2 = ',t_MP,' seconds'
@ -57,7 +57,7 @@ subroutine MP(doMP2,doMP3,unrestricted,regularize,nBas,nC,nO,nV,nR,ERI,ERI_aaaa,
if(doMP3) then
call cpu_time(start_MP)
call wall_time(start_MP)
if(unrestricted) then
write(*,*) 'MP3 NYI for UHF reference'
@ -65,7 +65,7 @@ subroutine MP(doMP2,doMP3,unrestricted,regularize,nBas,nC,nO,nV,nR,ERI,ERI_aaaa,
else
call MP3(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
end if
call cpu_time(end_MP)
call wall_time(end_MP)
t_MP = end_MP - start_MP
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2 = ',t_MP,' seconds'

16
src/RPA/RPA.f90
View File

@ -52,14 +52,14 @@ subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted,
if(dophRPA) then
call cpu_time(start_RPA)
call wall_time(start_RPA)
if(unrestricted) then
call phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EHF, &
ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
else
call phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if
call cpu_time(end_RPA)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPA = ',t_RPA,' seconds'
@ -73,14 +73,14 @@ subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted,
if(dophRPAx) then
call cpu_time(start_RPA)
call wall_time(start_RPA)
if(unrestricted) then
call phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EHF, &
ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
else
call phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if
call cpu_time(end_RPA)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPAx = ',t_RPA,' seconds'
@ -94,13 +94,13 @@ subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted,
if(docrRPA) then
call cpu_time(start_RPA)
call wall_time(start_RPA)
if(unrestricted) then
write(*,*) 'Unrestricted version of crRPA not yet implemented! Sorry.'
else
call crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if
call cpu_time(end_RPA)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds'
@ -114,13 +114,13 @@ subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted,
if(doppRPA) then
call cpu_time(start_RPA)
call wall_time(start_RPA)
if(unrestricted) then
call ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,epsHF)
else
call ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,epsHF)
end if
call cpu_time(end_RPA)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds'

33
src/RPA/crRPA.f90
View File

@ -35,8 +35,7 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS
double precision,allocatable :: XpY(:,:)
double precision,allocatable :: XmY(:,:)
double precision :: EcTr(nspin)
double precision :: EcAC(nspin)
double precision :: EcRPA(nspin)
! Hello world
@ -57,8 +56,8 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS
dRPA = .false.
EcTr(:) = 0d0
EcAC(:) = 0d0
EcRPA(:) = 0d0
EcRPA(:) = 0d0
! Memory allocation
@ -74,7 +73,7 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS
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 phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
call print_excitation_energies('crRPA@HF',ispin,nS,Om)
call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
@ -89,7 +88,7 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS
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 phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
call print_excitation_energies('crRPA@HF',ispin,nS,Om)
call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
@ -97,17 +96,17 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS
if(exchange_kernel) then
EcTr(1) = 0.5d0*EcTr(1)
EcTr(2) = 1.5d0*EcTr(2)
EcRPA(1) = 0.5d0*EcRPA(1)
EcRPA(2) = 1.5d0*EcRPA(2)
end if
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
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(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@crRPA total energy =',ENuc + EHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -120,14 +119,14 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS
write(*,*) '-------------------------------------------------------'
write(*,*)
call crACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
call crACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcRPA)
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 + EHF + EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'AC@crRPA total energy =',ENuc + EHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)

32
src/RPA/phRPA.f90
View File

@ -35,8 +35,7 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS
double precision,allocatable :: XpY(:,:)
double precision,allocatable :: XmY(:,:)
double precision :: EcTr(nspin)
double precision :: EcAC(nspin)
double precision :: EcRPA(nspin)
! Hello world
@ -57,8 +56,7 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS
dRPA = .true.
EcTr(:) = 0d0
EcAC(:) = 0d0
EcRPA(:) = 0d0
! Memory allocation
@ -74,7 +72,7 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS
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 phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
call print_excitation_energies('phRPA@HF',ispin,nS,Om)
call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
@ -89,7 +87,7 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS
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 phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
call print_excitation_energies('phRPA@HF',ispin,nS,Om)
call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
@ -97,17 +95,17 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS
if(exchange_kernel) then
EcTr(1) = 0.5d0*EcTr(1)
EcTr(2) = 1.5d0*EcTr(2)
EcRPA(1) = 0.5d0*EcRPA(1)
EcRPA(2) = 1.5d0*EcRPA(2)
end if
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
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(*,'(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 + EHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -122,14 +120,14 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS
write(*,*) '--------------------------------------------------------'
write(*,*)
call phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
call phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcRPA)
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 + EHF + EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@phRPA correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'AC@phRPA correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'AC@phRPA correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'AC@phRPA total energy =',ENuc + EHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)

32
src/RPA/phRPAx.f90
View File

@ -35,8 +35,7 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,n
double precision,allocatable :: XpY(:,:)
double precision,allocatable :: XmY(:,:)
double precision :: EcTr(nspin)
double precision :: EcAC(nspin)
double precision :: EcRPA(nspin)
! Hello world
@ -58,8 +57,7 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,n
dRPA = .false.
EcTr(:) = 0d0
EcAC(:) = 0d0
EcRPA(:) = 0d0
! Memory allocation
@ -75,7 +73,7 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,n
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 phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
call print_excitation_energies('phRPAx@HF',ispin,nS,Om)
call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
@ -90,7 +88,7 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,n
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 phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
call print_excitation_energies('phRPAx@HF',ispin,nS,Om)
call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
@ -98,17 +96,17 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,n
if(exchange_kernel) then
EcTr(1) = 0.5d0*EcTr(1)
EcTr(2) = 1.5d0*EcTr(2)
EcRPA(1) = 0.5d0*EcRPA(1)
EcRPA(2) = 1.5d0*EcRPA(2)
end if
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
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(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@phRPAx total energy =',ENuc + EHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -125,14 +123,14 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,n
write(*,*) '---------------------------------------------------------'
write(*,*)
call phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
call phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcRPA)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
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(*,'(2X,A50,F20.10)') 'AC@phRPAx correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'AC@phRPAx correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'AC@phRPAx correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'AC@phRPAx total energy =',ENuc + EHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)

12
src/RPA/ppRPA.f90
View File

@ -38,7 +38,6 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipol
double precision,allocatable :: Y2(:,:)
double precision :: EcRPA(nspin)
double precision :: EcAC(nspin)
! Hello world
@ -51,7 +50,6 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipol
! Initialization
EcRPA(:) = 0d0
EcAC(:) = 0d0
! Singlet manifold
@ -135,14 +133,14 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipol
write(*,*) '--------------------------------------------------------'
write(*,*)
call ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ERI,e,EcAC)
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 (singlet) =',EcAC(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (triplet) =',EcAC(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy =',EcAC(1) + EcAC(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA total energy =',ENuc + EHF + EcAC(1) + EcAC(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (singlet) =',EcRPA(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (triplet) =',EcRPA(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy =',EcRPA(1) + EcRPA(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA total energy =',ENuc + EHF + EcRPA(1) + EcRPA(2),' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)

30
src/RPA/sort_ppRPA.f90
View File

@ -1,4 +1,4 @@
subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
! Compute the metric matrix for pp-RPA
@ -9,7 +9,7 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
integer,intent(in) :: nOO
integer,intent(in) :: nVV
double precision,intent(in) :: Omega(nOO+nVV)
double precision,intent(in) :: Om(nOO+nVV)
double precision,intent(in) :: Z(nOO+nVV,nOO+nVV)
! Local variables
@ -32,10 +32,10 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
! Output variables
double precision,intent(out) :: Omega1(nVV)
double precision,intent(out) :: Om1(nVV)
double precision,intent(out) :: X1(nVV,nVV)
double precision,intent(out) :: Y1(nOO,nVV)
double precision,intent(out) :: Omega2(nOO)
double precision,intent(out) :: Om2(nOO)
double precision,intent(out) :: X2(nVV,nOO)
double precision,intent(out) :: Y2(nOO,nOO)
@ -47,11 +47,11 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
! Initializatiom
Omega1(:) = 0d0
Om1(:) = 0d0
X1(:,:) = 0d0
Y1(:,:) = 0d0
Omega2(:) = 0d0
Om2(:) = 0d0
X2(:,:) = 0d0
Y2(:,:) = 0d0
@ -74,24 +74,24 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
do pq=1,nOO+nVV
if(Omega(pq) > 0d0) then
if(Om(pq) > 0d0) then
ab = ab + 1
Omega1(ab) = Omega(pq)
Om1(ab) = Om(pq)
Z1(1:nOO+nVV,ab) = Z(1:nOO+nVV,pq)
else
ij = ij + 1
Omega2(ij) = Omega(pq)
Om2(ij) = Om(pq)
Z2(1:nOO+nVV,ij) = Z(1:nOO+nVV,pq)
end if
end do
if(minval(Omega1(:)) < 0d0 .or. ab /= nVV) call print_warning('You may have instabilities in pp-RPA!!')
if(maxval(Omega2(:)) > 0d0 .or. ij /= nOO) call print_warning('You may have instabilities in pp-RPA!!')
if(minval(Om1(:)) < 0d0 .or. ab /= nVV) call print_warning('You may have instabilities in pp-RPA!!')
if(maxval(Om2(:)) > 0d0 .or. ij /= nOO) call print_warning('You may have instabilities in pp-RPA!!')
if(nVV > 0) then
@ -100,7 +100,7 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
order1(ab) = ab
end do
call quick_sort(Omega1(:),order1(:),nVV)
call quick_sort(Om1(:),order1(:),nVV)
call set_order(Z1(:,:),order1(:),nOO+nVV,nVV)
end if
@ -111,7 +111,7 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
order2(ij) = ij
end do
call quick_sort(Omega2(:),order2(:),nOO)
call quick_sort(Om2(:),order2(:),nOO)
call set_order(Z2(:,:),order2(:),nOO+nVV,nOO)
end if
@ -123,7 +123,7 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
! do ab=2,nVV
! if(abs(Omega1(ab-1) - Omega1(ab)) < 1d-6) then
! if(abs(Om1(ab-1) - Om1(ab)) < 1d-6) then
! deg1 = deg1 + 1
@ -164,7 +164,7 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
! do ij=2,nOO
! if(abs(Omega2(ij-1) - Omega2(ij)) < 1d-6) then
! if(abs(Om2(ij-1) - Om2(ij)) < 1d-6) then
! deg2 = deg2 + 1