PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-03 20:53:53 +01:00
Code Releases Activity

introduce test in MP and RPA

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-11-11 23:00:00 +01:00
parent ef528e6a3c
commit 436231df11
34 changed files with 346 additions and 192 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
input/methods
View File

@ -1,7 +1,7 @@
# RHF UHF GHF ROHF
T T T T
# MP2 MP3
T F
T T
# 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 F T
# G0F2 evGF2 qsGF2 G0F3 evGF3
F F F F F
# G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW
@ -17,4 +17,4 @@
# G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
F F F F F F
# Rtest Utest Gtest
F F F
T T T

6
src/HF/GHF.f90
View File

@ -1,4 +1,4 @@
subroutine GHF(doGtest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
subroutine GHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nBas2,nO,Ov,T,V,Hc,ERI,dipole_int,Or,EHF,e,c,P)
! Perform unrestricted Hartree-Fock calculation
@ -8,7 +8,7 @@ subroutine GHF(doGtest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
! Input variables
logical,intent(in) :: doGtest
logical,intent(in) :: dotest
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
@ -298,7 +298,7 @@ subroutine GHF(doGtest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
! Print test values
if(doGtest) then
if(dotest) then
call dump_test_value('G','GHF energy',EHF)

6
src/HF/RHF.f90
View File

@ -1,4 +1,4 @@
subroutine RHF(doRtest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,e,c,P)
! Perform restricted Hartree-Fock calculation
@ -8,7 +8,7 @@ subroutine RHF(doRtest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
! Input variables
logical,intent(in) :: doRtest
logical,intent(in) :: dotest
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
@ -203,7 +203,7 @@ subroutine RHF(doRtest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
! Print test values
if(doRtest) then
if(dotest) then
call dump_test_value('R','RHF energy',EHF)

6
src/HF/ROHF.f90
View File

@ -1,4 +1,4 @@
subroutine ROHF(doRtest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,e,c,Ptot)
! Perform restricted open-shell Hartree-Fock calculation
@ -8,7 +8,7 @@ subroutine ROHF(doRtest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,Z
! Input variables
logical,intent(in) :: doRtest
logical,intent(in) :: dotest
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
@ -245,7 +245,7 @@ subroutine ROHF(doRtest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,Z
! Print test values
if(doRtest) then
if(dotest) then
call dump_test_value('R','ROHF energy',EHF)

6
src/HF/UHF.f90
View File

@ -1,4 +1,4 @@
subroutine UHF(doUtest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,e,c,P)
! Perform unrestricted Hartree-Fock calculation
@ -8,7 +8,7 @@ subroutine UHF(doUtest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
! Input variables
logical,intent(in) :: doUtest
logical,intent(in) :: dotest
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
@ -256,7 +256,7 @@ subroutine UHF(doUtest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
! Print test values
if(doUtest) then
if(dotest) then
call dump_test_value('U','UHF energy',EHF)

6
src/MP/GMP.f90
View File

@ -1,4 +1,4 @@
subroutine GMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
subroutine GMP(dotest,doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
! Moller-Plesset module
@ -7,6 +7,8 @@ subroutine GMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: doMP2
logical,intent(in) :: doMP3
@ -35,7 +37,7 @@ subroutine GMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
if(doMP2) then
call wall_time(start_MP)
call GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF,Ec)
call GMP2(dotest,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF,Ec)
call wall_time(end_MP)
t_MP = end_MP - start_MP

39
src/MP/GMP2.f90
View File

@ -1,4 +1,4 @@
subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
subroutine GMP2(dotest,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF,EcMP2)
! Perform second-order Moller-Plesset calculation with and without regularizers
@ -6,6 +6,8 @@ subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nC
@ -13,8 +15,8 @@ subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
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)
! Local variables
@ -38,9 +40,9 @@ subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Moller-Plesset second-order calculation |'
write(*,*)'************************************************'
write(*,*)'*******************************'
write(*,*)'* Generalized MP2 Calculation |'
write(*,*)'*******************************'
write(*,*)
!---------------------------------------------!
@ -70,7 +72,7 @@ subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
Dijab = e(a) + e(b) - e(i) - e(j)
Dijab = eHF(a) + eHF(b) - eHF(i) - eHF(j)
! Second-order ring diagram
@ -114,8 +116,8 @@ subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
write(*,'(A32,1X,F16.10)') ' Direct part = ',E2d
write(*,'(A32,1X,F16.10)') ' Exchange part = ',E2x
write(*,'(A32)') '---------------------------'
write(*,'(A32,1X,F16.10)') ' GMP2 electronic energy = ',EHF + EcMP2
write(*,'(A32,1X,F16.10)') ' GMP2 total energy = ',ENuc + EHF + EcMP2
write(*,'(A32,1X,F16.10)') ' GMP2 electronic energy = ',EGHF + EcMP2
write(*,'(A32,1X,F16.10)') ' GMP2 total energy = ',ENuc + EGHF + EcMP2
write(*,'(A32)') '---------------------------'
write(*,*)
@ -133,8 +135,8 @@ subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
write(*,'(A32,1X,F16.10)') ' Direct part = ',E2ds
write(*,'(A32,1X,F16.10)') ' Exchange part = ',E2xs
write(*,'(A32)') '---------------------------'
write(*,'(A32,1X,F16.10)') ' GMP2 electronic energy = ',EHF + EcsMP2
write(*,'(A32,1X,F16.10)') ' GMP2 total energy = ',ENuc + EHF + EcsMP2
write(*,'(A32,1X,F16.10)') ' GMP2 electronic energy = ',EGHF + EcsMP2
write(*,'(A32,1X,F16.10)') ' GMP2 total energy = ',ENuc + EGHF + EcsMP2
write(*,'(A32)') '---------------------------'
write(*,*)
@ -150,8 +152,8 @@ subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
write(*,'(A32,1X,F16.10)') ' Direct part = ',E2ds2
write(*,'(A32,1X,F16.10)') ' Exchange part = ',E2xs2
write(*,'(A32)') '---------------------------'
write(*,'(A32,1X,F16.10)') ' GMP2 electronic energy = ',EHF + Ecs2MP2
write(*,'(A32,1X,F16.10)') ' GMP2 total energy = ',ENuc + EHF + Ecs2MP2
write(*,'(A32,1X,F16.10)') ' GMP2 electronic energy = ',EGHF + Ecs2MP2
write(*,'(A32,1X,F16.10)') ' GMP2 total energy = ',ENuc + EGHF + Ecs2MP2
write(*,'(A32)') '---------------------------'
write(*,*)
@ -167,11 +169,18 @@ subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
write(*,'(A32,1X,F16.10)') ' Direct part = ',E2dk
write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2xk
write(*,'(A32)') '---------------------------'
write(*,'(A32,1X,F16.10)') ' GMP2 electronic energy = ',EHF + EckMP2
write(*,'(A32,1X,F16.10)') ' GMP2 total energy = ',ENuc + EHF + EckMP2
write(*,'(A32,1X,F16.10)') ' GMP2 electronic energy = ',EGHF + EckMP2
write(*,'(A32,1X,F16.10)') ' GMP2 total energy = ',ENuc + EGHF + EckMP2
write(*,'(A32)') '---------------------------'
write(*,*)
end if
if(dotest) then
call dump_test_value('G','GMP2 correlation energy',EcMP2)
end if
end subroutine

6
src/MP/GMP3.f90
View File

@ -38,9 +38,9 @@ subroutine GMP3(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Moller-Plesset third-order calculation |'
write(*,*)'************************************************'
write(*,*)'*******************************'
write(*,*)'* Generalized MP3 Calculation *'
write(*,*)'*******************************'
write(*,*)
! Antysymmetrize ERIs

6
src/MP/RMP.f90
View File

@ -1,4 +1,4 @@
subroutine RMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
subroutine RMP(dotest,doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
! Moller-Plesset module
@ -7,6 +7,8 @@ subroutine RMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: doMP2
logical,intent(in) :: doMP3
@ -35,7 +37,7 @@ subroutine RMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
if(doMP2) then
call wall_time(start_MP)
call RMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF,Ec)
call RMP2(dotest,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF,Ec)
call wall_time(end_MP)
t_MP = end_MP - start_MP

16
src/MP/RMP2.f90
View File

@ -1,4 +1,4 @@
subroutine RMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF,EcMP2)
subroutine RMP2(dotest,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF,EcMP2)
! Perform second-order Moller-Plesset calculation with and without regularizers
@ -6,6 +6,8 @@ subroutine RMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF,EcMP2)
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nC
@ -37,9 +39,9 @@ subroutine RMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF,EcMP2)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Moller-Plesset second-order calculation |'
write(*,*)'************************************************'
write(*,*)'******************************'
write(*,*)'* Restricted MP2 Calculation *'
write(*,*)'******************************'
write(*,*)
!---------------------------------------------!
@ -171,4 +173,10 @@ subroutine RMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF,EcMP2)
end if
if(dotest) then
call dump_test_value('R','RMP2 correlation energy',EcMP2)
end if
end subroutine

6
src/MP/RMP3.f90
View File

@ -45,9 +45,9 @@ subroutine RMP3(nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,EHF,e)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Moller-Plesset third-order calculation |'
write(*,*)'************************************************'
write(*,*)'******************************'
write(*,*)'* Restricted MP3 Calculation *'
write(*,*)'******************************'
write(*,*)
! Spatial to spin orbitals

6
src/MP/UMP.f90
View File

@ -1,4 +1,4 @@
subroutine UMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF)
subroutine UMP(dotest,doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF)
! Moller-Plesset module
@ -7,6 +7,8 @@ subroutine UMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbb
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: doMP2
logical,intent(in) :: doMP3
@ -37,7 +39,7 @@ subroutine UMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbb
if(doMP2) then
call wall_time(start_MP)
call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF,Ec)
call UMP2(dotest,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF,Ec)
call wall_time(end_MP)
t_MP = end_MP - start_MP

31
src/MP/UMP2.f90
View File

@ -1,24 +1,25 @@
subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec)
subroutine UMP2(dotest,nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EUHF,eHF,Ec)
! Perform unrestricted second-order Moller-Plesset calculation
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dotest
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
double precision,intent(in) :: ENuc
double precision,intent(in) :: EHF
double precision,intent(in) :: EUHF
double precision,intent(in) :: ERI_aa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_ab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: e(nBas,nspin)
double precision,intent(in) :: eHF(nBas,nspin)
! Local variables
@ -37,9 +38,9 @@ subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec)
! Hello world
write(*,*)
write(*,*)'********************************************************'
write(*,*)'| Unrestricted second-order Moller-Plesset calculation |'
write(*,*)'********************************************************'
write(*,*)'********************************'
write(*,*)'* Unrestricted MP2 Calculation *'
write(*,*)'********************************'
write(*,*)
!---------------------!
@ -60,7 +61,7 @@ subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec)
do j=nC(ket)+1,nO(ket)
do b=nO(ket)+1,nBas-nR(ket)
eps = e(i,bra) + e(j,ket) - e(a,bra) - e(b,ket)
eps = eHF(i,bra) + eHF(j,ket) - eHF(a,bra) - eHF(b,ket)
Edaa = Edaa + 0.5d0*ERI_aa(i,j,a,b)*ERI_aa(i,j,a,b)/eps
Exaa = Exaa - 0.5d0*ERI_aa(i,j,a,b)*ERI_aa(i,j,b,a)/eps
@ -88,7 +89,7 @@ subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec)
do j=nC(ket)+1,nO(ket)
do b=nO(ket)+1,nBas-nR(ket)
eps = e(i,bra) + e(j,ket) - e(a,bra) - e(b,ket)
eps = eHF(i,bra) + eHF(j,ket) - eHF(a,bra) - eHF(b,ket)
Edab = Edab + ERI_ab(i,j,a,b)*ERI_ab(i,j,a,b)/eps
@ -114,7 +115,7 @@ subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec)
do j=nC(ket)+1,nO(ket)
do b=nO(ket)+1,nBas-nR(ket)
eps = e(i,bra) + e(j,ket) - e(a,bra) - e(b,ket)
eps = eHF(i,bra) + eHF(j,ket) - eHF(a,bra) - eHF(b,ket)
Edbb = Edbb + 0.5d0*ERI_bb(i,j,a,b)*ERI_bb(i,j,a,b)/eps
Exbb = Exbb - 0.5d0*ERI_bb(i,j,a,b)*ERI_bb(i,j,b,a)/eps
@ -152,9 +153,15 @@ subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec)
write(*,'(A32,1X,F16.10)') ' alpha-beta = ',Exab
write(*,'(A32,1X,F16.10)') ' beta-beta = ',Exbb
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ', EHF + sum(Ec(:))
write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + sum(Ec(:))
write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ', EUHF + sum(Ec(:))
write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EUHF + sum(Ec(:))
write(*,'(A32)') '--------------------------'
write(*,*)
if(dotest) then
call dump_test_value('U','UMP2 correlation energy',sum(Ec))
end if
end subroutine

