LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-01-03 01:56:09 +01:00
Code Releases Activity

test in GT

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-11-13 18:36:09 +01:00
parent 976d2b0c54
commit 594805d429
15 changed files with 244 additions and 125 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
input/methods
View File

@ -1,5 +1,5 @@
# RHF UHF GHF ROHF # RHF UHF GHF ROHF
T F T F T F F F
# MP2 MP3 # MP2 MP3
T T T T
# CCD pCCD DCD CCSD CCSD(T) # CCD pCCD DCD CCSD CCSD(T)
@ -13,8 +13,8 @@
# G0F2 evGF2 qsGF2 G0F3 evGF3 # G0F2 evGF2 qsGF2 G0F3 evGF3
T T F F F T T F F F
# G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW # G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW
F F F F F F T T T F F F
# G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
F F F F F F F F F F F F
# Rtest Utest Gtest # Rtest Utest Gtest
T F T T F F

4
input/options
View File

@ -9,9 +9,9 @@
# GF: maxSCF thresh DIIS lin eta renorm reg # GF: maxSCF thresh DIIS lin eta renorm reg
256 0.00001 5 F 0.0 0 F 256 0.00001 5 F 0.0 0 F
# GW: maxSCF thresh DIIS lin eta TDA_W reg # GW: maxSCF thresh DIIS lin eta TDA_W reg
10 0.00001 5 F 0.0 F F 256 0.00001 5 F 0.0 F F
# GT: maxSCF thresh DIIS lin eta TDA_T reg # GT: maxSCF thresh DIIS lin eta TDA_T reg
256 0.00001 5 F 0.0 F F 256 0.00001 5 F 0.0 F F
# ACFDT: AC Kx XBS # ACFDT: AC Kx XBS
F F T F F T
# BSE: phBSE phBSE2 ppBSE dBSE dTDA # BSE: phBSE phBSE2 ppBSE dBSE dTDA

23
src/GT/G0T0eh.f90 → src/GT/RG0T0eh.f90
View File

@ -1,6 +1,5 @@
subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & subroutine RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
ERI,dipole_int,eHF)
! Perform ehG0T0 calculation ! Perform ehG0T0 calculation
@ -10,6 +9,8 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
! Input variables ! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: doACFDT logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel logical,intent(in) :: exchange_kernel
logical,intent(in) :: doXBS logical,intent(in) :: doXBS
@ -73,9 +74,9 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
! Hello world ! Hello world
write(*,*) write(*,*)
write(*,*)'************************************************' write(*,*)'*********************************'
write(*,*)'| One-shot G0T0eh calculation |' write(*,*)'* Restricted G0T0eh Calculation *'
write(*,*)'************************************************' write(*,*)'*********************************'
write(*,*) write(*,*)
! Initialization ! Initialization
@ -167,4 +168,14 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
call print_G0T0eh(nBas,nO,eHF,ENuc,ERHF,Sig,Z,eGT,EcRPA,EcGM) call print_G0T0eh(nBas,nO,eHF,ENuc,ERHF,Sig,Z,eGT,EcRPA,EcGM)
! Testing zone
if(dotest) then
call dump_test_value('R','RG0T0eh correlation energy',EcRPA)
call dump_test_value('R','RG0T0eh HOMO energy',eGT(nO))
call dump_test_value('R','RG0T0eh LUMO energy',eGT(nO+1))
end if
end subroutine end subroutine

34
src/GT/G0T0pp.f90 → src/GT/RG0T0pp.f90
View File

