diff --git a/src/GF/RG0F2.f90 b/src/GF/RG0F2.f90 index 546a7d0..8ed0a38 100644 --- a/src/GF/RG0F2.f90 +++ b/src/GF/RG0F2.f90 @@ -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) diff --git a/src/GF/RG0F3.f90 b/src/GF/RG0F3.f90 index 6bc6c7b..8e587b2 100644 --- a/src/GF/RG0F3.f90 +++ b/src/GF/RG0F3.f90 @@ -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 diff --git a/src/GF/RGF.f90 b/src/GF/RGF.f90 index 0f1deaf..3493ff4 100644 --- a/src/GF/RGF.f90 +++ b/src/GF/RGF.f90 @@ -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 diff --git a/src/GF/evRGF2.f90 b/src/GF/evRGF2.f90 index 4af5c8b..ff5bc50 100644 --- a/src/GF/evRGF2.f90 +++ b/src/GF/evRGF2.f90 @@ -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(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GF/evRGF3.f90 b/src/GF/evRGF3.f90 index c056076..d9a02e6 100644 --- a/src/GF/evRGF3.f90 +++ b/src/GF/evRGF3.f90 @@ -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 diff --git a/src/GF/qsRGF2.f90 b/src/GF/qsRGF2.f90 index 1a0afc0..b2b10b9 100644 --- a/src/GF/qsRGF2.f90 +++ b/src/GF/qsRGF2.f90 @@ -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 diff --git a/src/GF/ufRG0F02.f90 b/src/GF/ufRG0F02.f90 index 5719022..2683e9c 100644 --- a/src/GF/ufRG0F02.f90 +++ b/src/GF/ufRG0F02.f90 @@ -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) & diff --git a/src/GT/RG0T0eh.f90 b/src/GT/RG0T0eh.f90 index 694017f..4d7710c 100644 --- a/src/GT/RG0T0eh.f90 +++ b/src/GT/RG0T0eh.f90 @@ -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 diff --git a/src/GT/RG0T0pp.f90 b/src/GT/RG0T0pp.f90 index cace39e..67dca13 100644 --- a/src/GT/RG0T0pp.f90 +++ b/src/GT/RG0T0pp.f90 @@ -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) diff --git a/src/GT/RGT.f90 b/src/GT/RGT.f90 index e194cac..c18d794 100644 --- a/src/GT/RGT.f90 +++ b/src/GT/RGT.f90 @@ -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 diff --git a/src/GT/RGTpp_phBSE.f90 b/src/GT/RGTpp_phBSE.f90 index 778cd8c..dfe1d76 100644 --- a/src/GT/RGTpp_phBSE.f90 +++ b/src/GT/RGTpp_phBSE.f90 @@ -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 ! !---------------------------------------! diff --git a/src/GT/evRGTeh.f90 b/src/GT/evRGTeh.f90 index 6a71e98..107e915 100644 --- a/src/GT/evRGTeh.f90 +++ b/src/GT/evRGTeh.f90 @@ -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 diff --git a/src/GT/evRGTpp.f90 b/src/GT/evRGTpp.f90 index cd2ea89..87b6421 100644 --- a/src/GT/evRGTpp.f90 +++ b/src/GT/evRGTpp.f90 @@ -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) diff --git a/src/GT/qsRGTpp.f90 b/src/GT/qsRGTpp.f90 index 89841bf..648a825 100644 --- a/src/GT/qsRGTpp.f90 +++ b/src/GT/qsRGTpp.f90 @@ -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))