LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 12:23:50 +01:00
Code Releases Activity

more cleaning

Browse Source
This commit is contained in:
Pierre-Francois Loos 2024-09-11 10:13:11 +02:00
parent f04ea269da
commit c1107474ff
14 changed files with 305 additions and 310 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
src/GF/RG0F2.f90
View File

@ -1,5 +1,5 @@
subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform a one-shot second-order Green function calculation
@ -20,7 +20,9 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,
logical,intent(in) :: linearize
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nO
integer,intent(in) :: nC
integer,intent(in) :: nV
@ -28,9 +30,9 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart)
! Local variables
@ -51,17 +53,17 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,
! Memory allocation
allocate(SigC(nBas), Z(nBas), eGFlin(nBas), eGF(nBas))
allocate(SigC(nOrb), Z(nOrb), eGFlin(nOrb), eGF(nOrb))
! Frequency-dependent second-order contribution
if(regularize) then
call RGF2_reg_self_energy_diag(eta,nBas,nC,nO,nV,nR,eHF,ERI,SigC,Z)
call RGF2_reg_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eHF,ERI,SigC,Z)
else
call RGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,eHF,ERI,SigC,Z)
call RGF2_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eHF,ERI,SigC,Z)
end if
@ -78,20 +80,20 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,
write(*,*) ' *** Quasiparticle energies obtained by root search *** '
write(*,*)
call RGF2_QP_graph(eta,nBas,nC,nO,nV,nR,eHF,ERI,eGFlin,eHF,eGF,Z)
call RGF2_QP_graph(eta,nOrb,nC,nO,nV,nR,eHF,ERI,eGFlin,eHF,eGF,Z)
end if
! Print results
call RMP2(.false.,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
call print_RG0F2(nBas,nO,eHF,SigC,eGF,Z,ENuc,ERHF,Ec)
call RMP2(.false.,regularize,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
call print_RG0F2(nOrb,nO,eHF,SigC,eGF,Z,ENuc,ERHF,Ec)
! Perform BSE@GF2 calculation
if(dophBSE) then
call RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
call RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
@ -108,7 +110,7 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,
if(doppBSE) then
call RGF2_ppBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
call RGF2_ppBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
EcBSE(2) = 3d0*EcBSE(2)

120
src/GF/RG0F3.f90
View File

@ -1,4 +1,4 @@
subroutine RG0F3(dotest,renormalization,nBas,nC,nO,nV,nR,V,e0)
subroutine RG0F3(dotest,renormalization,nBas,nOrb,nC,nO,nV,nR,V,e0)
! Perform third-order Green function calculation in diagonal approximation
@ -9,8 +9,8 @@
logical,intent(in) :: dotest
integer,intent(in) :: renormalization
integer,intent(in) :: nBas,nC,nO,nV,nR
double precision,intent(in) :: e0(nBas),V(nBas,nBas,nBas,nBas)
integer,intent(in) :: nBas,nOrb,nC,nO,nV,nR
double precision,intent(in) :: e0(nOrb),V(nOrb,nOrb,nOrb,nOrb)
! Local variables
@ -31,9 +31,9 @@
! Memory allocation
allocate(eGF3(nBas),Sig2(nBas),SigInf(nBas),Sig3(nBas), &
App(nBas,6),Bpp(nBas,2),Cpp(nBas,6),Dpp(nBas,6), &
Z(nBas),X2h1p(nBas),X1h2p(nBas),Sig2h1p(nBas),Sig1h2p(nBas))
allocate(eGF3(nOrb),Sig2(nOrb),SigInf(nOrb),Sig3(nOrb), &
App(nOrb,6),Bpp(nOrb,2),Cpp(nOrb,6),Dpp(nOrb,6), &
Z(nOrb),X2h1p(nOrb),X1h2p(nOrb),Sig2h1p(nOrb),Sig1h2p(nOrb))
!------------------------------------------------------------------------
! Compute third-order frequency-independent contribution
@ -41,12 +41,12 @@
App(:,:) = 0d0
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
eps1 = e0(j) + e0(i) - e0(a) - e0(b)
eps2 = e0(k) + e0(i) - e0(a) - e0(b)
@ -61,12 +61,12 @@
end do
end do
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do c=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
do c=nO+1,nOrb-nR
eps1 = e0(j) + e0(i) - e0(a) - e0(b)
eps2 = e0(j) + e0(i) - e0(a) - e0(c)
@ -81,12 +81,12 @@
end do
end do
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do c=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
do c=nO+1,nOrb-nR
eps1 = e0(j) + e0(i) - e0(a) - e0(b)
eps2 = e0(j) - e0(c)
@ -103,12 +103,12 @@
App(:,4) = App(:,3)
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
eps1 = e0(j) + e0(i) - e0(a) - e0(b)
eps2 = e0(k) - e0(b)
@ -133,10 +133,10 @@
Bpp(:,:) = 0d0
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do a=nO+1,nOrb-nR
eps = eGF3(p) + e0(a) - e0(i) - e0(j)
@ -148,10 +148,10 @@
end do
end do
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
eps = eGF3(p) + e0(i) - e0(a) - e0(b)
@ -171,12 +171,12 @@
Cpp(:,:) = 0d0
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do c=nO+1,nBas-nR
do d=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
do c=nO+1,nOrb-nR
do d=nO+1,nOrb-nR
eps1 = eGF3(p) + e0(i) - e0(a) - e0(b)
eps2 = eGF3(p) + e0(i) - e0(c) - e0(d)
@ -191,12 +191,12 @@
end do
end do
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
eps1 = eGF3(p) + e0(i) - e0(a) - e0(b)
eps2 = e0(j) + e0(k) - e0(a) - e0(b)
@ -213,12 +213,12 @@
Cpp(:,3) = Cpp(:,2)
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do c=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
do c=nO+1,nOrb-nR
eps1 = eGF3(p) + e0(a) - e0(i) - e0(j)
eps2 = e0(i) + e0(j) - e0(b) - e0(c)
@ -234,12 +234,12 @@
Cpp(:,5) = Cpp(:,4)
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do l=nC+1,nO
do a=nO+1,nBas-nR
do a=nO+1,nOrb-nR
eps1 = eGF3(p) + e0(a) - e0(i) - e0(j)
eps2 = eGF3(p) + e0(a) - e0(k) - e0(l)
@ -257,12 +257,12 @@
Dpp(:,:) = 0d0
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do c=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
do c=nO+1,nOrb-nR
eps1 = eGF3(p) + e0(i) - e0(a) - e0(b)
eps2 = eGF3(p) + e0(j) - e0(b) - e0(c)
@ -282,12 +282,12 @@
end do
end do
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do c=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
do c=nO+1,nOrb-nR
eps1 = eGF3(p) + e0(i) - e0(a) - e0(c)
eps2 = e0(i) + e0(j) - e0(a) - e0(b)
@ -309,12 +309,12 @@
Dpp(:,3) = Dpp(:,2)
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
eps1 = eGF3(p) + e0(a) - e0(j) - e0(k)
eps2 = e0(i) + e0(j) - e0(a) - e0(b)
@ -336,12 +336,12 @@
Dpp(:,5) = Dpp(:,4)
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
eps1 = eGF3(p) + e0(a) - e0(i) - e0(k)
eps2 = eGF3(p) + e0(b) - e0(j) - e0(k)
@ -380,7 +380,7 @@
Z(:) = Cpp(:,2) + Cpp(:,3) + Cpp(:,4) + Cpp(:,5) &
+ Dpp(:,2) + Dpp(:,3) + Dpp(:,4) + Dpp(:,5)
Z(nC+1:nBas-nR) = Z(nC+1:nBas-nR)/Sig2(nC+1:nBas-nR)
Z(nC+1:nOrb-nR) = Z(nC+1:nOrb-nR)/Sig2(nC+1:nOrb-nR)
Z(:) = 1d0/(1d0 - Z(:))
Sig3(:) = Z(:)*Sig3(:)
@ -393,8 +393,8 @@
X2h1p(:) = Cpp(:,4) + Cpp(:,5) + Dpp(:,4) + Dpp(:,5)
X1h2p(:) = Cpp(:,2) + Cpp(:,3) + Dpp(:,2) + Dpp(:,3)
X2h1p(nC+1:nBas-nR) = X2h1p(nC+1:nBas-nR)/Bpp(nC+1:nBas-nR,1)
X1h2p(nC+1:nBas-nR) = X1h2p(nC+1:nBas-nR)/Bpp(nC+1:nBas-nR,2)
X2h1p(nC+1:nOrb-nR) = X2h1p(nC+1:nOrb-nR)/Bpp(nC+1:nOrb-nR,1)
X1h2p(nC+1:nOrb-nR) = X1h2p(nC+1:nOrb-nR)/Bpp(nC+1:nOrb-nR,2)
Sig3(:) = SigInf(:) + &
+ 1d0/(1d0 - X2h1p(:))*Sig2h1p(:) &
@ -412,11 +412,11 @@
X2h1p(:) = Cpp(:,4) + Cpp(:,5) + Dpp(:,4) + Dpp(:,5)
X1h2p(:) = Cpp(:,2) + Cpp(:,3) + Dpp(:,2) + Dpp(:,3)
X2h1p(nC+1:nBas-nR) = X2h1p(nC+1:nBas-nR)/Bpp(nC+1:nBas-nR,1)
X1h2p(nC+1:nBas-nR) = X1h2p(nC+1:nBas-nR)/Bpp(nC+1:nBas-nR,2)
X2h1p(nC+1:nOrb-nR) = X2h1p(nC+1:nOrb-nR)/Bpp(nC+1:nOrb-nR,1)
X1h2p(nC+1:nOrb-nR) = X1h2p(nC+1:nOrb-nR)/Bpp(nC+1:nOrb-nR,2)
Z(:) = X2h1p(:)*Sig2h1p(:) + X1h2p(:)*Sig1h2p(:)
Z(nC+1:nBas-nR) = Z(nC+1:nBas-nR)/(Sig3(nC+1:nBas-nR) - SigInf(nC+1:nBas-nR))
Z(nC+1:nOrb-nR) = Z(nC+1:nOrb-nR)/(Sig3(nC+1:nOrb-nR) - SigInf(nC+1:nOrb-nR))
Z(:) = 1d0/(1d0 - Z(:))
Sig3(:) = Z(:)*Sig3(:)
@ -429,6 +429,6 @@
! Print results
call print_G0F3(nBas,nO,e0,Z,eGF3)
call print_G0F3(nOrb,nO,e0,Z,eGF3)
end subroutine

44
src/GF/RGF.f90
View File

@ -1,10 +1,7 @@
! ---
subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, renorm, maxSCF, &
thresh, max_diis, dophBSE, doppBSE, TDA, dBSE, dTDA, singlet, triplet, linearize, &
eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, EHF, &
S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, epsHF)
subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,maxSCF, &
thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize, &
eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF, &
S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
! Green's function module
@ -42,15 +39,16 @@ subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, ren
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: EHF
double precision,intent(in) :: epsHF(nOrb)
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: cHF(nBas,nOrb)
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
@ -74,9 +72,9 @@ subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, ren
if(doG0F2) then
call wall_time(start_GF)
call RG0F2(dotest, dophBSE, doppBSE, TDA, dBSE, dTDA, singlet, triplet, &
linearize, eta, regularize, nOrb, nC, nO, nV, nR, nS, &
ENuc, EHF, ERI_MO, dipole_int_MO, epsHF)
call RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS, &
ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
@ -92,9 +90,9 @@ subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, ren
if(doevGF2) then
call wall_time(start_GF)
call evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, &
singlet,triplet,linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,EHF, &
ERI_MO,dipole_int_MO,epsHF)
call evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, &
singlet,triplet,linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS, &
ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
@ -110,10 +108,10 @@ subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, ren
if(doqsGF2) then
call wall_time(start_GF)
call qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, &
dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, &
rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, EHF, S, &
X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, epsHF)
call qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA, &
dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc, &
rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, &
ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
@ -129,7 +127,7 @@ subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, ren
if(doufG0F02) then
call wall_time(start_GF)
call ufRG0F02(dotest, nOrb, nC, nO, nV, nR, nS, ENuc, EHF, ERI_MO, epsHF)
call ufRG0F02(dotest,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
@ -145,7 +143,7 @@ subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, ren
if(doG0F3) then
call wall_time(start_GF)
call RG0F3(dotest, renorm, nOrb, nC, nO, nV, nR, ERI_MO, epsHF)
call RG0F3(dotest,renorm,nBas,nOrb,nC,nO,nV,nR,ERI_MO, eHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
@ -161,7 +159,7 @@ subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, ren
if(doevGF3) then
call wall_time(start_GF)
call evRGF3(dotest, maxSCF, thresh, max_diis, renorm, nOrb, nC, nO, nV, nR, ERI_MO, epsHF)
call evRGF3(dotest,maxSCF,thresh,max_diis,renorm,nBas,nOrb,nC,nO,nV,nR,ERI_MO,eHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF

27
src/GF/evRGF2.f90
View File

@ -1,5 +1,5 @@
subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform eigenvalue self-consistent second-order Green function calculation
@ -25,6 +25,7 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nO
integer,intent(in) :: nC
integer,intent(in) :: nV
@ -32,9 +33,9 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart)
! Local variables
@ -62,7 +63,7 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
! Memory allocation
allocate(SigC(nBas), Z(nBas), eGF(nBas), eOld(nBas), error_diis(nBas,max_diis), e_diis(nBas,max_diis))
allocate(SigC(nOrb), Z(nOrb), eGF(nOrb), eOld(nOrb), error_diis(nOrb,max_diis), e_diis(nOrb,max_diis))
! Initialization
@ -85,11 +86,11 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
if(regularize) then
call RGF2_reg_self_energy_diag(eta,nBas,nC,nO,nV,nR,eGF,ERI,SigC,Z)
call RGF2_reg_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eGF,ERI,SigC,Z)
else
call RGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,eGF,ERI,SigC,Z)
call RGF2_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eGF,ERI,SigC,Z)
end if
@ -102,7 +103,7 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
write(*,*) ' *** Quasiparticle energies obtained by root search *** '
write(*,*)
call RGF2_QP_graph(eta,nBas,nC,nO,nV,nR,eHF,ERI,eOld,eOld,eGF,Z)
call RGF2_QP_graph(eta,nOrb,nC,nO,nV,nR,eHF,ERI,eOld,eOld,eGF,Z)
end if
@ -110,13 +111,13 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
! Print results
call RMP2(.false.,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
call print_evRGF2(nBas,nO,nSCF,Conv,eHF,SigC,Z,eGF,ENuc,ERHF,Ec)
call RMP2(.false.,regularize,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
call print_evRGF2(nOrb,nO,nSCF,Conv,eHF,SigC,Z,eGF,ENuc,ERHF,Ec)
! DIIS extrapolation
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGF-eOld,eGF)
call DIIS_extrapolation(rcond,nOrb,nOrb,n_diis,error_diis,e_diis,eGF-eOld,eGF)
if(abs(rcond) < 1d-15) n_diis = 0
@ -149,7 +150,7 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
if(dophBSE) then
call RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
call RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
@ -166,7 +167,7 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
if(doppBSE) then
call RGF2_ppBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
call RGF2_ppBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

126
src/GF/evRGF3.f90
View File

@ -1,4 +1,4 @@
subroutine evRGF3(dotest,maxSCF,thresh,max_diis,renormalization,nBas,nC,nO,nV,nR,V,e0)
subroutine evRGF3(dotest,maxSCF,thresh,max_diis,renormalization,nBas,nOrb,nC,nO,nV,nR,V,e0)
! Perform third-order Green function calculation in diagonal approximation
@ -11,8 +11,8 @@
double precision,intent(in) :: thresh
integer,intent(in) :: maxSCF,max_diis,renormalization
integer,intent(in) :: nBas,nC,nO,nV,nR
double precision,intent(in) :: e0(nBas),V(nBas,nBas,nBas,nBas)
integer,intent(in) :: nBas,nOrb,nC,nO,nV,nR
double precision,intent(in) :: e0(nOrb),V(nOrb,nOrb,nOrb,nOrb)
! Local variables
@ -37,11 +37,11 @@
! Memory allocation
allocate(eGF3(nBas),eOld(nBas), &
Sig2(nBas),SigInf(nBas),Sig3(nBas), &
App(nBas,6),Bpp(nBas,2),Cpp(nBas,6),Dpp(nBas,6), &
Z(nBas),X2h1p(nBas),X1h2p(nBas),Sig2h1p(nBas),Sig1h2p(nBas), &
error_diis(nBas,max_diis),e_diis(nBas,max_diis))
allocate(eGF3(nOrb),eOld(nOrb), &
Sig2(nOrb),SigInf(nOrb),Sig3(nOrb), &
App(nOrb,6),Bpp(nOrb,2),Cpp(nOrb,6),Dpp(nOrb,6), &
Z(nOrb),X2h1p(nOrb),X1h2p(nOrb),Sig2h1p(nOrb),Sig1h2p(nOrb), &
error_diis(nOrb,max_diis),e_diis(nOrb,max_diis))
!------------------------------------------------------------------------
! Compute third-order frequency-independent contribution
@ -49,12 +49,12 @@
App(:,:) = 0d0
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
eps1 = e0(j) + e0(i) - e0(a) - e0(b)
eps2 = e0(k) + e0(i) - e0(a) - e0(b)
@ -69,12 +69,12 @@
end do
end do
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do c=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
do c=nO+1,nOrb-nR
eps1 = e0(j) + e0(i) - e0(a) - e0(b)
eps2 = e0(j) + e0(i) - e0(a) - e0(c)
@ -89,12 +89,12 @@
end do
end do
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do c=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
do c=nO+1,nOrb-nR
eps1 = e0(j) + e0(i) - e0(a) - e0(b)
eps2 = e0(j) - e0(c)
@ -111,12 +111,12 @@
App(:,4) = App(:,3)
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
eps1 = e0(j) + e0(i) - e0(a) - e0(b)
eps2 = e0(k) - e0(b)
@ -157,10 +157,10 @@
Bpp(:,:) = 0d0
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do a=nO+1,nOrb-nR
eps = eGF3(p) + e0(a) - e0(i) - e0(j)
@ -172,10 +172,10 @@
end do
end do
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
eps = eGF3(p) + e0(i) - e0(a) - e0(b)
@ -195,12 +195,12 @@
Cpp(:,:) = 0d0
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do c=nO+1,nBas-nR
do d=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
do c=nO+1,nOrb-nR
do d=nO+1,nOrb-nR
eps1 = eGF3(p) + e0(i) - e0(a) - e0(b)
eps2 = eGF3(p) + e0(i) - e0(c) - e0(d)
@ -215,12 +215,12 @@
end do
end do
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
eps1 = eGF3(p) + e0(i) - e0(a) - e0(b)
eps2 = e0(j) + e0(k) - e0(a) - e0(b)
@ -237,12 +237,12 @@
Cpp(:,3) = Cpp(:,2)
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do c=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
do c=nO+1,nOrb-nR
eps1 = eGF3(p) + e0(a) - e0(i) - e0(j)
eps2 = e0(i) + e0(j) - e0(b) - e0(c)
@ -258,12 +258,12 @@
Cpp(:,5) = Cpp(:,4)
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do l=nC+1,nO
do a=nO+1,nBas-nR
do a=nO+1,nOrb-nR
eps1 = eGF3(p) + e0(a) - e0(i) - e0(j)
eps2 = eGF3(p) + e0(a) - e0(k) - e0(l)
@ -281,12 +281,12 @@
Dpp(:,:) = 0d0
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do c=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
do c=nO+1,nOrb-nR
eps1 = eGF3(p) + e0(i) - e0(a) - e0(b)
eps2 = eGF3(p) + e0(j) - e0(b) - e0(c)
@ -306,12 +306,12 @@
end do
end do
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do c=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
do c=nO+1,nOrb-nR
eps1 = eGF3(p) + e0(i) - e0(a) - e0(c)
eps2 = e0(i) + e0(j) - e0(a) - e0(b)
@ -333,12 +333,12 @@
Dpp(:,3) = Dpp(:,2)
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
eps1 = eGF3(p) + e0(a) - e0(j) - e0(k)
eps2 = e0(i) + e0(j) - e0(a) - e0(b)
@ -360,12 +360,12 @@
Dpp(:,5) = Dpp(:,4)
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
eps1 = eGF3(p) + e0(a) - e0(i) - e0(k)
eps2 = eGF3(p) + e0(b) - e0(j) - e0(k)
@ -404,7 +404,7 @@
Z(:) = Cpp(:,2) + Cpp(:,3) + Cpp(:,4) + Cpp(:,5) &
+ Dpp(:,2) + Dpp(:,3) + Dpp(:,4) + Dpp(:,5)
Z(nC+1:nBas-nR) = Z(nC+1:nBas-nR)/Sig2(nC+1:nBas-nR)
Z(nC+1:nOrb-nR) = Z(nC+1:nOrb-nR)/Sig2(nC+1:nOrb-nR)
Z(:) = 1d0/(1d0 - Z(:))
Sig3(:) = Z(:)*Sig3(:)
@ -417,8 +417,8 @@
X2h1p(:) = Cpp(:,4) + Cpp(:,5) + Dpp(:,4) + Dpp(:,5)
X1h2p(:) = Cpp(:,2) + Cpp(:,3) + Dpp(:,2) + Dpp(:,3)
X2h1p(nC+1:nBas-nR) = X2h1p(nC+1:nBas-nR)/Bpp(nC+1:nBas-nR,1)
X1h2p(nC+1:nBas-nR) = X1h2p(nC+1:nBas-nR)/Bpp(nC+1:nBas-nR,2)
X2h1p(nC+1:nOrb-nR) = X2h1p(nC+1:nOrb-nR)/Bpp(nC+1:nOrb-nR,1)
X1h2p(nC+1:nOrb-nR) = X1h2p(nC+1:nOrb-nR)/Bpp(nC+1:nOrb-nR,2)
Sig3(:) = SigInf(:) + &
+ 1d0/(1d0 - X2h1p(:))*Sig2h1p(:) &
@ -436,11 +436,11 @@
X2h1p(:) = Cpp(:,4) + Cpp(:,5) + Dpp(:,4) + Dpp(:,5)
X1h2p(:) = Cpp(:,2) + Cpp(:,3) + Dpp(:,2) + Dpp(:,3)
X2h1p(nC+1:nBas-nR) = X2h1p(nC+1:nBas-nR)/Bpp(nC+1:nBas-nR,1)
X1h2p(nC+1:nBas-nR) = X1h2p(nC+1:nBas-nR)/Bpp(nC+1:nBas-nR,2)
X2h1p(nC+1:nOrb-nR) = X2h1p(nC+1:nOrb-nR)/Bpp(nC+1:nOrb-nR,1)
X1h2p(nC+1:nOrb-nR) = X1h2p(nC+1:nOrb-nR)/Bpp(nC+1:nOrb-nR,2)
Z(:) = X2h1p(:)*Sig2h1p(:) + X1h2p(:)*Sig1h2p(:)
Z(nC+1:nBas-nR) = Z(nC+1:nBas-nR)/(Sig3(nC+1:nBas-nR) - SigInf(nC+1:nBas-nR))
Z(nC+1:nOrb-nR) = Z(nC+1:nOrb-nR)/(Sig3(nC+1:nOrb-nR) - SigInf(nC+1:nOrb-nR))
Z(:) = 1d0/(1d0 - Z(:))
Sig3(:) = Z(:)*Sig3(:)
@ -457,12 +457,12 @@
! Print results
call print_evGF3(nBas,nO,nSCF,Conv,e0,Z,eGF3)
call print_evGF3(nOrb,nO,nSCF,Conv,e0,Z,eGF3)
! DIIS extrapolation
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGF3-eOld,eGF3)
call DIIS_extrapolation(rcond,nOrb,nOrb,n_diis,error_diis,e_diis,eGF3-eOld,eGF3)
if(abs(rcond) < 1d-15) n_diis = 0

11
src/GF/qsRGF2.f90
View File

@ -1,10 +1,7 @@
! ---
subroutine qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, &
dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, &
rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, &
S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA, &
dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc, &
rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, &
ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
! Perform a quasiparticle self-consistent GF2 calculation

31
src/GF/ufRG0F02.f90
View File

@ -1,4 +1,4 @@
subroutine ufRG0F02(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
subroutine ufRG0F02(dotest,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Unfold G0F02
@ -10,6 +10,7 @@ subroutine ufRG0F02(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
logical,intent(in) :: dotest
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
@ -17,8 +18,8 @@ subroutine ufRG0F02(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eHF(nOrb)
! Local variables
@ -103,7 +104,7 @@ subroutine ufRG0F02(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
ija = 0
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do a=nO+1,nOrb-nR
ija = ija + 1
H(1 ,1+ija) = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))
@ -119,8 +120,8 @@ subroutine ufRG0F02(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
iab = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
iab = iab + 1
H(1 ,1+n2h1p+iab) = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))
@ -137,13 +138,13 @@ subroutine ufRG0F02(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
ija = 0
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do a=nO+1,nOrb-nR
ija = ija + 1
klc = 0
do k=nC+1,nO
do l=nC+1,nO
do c=nO+1,nBas-nR
do c=nO+1,nOrb-nR
klc = klc + 1
H(1+ija,1+klc) &
@ -163,14 +164,14 @@ subroutine ufRG0F02(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
iab = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
iab = iab + 1
kcd = 0
do k=nC+1,nO
do c=nO+1,nBas-nR
do d=nO+1,nBas-nR
do c=nO+1,nOrb-nR
do d=nO+1,nOrb-nR
kcd = kcd + 1
H(1+n2h1p+iab,1+n2h1p+kcd) &
@ -260,7 +261,7 @@ subroutine ufRG0F02(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
ija = 0
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do a=nO+1,nOrb-nR
ija = ija + 1
if(abs(Reigv(1+ija,s)) > cutoff2) &
@ -274,8 +275,8 @@ subroutine ufRG0F02(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
iab = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
iab = iab + 1
if(abs(Reigv(1+n2h1p+iab,s)) > cutoff2) &

33
src/GT/RG0T0eh.f90
View File

@ -1,5 +1,5 @@
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,ERI,dipole_int,eHF)
singlet,triplet,linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform ehG0T0 calculation
@ -28,6 +28,7 @@ subroutine RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,T
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
@ -35,9 +36,9 @@ subroutine RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,T
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart)
double precision,intent(in) :: eHF(nOrb)
! Local variables
@ -99,8 +100,8 @@ subroutine RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,T
! Memory allocation
allocate(Aph(nS,nS),Bph(nS,nS),Sig(nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), &
rhoL(nBas,nBas,nS),rhoR(nBas,nBas,nS),eGT(nBas),eGTlin(nBas))
allocate(Aph(nS,nS),Bph(nS,nS),Sig(nOrb),Z(nOrb),Om(nS),XpY(nS,nS),XmY(nS,nS), &
rhoL(nOrb,nOrb,nS),rhoR(nOrb,nOrb,nS),eGT(nOrb),eGTlin(nOrb))
!---------------------------------
! Compute (triplet) RPA screening
@ -108,8 +109,8 @@ subroutine RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,T
ispin = 2
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
if(.not.TDA_T) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
if(.not.TDA_T) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
@ -119,15 +120,15 @@ subroutine RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,T
! Compute spectral weights !
!--------------------------!
call RGTeh_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,XmY,rhoL,rhoR)
call RGTeh_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,XmY,rhoL,rhoR)
!------------------------!
! Compute GW self-energy !
!------------------------!
if(regularize) call GTeh_regularization(nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR)
if(regularize) call GTeh_regularization(nOrb,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR)
call RGTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,EcGM,Sig,Z)
call RGTeh_self_energy_diag(eta,nOrb,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,EcGM,Sig,Z)
!-----------------------------------!
! Solve the quasi-particle equation !
@ -149,16 +150,16 @@ subroutine RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,T
write(*,*) ' *** Quasiparticle energies obtained by root search *** '
write(*,*)
call RGTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eHF,eGT,Z)
call RGTeh_QP_graph(eta,nOrb,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eGTlin,eHF,eGT,Z)
end if
! call RGTeh_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eHF,Om,rhoL,rhoR)
! call RGTeh_plot_self_energy(nOrb,nC,nO,nV,nR,nS,eHF,eHF,Om,rhoL,rhoR)
! Compute the RPA correlation energy based on the G0T0eh quasiparticle energies
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
if(.not.TDA_T) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
if(.not.TDA_T) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
@ -166,7 +167,7 @@ subroutine RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,T
! Dump results !
!--------------!
call print_RG0T0eh(nBas,nO,eHF,ENuc,ERHF,Sig,Z,eGT,EcRPA,EcGM)
call print_RG0T0eh(nOrb,nO,eHF,ENuc,ERHF,Sig,Z,eGT,EcRPA,EcGM)
! Testing zone