@ -1,5 +1,5 @@
subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform one-shot calculation with a T-matrix self-energy (G0T0) ! Perform one-shot calculation with a T-matrix self-energy (G0T0)
@ -8,6 +8,8 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
! Input variables ! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: doACFDT logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel logical,intent(in) :: exchange_kernel
logical,intent(in) :: doXBS logical,intent(in) :: doXBS
@ -65,9 +67,9 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
! Hello world ! Hello world
write(*,*) write(*,*)
write(*,*)'************************************************' write(*,*)'*********************************'
write(*,*)'| One-shot G0T0pp calculation |' write(*,*)'* Restricted G0T0pp Calculation *'
write(*,*)'************************************************' write(*,*)'*********************************'
write(*,*) write(*,*)
@ -260,8 +262,8 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp correlation energy (singlet) =',EcBSE(1),' au' write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp correlation energy (singlet) =',EcBSE(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp correlation energy (triplet) =',EcBSE(2),' au' write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp correlation energy (triplet) =',EcBSE(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp correlation energy =',EcBSE(1) + EcBSE(2),' au' write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp correlation energy =',sum(EcBSE),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2),' au' write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp total energy =',ENuc + ERHF + sum(EcBSE),' au'
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,*) write(*,*)
@ -296,8 +298,8 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp correlation energy (singlet) =',EcBSE(1),' au' write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp correlation energy (singlet) =',EcBSE(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp correlation energy (triplet) =',EcBSE(2),' au' write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp correlation energy (triplet) =',EcBSE(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp correlation energy =',EcBSE(1) + EcBSE(2),' au' write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp correlation energy =',sum(EcBSE),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2),' au' write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp total energy =',ENuc + ERHF + sum(EcBSE),' au'
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,*) write(*,*)
@ -315,11 +317,21 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp correlation energy (singlet) =',EcBSE(1),' au' write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp correlation energy (singlet) =',EcBSE(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp correlation energy (triplet) =',EcBSE(2),' au' write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp correlation energy (triplet) =',EcBSE(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp correlation energy =',EcBSE(1) + EcBSE(2),' au' write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp correlation energy =',sum(EcBSE),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2),' au' write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp total energy =',ENuc + ERHF + sum(EcBSE),' au'
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,*) write(*,*)
end if end if
! Testing zone
if(dotest) then
call dump_test_value('R','RG0T0pp correlation energy',sum(EcRPA))
call dump_test_value('R','RG0T0pp HOMO energy',eGT(nO))
call dump_test_value('R','RG0T0pp LUMO energy',eGT(nO+1))
end if
end subroutine end subroutine

42
src/GT/RGT.f90
View File

@ -1,4 +1,4 @@
subroutine RGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thresh,max_diis,doACFDT, & subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thresh,max_diis,doACFDT, &
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, & linearize,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) ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF)
@ -10,12 +10,14 @@ subroutine RGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
! Input variables ! Input variables
logical :: doG0T0pp logical,intent(in) :: dotest
logical :: doevGTpp
logical :: doqsGTpp logical,intent(in) :: doG0T0pp
logical :: doG0T0eh logical,intent(in) :: doevGTpp
logical :: doevGTeh logical,intent(in) :: doqsGTpp
logical :: doqsGTeh logical,intent(in) :: doG0T0eh
logical,intent(in) :: doevGTeh
logical,intent(in) :: doqsGTeh
integer,intent(in) :: maxSCF integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis integer,intent(in) :: max_diis
@ -74,7 +76,7 @@ subroutine RGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
if(doG0T0pp) then if(doG0T0pp) then
call wall_time(start_GT) call wall_time(start_GT)
call G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & call RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_GT) call wall_time(end_GT)
@ -91,8 +93,8 @@ subroutine RGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
if(doevGTpp) then if(doevGTpp) then
call wall_time(start_GT) call wall_time(start_GT)
call evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & call evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_GT) call wall_time(end_GT)
t_GT = end_GT - start_GT t_GT = end_GT - start_GT
@ -108,9 +110,9 @@ subroutine RGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
if(doqsGTpp) then if(doqsGTpp) then
call wall_time(start_GT) call wall_time(start_GT)
call qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & call qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,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, & 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) PHF,cHF,epsHF)
call wall_time(end_GT) call wall_time(end_GT)
t_GT = end_GT - start_GT t_GT = end_GT - start_GT
@ -126,8 +128,8 @@ subroutine RGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
if(doG0T0eh) then if(doG0T0eh) then
call wall_time(start_GT) call wall_time(start_GT)
call G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & call RG0T0eh(dotest,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) linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_GT) call wall_time(end_GT)
t_GT = end_GT - start_GT t_GT = end_GT - start_GT
@ -143,8 +145,8 @@ subroutine RGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
if(doevGTeh) then if(doevGTeh) then
call wall_time(start_GT) call wall_time(start_GT)
call evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & call evRGTeh(dotest,maxSCF,thresh,max_diis,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) singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_GT) call wall_time(end_GT)
t_GT = end_GT - start_GT t_GT = end_GT - start_GT
@ -160,9 +162,9 @@ subroutine RGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
if(doqsGTeh) then if(doqsGTeh) then
call wall_time(start_GT) call wall_time(start_GT)
call qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & call qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,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, & 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) PHF,cHF,epsHF)
call wall_time(end_GT) call wall_time(end_GT)
t_GT = end_GT - start_GT t_GT = end_GT - start_GT

37
src/GT/UG0T0pp.f90
View File

@ -1,4 +1,4 @@
subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & subroutine UG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV, & spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV, &
nR,nS,ENuc,EUHF,ERI_aaaa,ERI_aabb,ERI_bbbb, & nR,nS,ENuc,EUHF,ERI_aaaa,ERI_aabb,ERI_bbbb, &
dipole_int_aa,dipole_int_bb,cHF,eHF) dipole_int_aa,dipole_int_bb,cHF,eHF)
@ -10,6 +10,8 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
! Input variables ! Input variables
logical,intent(in) :: dotest
logical,intent(in) :: doACFDT logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel logical,intent(in) :: exchange_kernel
logical,intent(in) :: doXBS logical,intent(in) :: doXBS
@ -60,17 +62,16 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
double precision,allocatable :: rho2ab(:,:,:),rho2aa(:,:,:),rho2bb(:,:,:) double precision,allocatable :: rho2ab(:,:,:),rho2aa(:,:,:),rho2bb(:,:,:)
double precision,allocatable :: SigT(:,:) double precision,allocatable :: SigT(:,:)
double precision,allocatable :: Z(:,:) double precision,allocatable :: Z(:,:)
double precision,allocatable :: eG0T0(:,:) double precision,allocatable :: eGT(:,:)
! Output variables ! Output variables
! Hello world ! Hello world
write(*,*) write(*,*)
write(*,*)'************************************************' write(*,*)'***********************************'
write(*,*)'| One-shot G0T0 calculation |' write(*,*)'* Unrestricted G0T0pp Calculation *'
write(*,*)'| *** Unrestricted version *** |' write(*,*)'***********************************'
write(*,*)'************************************************'
write(*,*) write(*,*)
! Dimensions of the pp-URPA linear reponse matrices ! Dimensions of the pp-URPA linear reponse matrices
@ -101,7 +102,7 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
Om1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), & Om1bb(nPbb),X1bb(nPbb,nPbb),Y1bb(nHbb,nPbb), &
Om2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), & Om2bb(nPbb),X2bb(nPbb,nPbb),Y2bb(nHbb,nPbb), &
rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb), & rho1bb(nBas,nBas,nPbb),rho2bb(nBas,nBas,nHbb), &
SigT(nBas,nspin),Z(nBas,nspin),eG0T0(nBas,nspin)) SigT(nBas,nspin),Z(nBas,nspin),eGT(nBas,nspin))
!---------------------------------------------- !----------------------------------------------
! alpha-beta block ! alpha-beta block
@ -183,7 +184,7 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
if(linearize) then if(linearize) then
eG0T0(:,:) = eHF(:,:) + Z(:,:)*SigT(:,:) eGT(:,:) = eHF(:,:) + Z(:,:)*SigT(:,:)
else else
@ -203,7 +204,7 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
ispin = 1 ispin = 1
iblock = 3 iblock = 3
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPab,nHaa,nHab,nHbb,nHab,1d0,eG0T0,ERI_aaaa, & call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPab,nHaa,nHab,nHbb,nHab,1d0,eGT,ERI_aaaa, &
ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin)) ERI_aabb,ERI_bbbb,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin))
!alpha-alpha block !alpha-alpha block
@ -211,7 +212,7 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
ispin = 2 ispin = 2
iblock = 4 iblock = 4
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPaa,nHaa,nHab,nHbb,nHaa,1d0,eG0T0,ERI_aaaa, & call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPaa,nHaa,nHab,nHbb,nHaa,1d0,eGT,ERI_aaaa, &
ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin)) ERI_aabb,ERI_bbbb,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin))
Ecaa = EcRPA(2) Ecaa = EcRPA(2)
@ -220,7 +221,7 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
iblock = 7 iblock = 7
call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPbb,nHaa,nHab,nHbb,nHbb,1d0,eG0T0,ERI_aaaa, & call ppULR(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPbb,nHaa,nHab,nHbb,nHbb,1d0,eGT,ERI_aaaa, &
ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,EcRPA(ispin)) ERI_aabb,ERI_bbbb,Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,EcRPA(ispin))
Ecbb = EcRPA(2) Ecbb = EcRPA(2)
@ -228,7 +229,7 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
EcRPA(1) = EcRPA(1) - EcRPA(2) EcRPA(1) = EcRPA(1) - EcRPA(2)
EcRPA(2) = 3d0*EcRPA(2) EcRPA(2) = 3d0*EcRPA(2)
call print_UG0T0(nBas,nO,eHF,ENuc,EUHF,SigT,Z,eG0T0,EcGM,EcRPA) call print_UG0T0(nBas,nO,eHF,ENuc,EUHF,SigT,Z,eGT,EcGM,EcRPA)
! Free memory ! Free memory
@ -236,4 +237,16 @@ subroutine UG0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, &
Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, & Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, &
Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,rho1bb,rho2bb) Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,rho1bb,rho2bb)
! Testing zone
if(dotest) then
call dump_test_value('U','UG0T0pp correlation energy',sum(EcRPA))
call dump_test_value('U','UG0T0pp HOMOa energy',eGT(nO(1),1))
call dump_test_value('U','UG0T0pp LUMOa energy',eGT(nO(1)+1,1))
call dump_test_value('U','UG0T0pp HOMOa energy',eGT(nO(2),2))
call dump_test_value('U','UG0T0pp LUMOa energy',eGT(nO(2)+1,2))
end if
end subroutine end subroutine