10
src/QuAcK/GQuAcK.f90
View File

@ -1,4 +1,4 @@
subroutine GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,doppRPA, &
subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,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, &
@ -9,7 +9,7 @@ subroutine GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dop
implicit none
include 'parameters.h'
logical,intent(in) :: doGtest
logical,intent(in) :: dotest
logical,intent(in) :: doGHF
logical,intent(in) :: dostab
@ -103,7 +103,7 @@ subroutine GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dop
if(doGHF) then
call wall_time(start_HF)
call GHF(doGtest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
call GHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nBas2,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
call wall_time(end_HF)
@ -196,7 +196,7 @@ subroutine GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dop
if(doMP) then
call wall_time(start_MP)
call GMP(doMP2,doMP3,reg_MP,nBas2,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
call GMP(dotest,doMP2,doMP3,reg_MP,nBas2,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
call wall_time(end_MP)
t_MP = end_MP - start_MP
@ -214,7 +214,7 @@ subroutine GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,dop
if(doRPA) then
call wall_time(start_RPA)
call GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas2,nC,nO,nV,nR,nS,ENuc,EHF, &
call GRPA(dotest,dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas2,nC,nO,nV,nR,nS,ENuc,EHF, &
ERI_MO,dipole_int_MO,epsHF,cHF,S)
call wall_time(end_RPA)

13
src/QuAcK/RQuAcK.f90
View File

@ -1,4 +1,4 @@
subroutine RQuAcK(doRtest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, &
doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
@ -13,7 +13,7 @@ subroutine RQuAcK(doRtest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,
implicit none
include 'parameters.h'
logical,intent(in) :: doRtest
logical,intent(in) :: dotest
logical,intent(in) :: doRHF,doROHF
logical,intent(in) :: dostab
@ -117,11 +117,12 @@ subroutine RQuAcK(doRtest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,
if(doRHF) then
call wall_time(start_HF)
call RHF(doRtest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
call RHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
call wall_time(end_HF)
t_HF = end_HF - start_HF
write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RHF = ',t_HF,' seconds'
write(*,*)
@ -130,7 +131,7 @@ subroutine RQuAcK(doRtest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,
if(doROHF) then
call wall_time(start_HF)
call ROHF(doRtest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
call ROHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
call wall_time(end_HF)
@ -206,7 +207,7 @@ subroutine RQuAcK(doRtest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,
if(doMP) then
call wall_time(start_MP)
call RMP(doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
call RMP(dotest,doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
call wall_time(end_MP)
t_MP = end_MP - start_MP
@ -263,7 +264,7 @@ subroutine RQuAcK(doRtest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,
if(doRPA) then
call wall_time(start_RPA)
call RRPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, &
call RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, &
nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF,cHF,S)
call wall_time(end_RPA)

10
src/QuAcK/UQuAcK.f90
View File

@ -1,4 +1,4 @@
subroutine UQuAcK(doUtest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, &
doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
@ -11,7 +11,7 @@ subroutine UQuAcK(doUtest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,d
implicit none
include 'parameters.h'
logical,intent(in) :: doUtest
logical,intent(in) :: dotest
logical,intent(in) :: doUHF
logical,intent(in) :: dostab
@ -116,7 +116,7 @@ subroutine UQuAcK(doUtest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,d
if(doUHF) then
call wall_time(start_HF)
call UHF(doUtest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
call UHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
call wall_time(end_HF)
@ -223,7 +223,7 @@ subroutine UQuAcK(doUtest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,d
if(doMP) then
call wall_time(start_MP)
call UMP(doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF)
call UMP(dotest,doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF)
call wall_time(end_MP)
t_MP = end_MP - start_MP
@ -278,7 +278,7 @@ subroutine UQuAcK(doUtest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,d
if(doRPA) then
call wall_time(start_RPA)
call URPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip, &
call URPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,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)
call wall_time(end_RPA)

16
src/RPA/GRPA.f90
View File

@ -1,4 +1,4 @@
subroutine GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
subroutine GRPA(dotest,dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
! Random-phase approximation module
@ -7,9 +7,11 @@ subroutine GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas,nC,nO,
! Input variables
logical :: dophRPA
logical :: dophRPAx
logical :: doppRPA
logical,intent(in) :: dotest
logical,intent(in) :: dophRPA
logical,intent(in) :: dophRPAx
logical,intent(in) :: doppRPA
logical,intent(in) :: TDA
logical,intent(in) :: doACFDT
@ -37,7 +39,7 @@ subroutine GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas,nC,nO,
if(dophRPA) then
call wall_time(start_RPA)
call phGRPA(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
@ -53,7 +55,7 @@ subroutine GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas,nC,nO,
if(dophRPAx) then
call wall_time(start_RPA)
call phGRPAx(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_RPA)
@ -70,7 +72,7 @@ subroutine GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas,nC,nO,
if(doppRPA) then
call wall_time(start_RPA)
call ppGRPA(TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,epsHF)
call ppGRPA(dotest,TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA

20
src/RPA/RRPA.f90
View File

@ -1,4 +1,4 @@
subroutine RRPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, &
subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, &
nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF,cHF,S)
! Random-phase approximation module
@ -8,10 +8,12 @@ subroutine RRPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,sin
! Input variables
logical :: dophRPA
logical :: dophRPAx
logical :: docrRPA
logical :: doppRPA
logical,intent(in) :: dotest
logical,intent(in) :: dophRPA
logical,intent(in) :: dophRPAx
logical,intent(in) :: docrRPA
logical,intent(in) :: doppRPA
logical,intent(in) :: TDA
logical,intent(in) :: doACFDT
@ -43,7 +45,7 @@ subroutine RRPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,sin
if(dophRPA) then
call wall_time(start_RPA)
call phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
@ -59,7 +61,7 @@ subroutine RRPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,sin
if(dophRPAx) then
call wall_time(start_RPA)
call phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call phRRPAx(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
@ -75,7 +77,7 @@ subroutine RRPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,sin
if(docrRPA) then
call wall_time(start_RPA)
call crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
@ -91,7 +93,7 @@ subroutine RRPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,sin
if(doppRPA) then
call wall_time(start_RPA)
call ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,epsHF)
call ppRRPA(dotest,TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA

18
src/RPA/URPA.f90
View File

@ -1,4 +1,4 @@
subroutine URPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip, &
subroutine URPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,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)
! Random-phase approximation module
@ -8,10 +8,12 @@ subroutine URPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,spi
! Input variables
logical :: dophRPA
logical :: dophRPAx
logical :: docrRPA
logical :: doppRPA
logical,intent(in) :: dotest
logical,intent(in) :: dophRPA
logical,intent(in) :: dophRPAx
logical,intent(in) :: docrRPA
logical,intent(in) :: doppRPA
logical,intent(in) :: TDA
logical,intent(in) :: doACFDT
@ -46,7 +48,7 @@ subroutine URPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,spi
if(dophRPA) then
call wall_time(start_RPA)
call phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EHF, &
call phURPA(dotest,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)
call wall_time(end_RPA)
@ -63,7 +65,7 @@ subroutine URPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,spi
if(dophRPAx) then
call wall_time(start_RPA)
call phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EHF, &
call phURPAx(dotest,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)
call wall_time(end_RPA)
@ -96,7 +98,7 @@ subroutine URPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,spi
if(doppRPA) then
call wall_time(start_RPA)
call ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,epsHF)
call ppURPA(dotest,TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,epsHF)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA

4
src/RPA/crRPA.f90 → src/RPA/crRRPA.f90
View File

@ -1,4 +1,4 @@
subroutine crRPA(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,EHF,ERI,dipole_int,e)
! Crossed-ring channel of the random phase approximation
@ -8,6 +8,8 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: TDA
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel

20
src/RPA/phGRPA.f90
View File

@ -1,4 +1,4 @@
subroutine phGRPA(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
subroutine phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
! Perform a direct random phase approximation calculation
@ -8,6 +8,8 @@ subroutine phGRPA(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: TDA
integer,intent(in) :: nBas
integer,intent(in) :: nC
@ -36,9 +38,9 @@ subroutine phGRPA(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
! Hello world
write(*,*)
write(*,*)'***********************************************'
write(*,*)'| Random-phase approximation calculation |'
write(*,*)'***********************************************'
write(*,*)'**********************************'
write(*,*)'* Generalized ph-RPA Calculation |'
write(*,*)'**********************************'
write(*,*)
! TDA
@ -70,9 +72,15 @@ subroutine phGRPA(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy =',EcRPA
write(*,'(2X,A50,F20.10)') 'Tr@phRPA total energy =',ENuc + EHF + EcRPA
write(*,'(2X,A50,F20.10)') 'Tr@phGRPA correlation energy = ',EcRPA
write(*,'(2X,A50,F20.10)') 'Tr@phGRPA total energy = ',ENuc + EHF + EcRPA
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
if(dotest) then
call dump_test_value('G','phGRPA corrlation energy',EcRPA)
end if
end subroutine

20
src/RPA/phGRPAx.f90
View File

@ -1,4 +1,4 @@
subroutine phGRPAx(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
subroutine phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
! Perform random phase approximation calculation with exchange (aka TDHF)
@ -8,6 +8,8 @@ subroutine phGRPAx(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: TDA
integer,intent(in) :: nBas
integer,intent(in) :: nC
@ -36,9 +38,9 @@ subroutine phGRPAx(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
! Hello world
write(*,*)
write(*,*)'***********************************************************'
write(*,*)'| Random phase approximation calculation with exchange |'
write(*,*)'***********************************************************'
write(*,*)'************************************'
write(*,*)'* Generalized ph-RPAx Calculation *'
write(*,*)'************************************'
write(*,*)
! TDA
@ -71,9 +73,15 @@ subroutine phGRPAx(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy =',EcRPA
write(*,'(2X,A50,F20.10)') 'Tr@phRPAx total energy =',ENuc + EHF + EcRPA
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(*,*)'-------------------------------------------------------------------------------'
write(*,*)
if(dotest) then
call dump_test_value('G','phGRPAx correlation energy',EcRPA)
end if
end subroutine

26
src/RPA/phRPA.f90 → src/RPA/phRRPA.f90
View File

@ -1,4 +1,4 @@
subroutine phRPA(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,EHF,ERI,dipole_int,e)
! Perform a direct random phase approximation calculation
@ -8,6 +8,8 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: TDA
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
@ -102,10 +104,10 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS
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 + EHF + EcRPA(1) + EcRPA(2)
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(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -124,13 +126,19 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
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(*,'(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(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end if
if(dotest) then
call dump_test_value('R','phRRPA correlation energy',sum(EcRPA))
end if
end subroutine

26
src/RPA/phRPAx.f90 → src/RPA/phRRPAx.f90
View File

@ -1,4 +1,4 @@
subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
subroutine phRRPAx(dotest,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)
@ -8,6 +8,8 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,n
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: TDA
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
@ -103,10 +105,10 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,n
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
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(*,'(2X,A50,F20.10,A3)') 'Tr@phRRPAx correlation energy (singlet) = ',EcRPA(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@phRRPAx correlation energy (triplet) = ',EcRPA(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@phRRPAx correlation energy = ',sum(EcRPA),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@phRRPAx total energy = ',ENuc + EHF + sum(EcRPA),' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -127,13 +129,19 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,n
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
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(*,'(2X,A50,F20.10,A3)') 'AC@phRRPAx correlation energy (singlet) = ',EcRPA(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@phRRPAx correlation energy (triplet) = ',EcRPA(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@phRRPAx correlation energy = ',sum(EcRPA),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@phRRPAx total energy = ',ENuc + EHF + sum(EcRPA),' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end if
if(dotest) then
call dump_test_value('R','phRRPAx correlation energy',sum(EcRPA))
end if
end subroutine

38
src/RPA/phURPA.f90
View File

@ -1,5 +1,5 @@
subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,e,c,S)
subroutine phURPA(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,eHF,c,S)
! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism
@ -9,6 +9,8 @@ subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,n
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: TDA
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
@ -22,7 +24,7 @@ subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,n
integer,intent(in) :: nS(nspin)
double precision,intent(in) :: ENuc
double precision,intent(in) :: EUHF
double precision,intent(in) :: e(nBas,nspin)
double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: c(nBas,nBas,nspin)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
@ -81,7 +83,7 @@ subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,n
allocate(Aph(nSt,nSt),Bph(nSt,nSt),Om(nSt),XpY(nSt,nSt),XmY(nSt,nSt))
call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
call phULR_A(ispin,dRPA,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)
call phULR(TDA,nSa,nSb,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
@ -104,9 +106,9 @@ subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,n
nSb = (nO(2) - nC(2))*(nV(1) - nR(1))
nSt = nSa + nSb
allocate(Om(nSt),XpY(nSt,nSt),XmY(nSt,nSt))
allocate(Aph(nSt,nSt),Bph(nSt,nSt),Om(nSt),XpY(nSt,nSt),XmY(nSt,nSt))
call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
call phULR_A(ispin,dRPA,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)
call phULR(TDA,nSa,nSa,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
@ -126,10 +128,10 @@ subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,n
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@URPA correlation energy (spin-conserved) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'Tr@URPA correlation energy (spin-flip) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@URPA correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@URPA total energy =',ENuc + EUHF + EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10,A3)') 'Tr@phURPA correlation energy (spin-conserved) = ',EcRPA(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@phURPA correlation energy (spin-flip) = ',EcRPA(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@phURPA correlation energy = ',sum(EcRPA),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@phURPA total energy = ',ENuc + EUHF + sum(EcRPA),' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -143,7 +145,7 @@ subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,n
write(*,*)
call phUACFDT(exchange_kernel,.false.,.true.,.false.,TDA,.false.,spin_conserved,spin_flip, &
nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,e,e,EcRPA)
nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eHF,EcRPA)
if(exchange_kernel) then
@ -154,13 +156,19 @@ subroutine phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,n
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'AC@URPA correlation energy (spin-conserved) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'AC@URPA correlation energy (spin-flip) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'AC@URPA correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'AC@URPA total energy =',ENuc + EUHF + EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10,A3)') 'AC@phURPA correlation energy (spin-conserved) = ',EcRPA(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@phURPA correlation energy (spin-flip) = ',EcRPA(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@phURPA correlation energy = ',sum(EcRPA),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@phURPA total energy = ',ENuc + EUHF + sum(EcRPA),' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end if
if(dotest) then
call dump_test_value('U','phURPA correlation energy',sum(EcRPA))
end if
end subroutine

28
src/RPA/phURPAx.f90
View File

@ -1,4 +1,4 @@
subroutine phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
subroutine phURPAx(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,e,c,S)
! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism
@ -9,6 +9,8 @@ subroutine phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: TDA
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
@ -104,7 +106,7 @@ subroutine phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,
nSb = (nO(2) - nC(2))*(nV(1) - nR(1))
nSt = nSa + nSb
allocate(Om(nSt),XpY(nSt,nSt),XmY(nSt,nSt))
allocate(Aph(nSt,nSt),Bph(nSt,nSt),Om(nSt),XpY(nSt,nSt),XmY(nSt,nSt))
call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,e,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)
@ -137,10 +139,10 @@ subroutine phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@URPAx correlation energy (spin-conserved) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'Tr@URPAx correlation energy (spin-flip) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@URPAx correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@URPAx total energy =',ENuc + EUHF + EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10,A3)') 'Tr@URPAx correlation energy (spin-conserved) = ',EcRPA(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@URPAx correlation energy (spin-flip) = ',EcRPA(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@URPAx correlation energy = ',sum(EcRPA),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@URPAx total energy = ',ENuc + EUHF + sum(EcRPA),' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -158,13 +160,19 @@ subroutine phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'AC@URPAx correlation energy (spin-conserved) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'AC@URPAx correlation energy (spin-flip) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'AC@URPAx correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'AC@URPAx total energy =',ENuc + EUHF + EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10,A3)') 'AC@URPAx correlation energy (spin-conserved) =',EcRPA(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@URPAx correlation energy (spin-flip) =',EcRPA(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@URPAx correlation energy =',sum(EcRPA),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@URPAx total energy =',ENuc + EUHF + sum(EcRPA),' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end if
if(dotest) then
call dump_test_value('U','phURPAx correlation energy',sum(EcRPA))
end if
end subroutine

24
src/RPA/ppGRPA.f90
View File

@ -1,4 +1,4 @@
subroutine ppGRPA(TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
subroutine ppGRPA(dotest,TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
! Perform ppGRPA calculation
@ -7,6 +7,8 @@ subroutine ppGRPA(TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: TDA
logical,intent(in) :: doACFDT
integer,intent(in) :: nBas
@ -40,9 +42,9 @@ subroutine ppGRPA(TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
! Hello world
write(*,*)
write(*,*)'****************************************'
write(*,*)'| particle-particle GRPA calculation |'
write(*,*)'****************************************'
write(*,*)'**********************************'
write(*,*)'* Generalized pp-RPA Calculation *'
write(*,*)'**********************************'
write(*,*)
! Initialization
@ -70,8 +72,8 @@ subroutine ppGRPA(TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy =',EcRPA
write(*,'(2X,A50,F20.10)') 'Tr@ppRPA total energy =',ENuc + EHF + EcRPA
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(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -88,13 +90,17 @@ subroutine ppGRPA(TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
! 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 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)
end if
end subroutine

34
src/RPA/ppRPA.f90 → src/RPA/ppRRPA.f90
View File

@ -1,4 +1,4 @@
subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
subroutine ppRRPA(dotest,TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e)
! Perform ppRPA calculation
@ -7,6 +7,8 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipol
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: TDA
logical,intent(in) :: doACFDT
logical,intent(in) :: singlet
@ -42,9 +44,9 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipol
! Hello world
write(*,*)
write(*,*)'****************************************'
write(*,*)'| particle-particle RPA calculation |'
write(*,*)'****************************************'
write(*,*)'*********************************'
write(*,*)'* Restricted pp-RPA Calculation *'
write(*,*)'*********************************'
write(*,*)
! Initialization
@ -115,12 +117,14 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipol
endif
EcRPA(2) = 3d0*EcRPA(2)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (triplet) =',3d0*EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy =',EcRPA(1) + 3d0*EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@ppRPA total energy =',ENuc + EHF + EcRPA(1) + 3d0*EcRPA(2)
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppRRPA correlation energy (singlet) = ',EcRPA(1),'au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppRRPA correlation energy (triplet) = ',EcRPA(2),'au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppRRPA correlation energy = ',sum(EcRPA),'au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppRRPA total energy = ',ENuc + EHF + sum(EcRPA),'au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -137,13 +141,19 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipol
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
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(*,'(2X,A50,F20.10,A3)') 'AC@ppRRPA correlation energy (singlet) = ',EcRPA(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@ppRRPA correlation energy (triplet) = ',EcRPA(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@ppRRPA correlation energy = ',EcRPA(1) + EcRPA(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@ppRRPA total energy = ',ENuc + EHF + EcRPA(1) + EcRPA(2),' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end if
if(dotest) then
call dump_test_value('R','ppRRPA correlation energy',sum(EcRPA))
end if
end subroutine

36
src/RPA/ppURPA.f90
View File

@ -1,4 +1,4 @@
subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUHF,ERI_aaaa,ERI_aabb,ERI_bbbb,e)
subroutine ppURPA(dotest,TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUHF,ERI_aaaa,ERI_aabb,ERI_bbbb,e)
! Perform unrestricted pp-RPA calculations
@ -7,6 +7,8 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH
! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: TDA
logical,intent(in) :: doACFDT
logical,intent(in) :: spin_conserved
@ -35,20 +37,20 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH
double precision,allocatable :: X2sc(:,:),X2sf(:,:)
double precision,allocatable :: Y2sc(:,:),Y2sf(:,:)
double precision :: Ec_ppURPA(nspin)
double precision :: EcRPA(nspin)
double precision :: EcAC(nspin)
! Hello world
write(*,*)
write(*,*)'****************************************'
write(*,*)'| particle-particle URPA calculation |'
write(*,*)'****************************************'
write(*,*)'***********************************'
write(*,*)'* Unrestricted pp-RPA Calculation *'
write(*,*)'***********************************'
write(*,*)
! Initialization
Ec_ppURPA(:) = 0d0
EcRPA(:) = 0d0
EcAC(:) = 0d0
!alpha-beta block
@ -70,7 +72,7 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, &
nP_sc,nHaa,nHab,nHbb,nH_sc,1d0,e,ERI_aaaa, &
ERI_aabb,ERI_bbbb,Om1sc,X1sc,Y1sc, &
Om2sc,X2sc,Y2sc,Ec_ppURPA(ispin))
Om2sc,X2sc,Y2sc,EcRPA(ispin))
call print_excitation_energies('ppRPA@UHF (N+2)',iblock,nP_sc,Om1sc)
call print_excitation_energies('ppRPA@UHF (N-2)',iblock,nH_sc,Om2sc)
@ -96,7 +98,7 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, &
nP_sf,nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa, &
ERI_aabb,ERI_bbbb,Om1sf,X1sf,Y1sf, &
Om2sf,X2sf,Y2sf,Ec_ppURPA(ispin))
Om2sf,X2sf,Y2sf,EcRPA(ispin))
call print_excitation_energies('ppRPA@UHF (N+2)',iblock,nP_sf,Om1sf)
call print_excitation_energies('ppRPA@UHF (N-2)',iblock,nH_sf,Om2sf)
@ -116,17 +118,19 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,&
nP_sf,nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa,&
ERI_aabb,ERI_bbbb,Om1sf,X1sf,Y1sf,&
Om2sf,X2sf,Y2sf,Ec_ppURPA(ispin))
Om2sf,X2sf,Y2sf,EcRPA(ispin))
call print_excitation_energies('ppRPA@UHF (N+2)',iblock,nP_sf,Om1sf)
call print_excitation_energies('ppRPA@UHF (N-2)',iblock,nH_sf,Om2sf)
EcRPA(2) = 3d0*EcRPA(2)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (spin-conserved) =',Ec_ppURPA(1)
write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (spin-flip) =',3d0*Ec_ppURPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy =',Ec_ppURPA(1) + 3d0*Ec_ppURPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@ppRPA total energy =',ENuc + EUHF + Ec_ppURPA(1) + 3d0*Ec_ppURPA(2)
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppURPA correlation energy (spin-conserved) = ',EcRPA(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppURPA correlation energy (spin-flip) = ',EcRPA(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppURPA correlation energy = ',sum(EcRPA),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppURPA total energy = ',ENuc + EUHF + sum(EcRPA),' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -152,4 +156,10 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH
! end if
if(dotest) then
call dump_test_value('U','ppURPA correlation energy',sum(EcRPA))
end if
end subroutine

17
src/test/run_test.f90
View File

@ -10,6 +10,8 @@ subroutine run_test(doRtest,doUtest,doGtest)
! Local variables
double precision :: start_test ,end_test ,t_test
! Output variables
if(doRtest) then
@ -19,7 +21,12 @@ subroutine run_test(doRtest,doUtest,doGtest)
write(*,*) '****************************************'
write(*,*)
call wall_time(start_test)
call check_test_value('R')
call wall_time(end_test)
t_test = end_test - start_test
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for test of restricted branch = ',t_test,' seconds'
write(*,*)
write(*,*) '**************************'
@ -36,7 +43,12 @@ subroutine run_test(doRtest,doUtest,doGtest)
write(*,*) '******************************************'
write(*,*)
call wall_time(start_test)
call check_test_value('U')
call wall_time(end_test)
t_test = end_test - start_test
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for test of unrestricted branch = ',t_test,' seconds'
write(*,*)
write(*,*) '****************************'
@ -53,7 +65,12 @@ subroutine run_test(doRtest,doUtest,doGtest)
write(*,*) '*****************************************'
write(*,*)
call wall_time(start_test)
call check_test_value('G')
call wall_time(end_test)
t_test = end_test - start_test
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for test of generalized branch = ',t_test,' seconds'
write(*,*)
write(*,*) '***************************'

10
test/Gtest_ref.dat
View File

@ -1,2 +1,10 @@
# GHF energy
GHF energy
-85.160473883160876
GMP2 correlation energy
-0.128988144318866
phGRPA corrlation energy
-0.138552809810392
phGRPAx correlation energy
-0.368057033788489
ppGRPA correlation energy
-0.092561239023951

8
test/Rtest_ref.dat
View File

@ -2,3 +2,11 @@
-85.160473883160876
ROHF energy
-85.160473714509976
RMP2 correlation energy
-0.128988144386404
phRRPA correlation energy
-0.138552809856833
phRRPAx correlation energy
-0.197284981952336
ppRRPA correlation energy
-0.092561239071529

10
test/Utest_ref.dat
View File

@ -1,2 +1,10 @@
# UHF energy
UHF energy
-85.160473883160819
UMP2 correlation energy
-0.128988144318865
phURPA correlation energy
-0.138552809810790
phURPAx correlation energy
-0.197284981858218
ppURPA correlation energy
-0.103998858975444

4
test/methods.test
View File

@ -1,7 +1,7 @@
# RHF UHF GHF ROHF
T T T T
# MP2 MP3
F F
T T
# 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 F T
# G0F2 evGF2 qsGF2 G0F3 evGF3
F F F F F
# G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW