diff --git a/input/methods b/input/methods index 8848a96..91ecc0d 100644 --- a/input/methods +++ b/input/methods @@ -9,7 +9,7 @@ # CIS* CIS(D) CID CISD FCI F F F F F # phRPA* phRPAx* crRPA ppRPA - F F F F + F F F T # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW diff --git a/input/options b/input/options index ce371d3..dc7e3fb 100644 --- a/input/options +++ b/input/options @@ -5,7 +5,7 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip - F T F T T + F T T T T # GF: maxSCF thresh DIIS n_diis lin eta renorm reg 256 0.00001 T 5 T 0.0 0 F # GW: maxSCF thresh DIIS n_diis lin eta COHSEX TDA_W reg @@ -13,6 +13,6 @@ # GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg 256 0.00001 T 5 T 0.0 F F # ACFDT: AC Kx XBS - F T T + T F F # BSE: phBSE phBSE2 ppBSE dBSE dTDA evDyn F F F T F F diff --git a/src/GT/G0T0pp.f90 b/src/GT/G0T0pp.f90 index c5c149d..5c4eac9 100644 --- a/src/GT/G0T0pp.f90 +++ b/src/GT/G0T0pp.f90 @@ -108,8 +108,8 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp allocate(Bpp(nVVab,nOOab),Cpp(nVVab,nVVab),Dpp(nOOab,nOOab)) if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,ERI_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,eHF,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,eHF,ERI_MO,Dpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVab,1d0,eHF,ERI_MO,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOab,1d0,eHF,ERI_MO,Dpp) call ppLR(TDA_T,nOOab,nVVab,Bpp,Cpp,Dpp,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin)) @@ -130,8 +130,8 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp allocate(Bpp(nVVaa,nOOaa),Cpp(nVVaa,nVVaa),Dpp(nOOaa,nOOaa)) if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,ERI_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,eHF,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,eHF,ERI_MO,Dpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVaa,1d0,eHF,ERI_MO,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOaa,1d0,eHF,ERI_MO,Dpp) call ppLR(TDA_T,nOOaa,nVVaa,Bpp,Cpp,Dpp,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin)) @@ -220,8 +220,8 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp allocate(Bpp(nVVab,nOOab),Cpp(nVVab,nVVab),Dpp(nOOab,nOOab)) if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,ERI_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,eGT,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOab,nVVab,1d0,eGT,ERI_MO,Dpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVab,1d0,eGT,ERI_MO,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOab,1d0,eGT,ERI_MO,Dpp) call ppLR(TDA_T,nOOab,nVVab,Bpp,Cpp,Dpp,Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,EcRPA(ispin)) @@ -233,8 +233,8 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp allocate(Bpp(nVVaa,nOOaa),Cpp(nVVaa,nVVaa),Dpp(nOOaa,nOOaa)) if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,ERI_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,eGT,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOaa,nVVaa,1d0,eGT,ERI_MO,Dpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVaa,1d0,eGT,ERI_MO,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOaa,1d0,eGT,ERI_MO,Dpp) call ppLR(TDA_T,nOOaa,nVVaa,Bpp,Cpp,Dpp,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,EcRPA(ispin)) diff --git a/src/GT/evGTpp.f90 b/src/GT/evGTpp.f90 index f821c3e..3cf172c 100644 --- a/src/GT/evGTpp.f90 +++ b/src/GT/evGTpp.f90 @@ -138,8 +138,8 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T 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_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eHF,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eHF,ERI_MO,Dpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eHF,ERI_MO,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eHF,ERI_MO,Dpp) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) @@ -157,8 +157,8 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T 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_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eHF,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eHF,ERI_MO,Dpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eHF,ERI_MO,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eHF,ERI_MO,Dpp) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) diff --git a/src/GT/qsGTpp.f90 b/src/GT/qsGTpp.f90 index 0f8da08..33ddb83 100644 --- a/src/GT/qsGTpp.f90 +++ b/src/GT/qsGTpp.f90 @@ -194,8 +194,8 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T 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_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,ERI_MO,Dpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) @@ -207,8 +207,8 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T 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_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO,Dpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) diff --git a/src/LR/ppLR.f90 b/src/LR/ppLR.f90 index e63dbda..b921c74 100644 --- a/src/LR/ppLR.f90 +++ b/src/LR/ppLR.f90 @@ -1,6 +1,6 @@ subroutine ppLR(TDA,nOO,nVV,B,C,D,Om1,X1,Y1,Om2,X2,Y2,EcRPA) -! Compute the p-p channel of the linear response: see Scuseria et al. JCP 139, 104113 (2013) +! Compute the pp channel of the linear response: see Scuseria et al. JCP 139, 104113 (2013) implicit none include 'parameters.h' @@ -16,8 +16,6 @@ subroutine ppLR(TDA,nOO,nVV,B,C,D,Om1,X1,Y1,Om2,X2,Y2,EcRPA) ! Local variables - integer :: ab,cd,ij,kl - integer :: p,q,r,s double precision :: trace_matrix double precision :: EcRPA1 double precision :: EcRPA2 @@ -77,15 +75,15 @@ subroutine ppLR(TDA,nOO,nVV,B,C,D,Om1,X1,Y1,Om2,X2,Y2,EcRPA) ! Split the various quantities in p-p and h-h parts - call sort_ppRPA(nOO,nVV,Om(:),Z(:,:),Om1(:),X1(:,:),Y1(:,:),Om2(:),X2(:,:),Y2(:,:)) + call sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) end if ! Compute the RPA correlation energy - EcRPA = 0.5d0*( sum(Om1(:)) - sum(Om2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) ) - EcRPA1 = +sum(Om1(:)) - trace_matrix(nVV,C(:,:)) - EcRPA2 = -sum(Om2(:)) - trace_matrix(nOO,D(:,:)) + EcRPA = 0.5d0*( sum(Om1) - sum(Om2) - trace_matrix(nVV,C) - trace_matrix(nOO,D) ) + EcRPA1 = +sum(Om1) - trace_matrix(nVV,C) + EcRPA2 = -sum(Om2) - trace_matrix(nOO,D) if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) & print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!' diff --git a/src/RPA/ppACFDT.f90 b/src/RPA/ppACFDT.f90 index adaac80..ddb606e 100644 --- a/src/RPA/ppACFDT.f90 +++ b/src/RPA/ppACFDT.f90 @@ -1,4 +1,4 @@ -subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) +subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ERI,e,EcAC) ! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem for pp sector @@ -17,7 +17,6 @@ subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR - integer,intent(in) :: nS double precision,intent(in) :: e(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) @@ -32,6 +31,9 @@ subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) integer :: nOO integer :: nVV + double precision,allocatable :: Bpp(:,:) + double precision,allocatable :: Cpp(:,:) + double precision,allocatable :: Dpp(:,:) double precision,allocatable :: Om1(:) double precision,allocatable :: X1(:,:) double precision,allocatable :: Y1(:,:) @@ -64,7 +66,8 @@ subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) nVV = nV*(nV+1)/2 allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),rho1(nBas,nBas,nVV), & - Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO),rho2(nBas,nBas,nOO)) + Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO),rho2(nBas,nBas,nOO), & + Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)) write(*,*) '--------------' write(*,*) 'Singlet states' @@ -79,9 +82,13 @@ subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) lambda = rAC(iAC) - call ppLR(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Om1,X1,Y1,Om2,X2,Y2,EcAC(ispin)) + if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp) + call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp) + call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp) - call ppACFDT_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOO,nVV,X1,Y1,X2,Y2,Ec(iAC,ispin)) + call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcAC(ispin)) + + call ppACFDT_correlation_energy(ispin,nBas,nC,nO,nV,nR,ERI,nOO,nVV,X1,Y1,X2,Y2,Ec(iAC,ispin)) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) @@ -94,7 +101,7 @@ subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) write(*,*) '-----------------------------------------------------------------------------------' write(*,*) - deallocate(Om1,X1,Y1,rho1,Om2,X2,Y2,rho2) + deallocate(Om1,X1,Y1,rho1,Om2,X2,Y2,rho2,Bpp,Cpp,Dpp) end if @@ -108,7 +115,8 @@ subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) nVV = nV*(nV-1)/2 allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),rho1(nBas,nBas,nVV), & - Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO),rho2(nBas,nBas,nOO)) + Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO),rho2(nBas,nBas,nOO), & + Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)) write(*,*) '--------------' write(*,*) 'Triplet states' @@ -123,11 +131,13 @@ subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) lambda = rAC(iAC) - ! Initialize T matrix + if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp) + call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp) + call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp) - call ppLR(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Om1,X1,Y1,Om2,X2,Y2,EcAC(ispin)) + call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcAC(ispin)) - call ppACFDT_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOO,nVV,X1,Y1,X2,Y2,Ec(iAC,ispin)) + call ppACFDT_correlation_energy(ispin,nBas,nC,nO,nV,nR,ERI,nOO,nVV,X1,Y1,X2,Y2,Ec(iAC,ispin)) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) @@ -140,7 +150,7 @@ subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) write(*,*) '-----------------------------------------------------------------------------------' write(*,*) - deallocate(Om1,X1,Y1,rho1,Om2,X2,Y2,rho2) + deallocate(Om1,X1,Y1,rho1,Om2,X2,Y2,rho2,Bpp,Cpp,Dpp) end if diff --git a/src/RPA/ppACFDT_correlation_energy.f90 b/src/RPA/ppACFDT_correlation_energy.f90 index 18a9185..2f26fc1 100644 --- a/src/RPA/ppACFDT_correlation_energy.f90 +++ b/src/RPA/ppACFDT_correlation_energy.f90 @@ -1,4 +1,4 @@ -subroutine ppACFDT_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOO,nVV,X1,Y1,X2,Y2,EcAC) +subroutine ppACFDT_correlation_energy(ispin,nBas,nC,nO,nV,nR,ERI,nOO,nVV,X1,Y1,X2,Y2,EcAC) ! Compute the correlation energy via the adiabatic connection formula for the pp sector @@ -13,7 +13,6 @@ subroutine ppACFDT_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOO,nVV,X1,Y integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR - integer,intent(in) :: nS double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) integer,intent(in) :: nOO @@ -53,4 +52,3 @@ subroutine ppACFDT_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOO,nVV,X1,Y - trace_matrix(nVV,C) - trace_matrix(nOO,D) end subroutine - diff --git a/src/RPA/ppRPA.f90 b/src/RPA/ppRPA.f90 index 9c25bf4..555d964 100644 --- a/src/RPA/ppRPA.f90 +++ b/src/RPA/ppRPA.f90 @@ -1,6 +1,6 @@ -subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipole_int,e) +subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e) -! Perform pp-RPA calculation +! Perform ppRPA calculation implicit none include 'parameters.h' @@ -17,7 +17,7 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipo integer,intent(in) :: nV integer,intent(in) :: nR double precision,intent(in) :: ENuc - double precision,intent(in) :: ERHF + double precision,intent(in) :: EHF double precision,intent(in) :: e(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) @@ -25,9 +25,11 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipo ! Local variables integer :: ispin - integer :: nS integer :: nOO integer :: nVV + double precision,allocatable :: Bpp(:,:) + double precision,allocatable :: Cpp(:,:) + double precision,allocatable :: Dpp(:,:) double precision,allocatable :: Om1(:) double precision,allocatable :: X1(:,:) double precision,allocatable :: Y1(:,:) @@ -35,7 +37,7 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipo double precision,allocatable :: X2(:,:) double precision,allocatable :: Y2(:,:) - double precision :: Ec_ppRPA(nspin) + double precision :: EcRPA(nspin) double precision :: EcAC(nspin) ! Hello world @@ -48,12 +50,8 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipo ! Initialization - Ec_ppRPA(:) = 0d0 - EcAC(:) = 0d0 - -! Useful quantities - - nS = nO*nV + EcRPA(:) = 0d0 + EcAC(:) = 0d0 ! Singlet manifold @@ -69,16 +67,21 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipo nOO = nO*(nO+1)/2 nVV = nV*(nV+1)/2 - allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO)) + allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & + Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)) - call ppLR(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,e,ERI,Om1,X1,Y1,Om2,X2,Y2,Ec_ppRPA(ispin)) + if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) + call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,e,ERI,Cpp) + call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,e,ERI,Dpp) + + call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA(ispin)) ! call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) - call print_excitation('pp-BSE (N+2)',ispin,nVV,Om1) - call print_excitation('pp-BSE (N-2)',ispin,nOO,Om2) + call print_excitation('ppRPA (N+2) ',ispin,nVV,Om1) + call print_excitation('ppRPA (N-2) ',ispin,nOO,Om2) - deallocate(Om1,X1,Y1,Om2,X2,Y2) + deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp) endif @@ -96,25 +99,30 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipo nOO = nO*(nO-1)/2 nVV = nV*(nV-1)/2 - allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO)) + allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & + Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)) - call ppLR(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,e,ERI,Om1,X1,Y1,Om2,X2,Y2,Ec_ppRPA(ispin)) + if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) + call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,e,ERI,Cpp) + call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,e,ERI,Dpp) + + call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA(ispin)) ! call print_transition_vectors_pp(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) - call print_excitation('pp-BSE (N+2)',ispin,nVV,Om1) - call print_excitation('pp-BSE (N-2)',ispin,nOO,Om2) + call print_excitation('ppRPA (N+2) ',ispin,nVV,Om1) + call print_excitation('ppRPA (N-2) ',ispin,nOO,Om2) - deallocate(Om1,X1,Y1,Om2,X2,Y2) + deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp) endif write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (singlet) =',Ec_ppRPA(1) - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (triplet) =',3d0*Ec_ppRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy =',Ec_ppRPA(1) + 3d0*Ec_ppRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA total energy =',ENuc + ERHF + Ec_ppRPA(1) + 3d0*Ec_ppRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (singlet) =',EcRPA(1) + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (triplet) =',3d0*EcRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy =',EcRPA(1) + 3d0*EcRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA total energy =',ENuc + EHF + EcRPA(1) + 3d0*EcRPA(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) @@ -122,22 +130,22 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipo if(doACFDT) then - write(*,*) '---------------------------------------------------------' - write(*,*) 'Adiabatic connection version of pp-RPA correlation energy' - write(*,*) '---------------------------------------------------------' + write(*,*) '--------------------------------------------------------' + write(*,*) 'Adiabatic connection version of ppRPA correlation energy' + write(*,*) '--------------------------------------------------------' write(*,*) - call ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) + call ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ERI,e,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (singlet) =',EcAC(1),' au' write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (triplet) =',EcAC(2),' au' write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy =',EcAC(1) + EcAC(2),' au' - write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA total energy =',ENuc + ERHF + EcAC(1) + EcAC(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA total energy =',ENuc + EHF + EcAC(1) + EcAC(2),' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) end if -end subroutine ppRPA +end subroutine