22
src/GT/UGT.f90
View File

@ -1,4 +1,4 @@
subroutine UGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thresh,max_diis,doACFDT, & subroutine UGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thresh,max_diis,doACFDT, &
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, & exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, & linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, &
ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF)
@ -10,12 +10,14 @@ subroutine UGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
! Input variables ! Input variables
logical :: doG0T0pp logical,intent(in) :: dotest
logical :: doevGTpp
logical :: doqsGTpp logical,intent(in) :: doG0T0pp
logical :: doG0T0eh logical,intent(in) :: doevGTpp
logical :: doevGTeh logical,intent(in) :: doqsGTpp
logical :: doqsGTeh logical,intent(in) :: doG0T0eh
logical,intent(in) :: doevGTeh
logical,intent(in) :: doqsGTeh
integer,intent(in) :: maxSCF integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis integer,intent(in) :: max_diis
@ -77,7 +79,7 @@ subroutine UGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
if(doG0T0pp) then if(doG0T0pp) then
call wall_time(start_GT) call wall_time(start_GT)
call UG0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, & call UG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb, & linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb, &
dipole_int_aa,dipole_int_bb,cHF,epsHF) dipole_int_aa,dipole_int_bb,cHF,epsHF)
call wall_time(end_GT) call wall_time(end_GT)
@ -95,7 +97,7 @@ subroutine UGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
if(doevGTpp) then if(doevGTpp) then
call wall_time(start_GT) call wall_time(start_GT)
call evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, & call evUGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, & linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, &
cHF,epsHF) cHF,epsHF)
call wall_time(end_GT) call wall_time(end_GT)
@ -113,7 +115,7 @@ subroutine UGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thre
if(doqsGTpp) then if(doqsGTpp) then
call wall_time(start_GT) call wall_time(start_GT)
call qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, & call qsUGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
eta,regularize,nBas,nC,nO,nV,nR,nS,nNuc,ZNuc,rNuc,ENuc,EHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, & eta,regularize,nBas,nC,nO,nV,nR,nS,nNuc,ZNuc,rNuc,ENuc,EHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, &
dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF)
call wall_time(end_GT) call wall_time(end_GT)

22
src/GT/evGTeh.f90 → src/GT/evRGTeh.f90
View File

@ -1,5 +1,5 @@
subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & subroutine evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform self-consistent eigenvalue-only ehGT calculation ! Perform self-consistent eigenvalue-only ehGT calculation
@ -8,6 +8,8 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
! Input variables ! Input variables
logical,intent(in) :: dotest
integer,intent(in) :: maxSCF integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis integer,intent(in) :: max_diis
double precision,intent(in) :: thresh double precision,intent(in) :: thresh
@ -70,9 +72,9 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
! Hello world ! Hello world
write(*,*) write(*,*)
write(*,*)'************************************************' write(*,*)'**********************************'
write(*,*)'| Self-consistent evGTeh calculation |' write(*,*)'* Restricted evRGTeh Calculation *'
write(*,*)'************************************************' write(*,*)'**********************************'
write(*,*) write(*,*)
! TDA for T ! TDA for T
@ -296,4 +298,14 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
! end if ! end if
! Testing zone
if(dotest) then
call dump_test_value('R','evRGTeh correlation energy',EcRPA)
call dump_test_value('R','evRGTeh HOMO energy',eGT(nO))
call dump_test_value('R','evRGTeh LUMO energy',eGT(nO+1))
end if
end subroutine end subroutine

31
src/GT/evGTpp.f90 → src/GT/evRGTpp.f90
View File

@ -1,5 +1,5 @@
subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform eigenvalue self-consistent calculation with a T-matrix self-energy (evGT) ! Perform eigenvalue self-consistent calculation with a T-matrix self-energy (evGT)
@ -8,6 +8,8 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
! Input variables ! Input variables
logical,intent(in) :: dotest
integer,intent(in) :: maxSCF integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis integer,intent(in) :: max_diis
double precision,intent(in) :: thresh double precision,intent(in) :: thresh
@ -74,9 +76,9 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
! Hello world ! Hello world
write(*,*) write(*,*)
write(*,*)'************************************************' write(*,*)'*********************************'
write(*,*)'| Self-consistent evGTpp calculation |' write(*,*)'* Restricted evGTpp Calculation *'
write(*,*)'************************************************' write(*,*)'*********************************'
write(*,*) write(*,*)
! TDA for T ! TDA for T
@ -265,8 +267,8 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGTpp correlation energy (singlet) =',EcBSE(1) write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGTpp correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGTpp correlation energy (triplet) =',EcBSE(2) write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGTpp correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGTpp correlation energy =',EcBSE(1) + EcBSE(2) write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGTpp correlation energy =',sum(EcBSE)
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGTpp total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGTpp total energy =',ENuc + ERHF + sum(EcBSE)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,*) write(*,*)
@ -301,8 +303,8 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'AC@phBSE@evGTpp correlation energy (singlet) =',EcBSE(1) write(*,'(2X,A50,F20.10)') 'AC@phBSE@evGTpp correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A50,F20.10)') 'AC@phBSE@evGTpp correlation energy (triplet) =',EcBSE(2) write(*,'(2X,A50,F20.10)') 'AC@phBSE@evGTpp correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A50,F20.10)') 'AC@phBSE@evGTpp correlation energy =',EcBSE(1) + EcBSE(2) write(*,'(2X,A50,F20.10)') 'AC@phBSE@evGTpp correlation energy =',sum(EcBSE)
write(*,'(2X,A50,F20.10)') 'AC@phBSE@evGTpp total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) write(*,'(2X,A50,F20.10)') 'AC@phBSE@evGTpp total energy =',ENuc + ERHF + sum(EcBSE)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,*) write(*,*)
@ -310,4 +312,15 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T
end if end if
! Testing zone
if(dotest) then
call dump_test_value('R','evRGTpp correlation energy',sum(EcRPA))
call dump_test_value('R','evRGTpp HOMO energy',eGT(nO))
call dump_test_value('R','evRGTpp LUMO energy',eGT(nO+1))
end if
end subroutine end subroutine

23
src/GT/evUGTpp.f90
View File