56
src/GT/RG0T0pp.f90
View File

@ -1,5 +1,5 @@
subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform one-shot calculation with a T-matrix self-energy (G0T0)
@ -25,6 +25,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
@ -40,8 +41,9 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
! Local variables
logical :: print_T = .false.
double precision :: lambda
integer :: ispin
integer :: isp_T
integer :: iblock
integer :: nOOs,nOOt
integer :: nVVs,nVVt
@ -79,6 +81,9 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
write(*,*)'*********************************'
write(*,*)
! Initialization
lambda = 1d0
! TDA for T
@ -87,13 +92,6 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
write(*,*)
end if
! TDA
if(TDA) then
write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*)
end if
! Dimensions of the pp-RPA linear reponse matrices
!nOOs = nO*(nO + 1)/2
@ -119,7 +117,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
! alpha-beta block
!----------------------------------------------
ispin = 1
isp_T = 1
!iblock = 1
iblock = 3
@ -127,11 +125,11 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp)
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,lambda,eHF,ERI,Cpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,lambda,eHF,ERI,Dpp)
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,lambda,ERI,Bpp)
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T))
deallocate(Bpp,Cpp,Dpp)
!print*, 'LAPACK:'
!print*, Om2s
@ -162,7 +160,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
! alpha-alpha block
!----------------------------------------------
ispin = 2
isp_T = 2
! iblock = 2
iblock = 4
@ -170,11 +168,11 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp)
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,lambda,eHF,ERI,Cpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,lambda,eHF,ERI,Dpp)
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,lambda,ERI,Bpp)
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
deallocate(Bpp,Cpp,Dpp)
!print*, 'LAPACK:'
!print*, Om2t
@ -259,31 +257,31 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
! Compute the ppRPA correlation energy
ispin = 1
isp_T = 1
! iblock = 1
iblock = 3
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,lambda,eGT,ERI,Cpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,lambda,eGT,ERI,Dpp)
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,lambda,ERI,Bpp)
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T))
deallocate(Bpp,Cpp,Dpp)
ispin = 2
isp_T = 2
! iblock = 2
iblock = 4
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp)
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,lambda,eGT,ERI,Cpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,lambda,eGT,ERI,Dpp)
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,lambda,ERI,Bpp)
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
deallocate(Bpp,Cpp,Dpp)

22
src/GT/RGT.f90
View File

@ -82,7 +82,7 @@ subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevG
call wall_time(start_GT)
call RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
@ -99,7 +99,7 @@ subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevG
call wall_time(start_GT)
call evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
@ -115,9 +115,9 @@ subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevG
if(doqsGTpp) then
call wall_time(start_GT)
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, nOrb, nC, nO, &
nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
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,nOrb,nC,nO, &
nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
@ -150,7 +150,7 @@ subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevG
call wall_time(start_GT)
call RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
@ -167,7 +167,7 @@ subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevG
call wall_time(start_GT)
call evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
singlet,triplet,linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
@ -183,10 +183,10 @@ subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevG
if(doqsGTeh) then
call wall_time(start_GT)
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, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, V, &
Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
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,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T,V, &
Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT

9
src/GT/RGTpp_phBSE.f90
View File

@ -81,6 +81,15 @@ subroutine RGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,
allocate(Aph(nS,nS),Bph(nS,nS),TAab(nS,nS),TBab(nS,nS),TAaa(nS,nS),TBaa(nS,nS), &
OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS))
!-----!
! TDA !
!-----!
if(TDA) then
write(*,*) 'Tamm-Dancoff approximation activated in phBSE!'
write(*,*)
end if
!---------------------------------------!
! Compute T-matrix for alpha-beta block !
!---------------------------------------!