@ -1,4 +1,4 @@
subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & subroutine evUGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,& TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,&
eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,ERI_aaaa, & eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,ERI_aaaa, &
ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF) ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF)
@ -9,6 +9,9 @@ subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
include 'parameters.h' include 'parameters.h'
! Input variables ! Input variables
logical,intent(in) :: dotest
integer,intent(in) :: maxSCF integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis integer,intent(in) :: max_diis
double precision,intent(in) :: thresh double precision,intent(in) :: thresh
@ -73,9 +76,9 @@ subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
! Hello world ! Hello world
write(*,*) write(*,*)
write(*,*)'************************************************' write(*,*)'***********************************'
write(*,*)'| Self-consistent evUGT calculation |' write(*,*)'* Unrestricted evGTpp Calculation *'
write(*,*)'************************************************' write(*,*)'***********************************'
write(*,*) write(*,*)
! Dimensions of the pp-URPA linear reponse matrices ! Dimensions of the pp-URPA linear reponse matrices
@ -284,4 +287,16 @@ subroutine evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, & Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, &
Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,rho1bb,rho2bb) Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,rho1bb,rho2bb)
! Testing zone
if(dotest) then
call dump_test_value('U','evUGTpp correlation energy',sum(EcRPA))
call dump_test_value('U','evUGTpp HOMOa energy',eGT(nO(1),1))
call dump_test_value('U','evUGTpp LUMOa energy',eGT(nO(1)+1,1))
call dump_test_value('U','evUGTpp HOMOa energy',eGT(nO(2),2))
call dump_test_value('U','evUGTpp LUMOa energy',eGT(nO(2)+1,2))
end if
end subroutine end subroutine

60
src/GT/qsGTeh.f90 → src/GT/qsRGTeh.f90
View File

@ -1,6 +1,6 @@
subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA, & subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA, &
dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
! Perform a quasiparticle self-consistent GTeh calculation ! Perform a quasiparticle self-consistent GTeh calculation
@ -9,6 +9,8 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
! Input variables ! Input variables
logical,intent(in) :: dotest
integer,intent(in) :: maxSCF integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis integer,intent(in) :: max_diis
double precision,intent(in) :: thresh double precision,intent(in) :: thresh
@ -73,7 +75,7 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
double precision :: dipole(ncart) double precision :: dipole(ncart)
logical :: print_T = .true. logical :: print_T = .true.
double precision,allocatable :: error_diis(:,:) double precision,allocatable :: err_diis(:,:)
double precision,allocatable :: F_diis(:,:) double precision,allocatable :: F_diis(:,:)
double precision,allocatable :: Aph(:,:) double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:) double precision,allocatable :: Bph(:,:)
@ -94,14 +96,14 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
double precision,allocatable :: Sig(:,:) double precision,allocatable :: Sig(:,:)
double precision,allocatable :: Sigp(:,:) double precision,allocatable :: Sigp(:,:)
double precision,allocatable :: Z(:) double precision,allocatable :: Z(:)
double precision,allocatable :: error(:,:) double precision,allocatable :: err(:,:)
! Hello world ! Hello world
write(*,*) write(*,*)
write(*,*)'************************************************' write(*,*)'*********************************'
write(*,*)'| Self-consistent qsGTeh calculation |' write(*,*)'* Restricted qsGTeh Calculation *'
write(*,*)'************************************************' write(*,*)'*********************************'
write(*,*) write(*,*)
! Warning ! Warning
@ -131,21 +133,21 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
allocate(Aph(nS,nS),Bph(nS,nS),eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & allocate(Aph(nS,nS),Bph(nS,nS),eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), & J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), &
rhoL(nBas,nBas,nS),rhoR(nBas,nBas,nS),error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) rhoL(nBas,nBas,nS),rhoR(nBas,nBas,nS),err(nBas,nBas),err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Initialization ! Initialization
nSCF = -1 nSCF = -1
n_diis = 0 n_diis = 0
ispin = 2 ispin = 2
Conv = 1d0 Conv = 1d0
P(:,:) = PHF(:,:) P(:,:) = PHF(:,:)
eGT(:) = eHF(:) eGT(:) = eHF(:)
eOld(:) = eHF(:) eOld(:) = eHF(:)
c(:,:) = cHF(:,:) c(:,:) = cHF(:,:)
F_diis(:,:) = 0d0 F_diis(:,:) = 0d0
error_diis(:,:) = 0d0 err_diis(:,:) = 0d0
rcond = 0d0 rcond = 0d0
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Main loop ! Main loop
@ -198,14 +200,14 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
! Compute commutator and convergence criteria ! Compute commutator and convergence criteria
error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
! DIIS extrapolation ! DIIS extrapolation
if(max_diis > 1) then if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis) n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis,F_diis,err,F)
end if end if
@ -223,7 +225,7 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
! Save quasiparticles energy for next cycle ! Save quasiparticles energy for next cycle
Conv = maxval(abs(error)) Conv = maxval(abs(err))
eOld(:) = eGT(:) eOld(:) = eGT(:)
!------------------------------------------------------------------------ !------------------------------------------------------------------------
@ -276,7 +278,7 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
! Deallocate memory ! Deallocate memory
deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Z,Om,XpY,XmY,rhoL,rhoR,error,error_diis,F_diis) deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Z,Om,XpY,XmY,rhoL,rhoR,err,err_diis,F_diis)
! Perform BSE calculation ! Perform BSE calculation
@ -332,4 +334,14 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
! end if ! end if
! Testing zone
if(dotest) then
call dump_test_value('R','qsRGTeh correlation energy',EcRPA)
call dump_test_value('R','qsRGTeh HOMO energy',eGT(nO))
call dump_test_value('R','qsRGTeh LUMO energy',eGT(nO+1))
end if
end subroutine end subroutine

32
src/GT/qsGTpp.f90 → src/GT/qsRGTpp.f90
View File

@ -1,6 +1,6 @@
subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA, & subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA, &
dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
! Perform a quasiparticle self-consistent GT calculation ! Perform a quasiparticle self-consistent GT calculation
@ -9,6 +9,7 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T
! Input variables ! Input variables
logical,intent(in) :: dotest
integer,intent(in) :: maxSCF integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis integer,intent(in) :: max_diis
double precision,intent(in) :: thresh double precision,intent(in) :: thresh
@ -98,9 +99,9 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T
! Hello world ! Hello world
write(*,*) write(*,*)
write(*,*)'************************************************' write(*,*)'*********************************'
write(*,*)'| Self-consistent qsGTpp calculation |' write(*,*)'* Restricted qsGTpp Calculation *'
write(*,*)'************************************************' write(*,*)'*********************************'
write(*,*) write(*,*)
! Dimensions of the pp-RPA linear reponse matrices ! Dimensions of the pp-RPA linear reponse matrices
@ -341,8 +342,8 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGTpp correlation energy (singlet) =',EcBSE(1) write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGTpp correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A50,F20.10)') 'Tr@BphSE@qsGTpp correlation energy (triplet) =',EcBSE(2) write(*,'(2X,A50,F20.10)') 'Tr@BphSE@qsGTpp correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGTpp correlation energy =',EcBSE(1) + EcBSE(2) write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGTpp correlation energy =',sum(EcBSE)
write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGTpp total energy =',ENuc + EqsGT + EcBSE(1) + EcBSE(2) write(*,'(2X,A50,F20.10)') 'Tr@phBSE@qsGTpp total energy =',ENuc + EqsGT + sum(EcBSE)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,*) write(*,*)
@ -370,8 +371,8 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'AC@phBSE@qsGTpp correlation energy (singlet) =',EcBSE(1) write(*,'(2X,A50,F20.10)') 'AC@phBSE@qsGTpp correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A50,F20.10)') 'AC@phBSE@qsGTpp correlation energy (triplet) =',EcBSE(2) write(*,'(2X,A50,F20.10)') 'AC@phBSE@qsGTpp correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A50,F20.10)') 'AC@phBSE@qsGTpp correlation energy =',EcBSE(1) + EcBSE(2) write(*,'(2X,A50,F20.10)') 'AC@phBSE@qsGTpp correlation energy =',sum(EcBSE)
write(*,'(2X,A50,F20.10)') 'AC@phBSE@qsGTpp total energy =',ENuc + EqsGT + EcBSE(1) + EcBSE(2) write(*,'(2X,A50,F20.10)') 'AC@phBSE@qsGTpp total energy =',ENuc + EqsGT + sum(EcBSE)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,*) write(*,*)
@ -379,4 +380,15 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T
end if end if
! Testing zone
if(dotest) then
call dump_test_value('R','qsRGTpp correlation energy',sum(EcRPA))
call dump_test_value('R','qsRGTpp HOMO energy',eGT(nO))
call dump_test_value('R','qsRGTpp LUMO energy',eGT(nO+1))
end if
end subroutine end subroutine