47
src/GT/evRGTeh.f90
View File

@ -1,5 +1,5 @@
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,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform self-consistent eigenvalue-only ehGT calculation
@ -32,14 +32,15 @@ subroutine evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart)
! Local variables
@ -98,8 +99,8 @@ subroutine evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Memory allocation
allocate(Aph(nS,nS),Bph(nS,nS),eGT(nBas),eOld(nBas),Z(nBas),Sig(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), &
rhoL(nBas,nBas,nS),rhoR(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis))
allocate(Aph(nS,nS),Bph(nS,nS),eGT(nOrb),eOld(nOrb),Z(nOrb),Sig(nOrb),Om(nS),XpY(nS,nS),XmY(nS,nS), &
rhoL(nOrb,nOrb,nS),rhoR(nOrb,nOrb,nS),error_diis(nOrb,max_diis),e_diis(nOrb,max_diis))
! Initialization
@ -122,8 +123,8 @@ subroutine evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Compute screening
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
if(.not.TDA_T) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI,Aph)
if(.not.TDA_T) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
@ -131,13 +132,13 @@ subroutine evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Compute spectral weights
call RGTeh_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,XmY,rhoL,rhoR)
call RGTeh_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,XmY,rhoL,rhoR)
! Compute correlation part of the self-energy
if(regularize) call GTeh_regularization(nBas,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR)
if(regularize) call GTeh_regularization(nOrb,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR)
call RGTeh_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR,EcGM,Sig,Z)
call RGTeh_self_energy_diag(eta,nOrb,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR,EcGM,Sig,Z)
! Solve the quasi-particle equation
@ -155,7 +156,7 @@ subroutine evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
write(*,*) ' *** Quasiparticle energies obtained by root search *** '
write(*,*)
call RGTeh_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eOld,eOld,eGT,Z)
call RGTeh_QP_graph(eta,nOrb,nC,nO,nV,nR,nS,eHF,Om,rhoL,rhoR,eOld,eOld,eGT,Z)
end if
@ -165,7 +166,7 @@ subroutine evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Print results
call print_evRGTeh(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,Sig,Z,eGT,EcRPA,EcGM)
call print_evRGTeh(nOrb,nO,nSCF,Conv,eHF,ENuc,ERHF,Sig,Z,eGT,EcRPA,EcGM)
! Linear mixing or DIIS extrapolation
@ -177,7 +178,7 @@ subroutine evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
n_diis = min(n_diis+1,max_diis)
if(abs(rcond) > 1d-7) then
call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGT-eOld,eGT)
call DIIS_extrapolation(rcond,nOrb,nOrb,n_diis,error_diis,e_diis,eGT-eOld,eGT)
else
n_diis = 0
end if
@ -219,7 +220,7 @@ subroutine evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! if(BSE) then
! call Bethe_Salpeter(BSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGW,eGW,EcBSE)
! call Bethe_Salpeter(BSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eGW,eGW,EcBSE)
! if(exchange_kernel) then
@ -253,7 +254,7 @@ subroutine evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! end if
! call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcBSE)
! call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcBSE)
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
@ -270,7 +271,7 @@ subroutine evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! if(ppBSE) then
! call Bethe_Salpeter_pp(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcppBSE)
! call Bethe_Salpeter_pp(TDA_W,TDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcppBSE)
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
@ -281,20 +282,20 @@ subroutine evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
! nBas2 = 2*nBas
! nOrb2 = 2*nOrb
! nO2 = 2*nO
! nV2 = 2*nV
! nC2 = 2*nC
! nR2 = 2*nR
! nS2 = nO2*nV2
!
! allocate(seHF(nBas2),seGW(nBas2),sERI(nBas2,nBas2,nBas2,nBas2))
! allocate(seHF(nOrb2),seGW(nOrb2),sERI(nOrb2,nOrb2,nOrb2,nOrb2))
!
! call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF)
! call spatial_to_spin_MO_energy(nBas,eGW,nBas2,seGW)
! call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI)
! call spatial_to_spin_MO_energy(nOrb,eHF,nOrb2,seHF)
! call spatial_to_spin_MO_energy(nOrb,eGW,nOrb2,seGW)
! call spatial_to_spin_ERI(nOrb,ERI,nOrb2,sERI)
!
! call Bethe_Salpeter_pp_so(TDA_W,TDA,singlet,triplet,eta,nBas2,nC2,nO2,nV2,nR2,nS2,sERI,dipole_int,seHF,seGW,EcppBSE)
! call Bethe_Salpeter_pp_so(TDA_W,TDA,singlet,triplet,eta,nOrb2,nC2,nO2,nV2,nR2,nS2,sERI,dipole_int,seHF,seGW,EcppBSE)
! end if