23
src/GT/qsUGTpp.f90
View File

@ -1,4 +1,4 @@
subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & subroutine qsUGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,& TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,&
eta,regularize,nBas,nC,nO,nV,nR,nS,nNuc,ZNuc,rNuc,ENuc,EUHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,& eta,regularize,nBas,nC,nO,nV,nR,nS,nNuc,ZNuc,rNuc,ENuc,EUHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,&
ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF)
@ -9,6 +9,9 @@ subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
include 'parameters.h' include 'parameters.h'
! Input variables ! Input variables
logical,intent(in) :: dotest
integer,intent(in) :: maxSCF integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis integer,intent(in) :: max_diis
double precision,intent(in) :: thresh double precision,intent(in) :: thresh
@ -103,9 +106,9 @@ subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
! Hello world ! Hello world
write(*,*) write(*,*)
write(*,*)'************************************************' write(*,*)'***********************************'
write(*,*)'| Self-consistent qsUGT calculation |' write(*,*)'* Unrestricted evGTpp Calculation *'
write(*,*)'************************************************' write(*,*)'***********************************'
write(*,*) write(*,*)
! Dimensions of the pp-URPA linear reponse matrices ! Dimensions of the pp-URPA linear reponse matrices
@ -409,4 +412,16 @@ write(*,*) 'EcGM', EcGM(1)
deallocate(c,cp,P,F,Fp,J,K,SigT,SigTp,Z,error,error_diis,F_diis) deallocate(c,cp,P,F,Fp,J,K,SigT,SigTp,Z,error,error_diis,F_diis)
! Testing zone
if(dotest) then
call dump_test_value('U','qsUGTpp correlation energy',sum(EcRPA))
call dump_test_value('U','qsUGTpp HOMOa energy',eGT(nO(1),1))
call dump_test_value('U','qsUGTpp LUMOa energy',eGT(nO(1)+1,1))
call dump_test_value('U','qsUGTpp HOMOa energy',eGT(nO(2),2))
call dump_test_value('U','qsUGTpp LUMOa energy',eGT(nO(2)+1,2))
end if
end subroutine end subroutine

6
test/methods.test
View File

@ -1,5 +1,5 @@
# RHF UHF GHF ROHF # RHF UHF GHF ROHF
T F T F T F F F
# MP2 MP3 # MP2 MP3
T T T T
# CCD pCCD DCD CCSD CCSD(T) # CCD pCCD DCD CCSD CCSD(T)
@ -13,8 +13,8 @@
# G0F2 evGF2 qsGF2 G0F3 evGF3 # G0F2 evGF2 qsGF2 G0F3 evGF3
T T F F F T T F F F
# G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW # G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW
F F F F F F T T T F F F
# G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
F F F F F F F F F F F F
# Rtest Utest Gtest # Rtest Utest Gtest
T F T T F F

4
test/options.test
View File

@ -9,9 +9,9 @@
# GF: maxSCF thresh DIIS lin eta renorm reg # GF: maxSCF thresh DIIS lin eta renorm reg
256 0.00001 5 F 0.0 0 F 256 0.00001 5 F 0.0 0 F
# GW: maxSCF thresh DIIS lin eta TDA_W reg # GW: maxSCF thresh DIIS lin eta TDA_W reg
10 0.00001 5 F 0.0 F F 256 0.00001 5 F 0.0 F F
# GT: maxSCF thresh DIIS lin eta TDA_T reg # GT: maxSCF thresh DIIS lin eta TDA_T reg
256 0.00001 5 F 0.0 F F 256 0.00001 5 F 0.0 F F
# ACFDT: AC Kx XBS # ACFDT: AC Kx XBS
F F T F F T
# BSE: phBSE phBSE2 ppBSE dBSE dTDA # BSE: phBSE phBSE2 ppBSE dBSE dTDA