56
src/GT/evRGTpp.f90
View File

@ -1,5 +1,5 @@
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,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform eigenvalue self-consistent calculation with a T-matrix self-energy (evGT)
@ -28,6 +28,7 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
@ -35,9 +36,9 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart)
! Local variables
@ -88,13 +89,6 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
write(*,*)
end if
! TDA
if(TDA) then
write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*)
end if
! Dimensions of the pp-RPA linear reponse matrices
nOOs = nO*nO
@ -107,12 +101,12 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), &
Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), &
rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), &
rho1s(nOrb,nOrb,nVVs),rho2s(nOrb,nOrb,nOOs), &
Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), &
Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), &
rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt), &
eGT(nBas),eOld(nBas),Z(nBas),Sig(nBas), &
error_diis(nBas,max_diis),e_diis(nBas,max_diis))
rho1t(nOrb,nOrb,nVVt),rho2t(nOrb,nOrb,nOOt), &
eGT(nOrb),eOld(nOrb),Z(nOrb),Sig(nOrb), &
error_diis(nOrb,max_diis),e_diis(nOrb,max_diis))
! Initialization
@ -143,9 +137,9 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
@ -162,9 +156,9 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp)
call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp)
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp)
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
@ -178,21 +172,21 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
!----------------------------------------------
iblock = 3
call RGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
call RGTpp_excitation_density(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
iblock = 4
call RGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
call RGTpp_excitation_density(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
!----------------------------------------------
! Compute T-matrix version of the self-energy
!----------------------------------------------
if(regularize) then
call GTpp_regularization(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s)
call GTpp_regularization(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t)
call GTpp_regularization(eta,nOrb,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s)
call GTpp_regularization(eta,nOrb,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t)
end if
call RGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nooS,nVVt,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s, &
call RGTpp_self_energy_diag(eta,nOrb,nC,nO,nV,nR,nooS,nVVt,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s, &
Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z)
!----------------------------------------------
@ -211,7 +205,7 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
write(*,*) ' *** Quasiparticle energies obtained by root search *** '
write(*,*)
call RGTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, &
call RGTpp_QP_graph(eta,nOrb,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, &
Om1t,rho1t,Om2t,rho2t,eOld,eOld,eGT,Z)
end if
@ -224,12 +218,12 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
! Dump results
!----------------------------------------------
call print_evRGTpp(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,Sig,Z,eGT,EcGM,EcRPA)
call print_evRGTpp(nOrb,nO,nSCF,Conv,eHF,ENuc,ERHF,Sig,Z,eGT,EcGM,EcRPA)
! DIIS extrapolation
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGT(:)-eOld(:),eGT(:))
call DIIS_extrapolation(rcond,nOrb,nOrb,n_diis,error_diis,e_diis,eGT(:)-eOld(:),eGT(:))
! Reset DIIS if required
@ -252,7 +246,7 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
if(BSE) then
call RGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
call RGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, &
ERI,dipole_int,eGT,eGT,EcBSE)
@ -288,7 +282,7 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
end if
call RGTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, &
call RGTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS, &
nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, &
Om2t,X2t,Y2t,rho1t,rho2t,ERI,eGT,eGT,EcBSE)

7
src/GT/qsRGTpp.f90
View File

@ -132,13 +132,6 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d
write(*,*)
end if
! TDA
if(TDA) then
write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*)
end if
! Memory allocation
allocate(eGT(nOrb))