diff --git a/src/CC/RCC.f90 b/src/CC/RCC.f90 index 142d8fc..fd271cc 100644 --- a/src/CC/RCC.f90 +++ b/src/CC/RCC.f90 @@ -25,10 +25,10 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d double precision,intent(in) :: thresh integer,intent(in) :: nBas - integer,intent(in) :: nC(nspin) - integer,intent(in) :: nO(nspin) - integer,intent(in) :: nV(nspin) - integer,intent(in) :: nR(nspin) + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF double precision,intent(in) :: eHF(nBas) diff --git a/src/CI/RCI.f90 b/src/CI/RCI.f90 index 84a7fe8..3762baf 100644 --- a/src/CI/RCI.f90 +++ b/src/CI/RCI.f90 @@ -19,11 +19,11 @@ subroutine RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,n logical,intent(in) :: singlet logical,intent(in) :: triplet integer,intent(in) :: nBas - integer,intent(in) :: nC(nspin) - integer,intent(in) :: nO(nspin) - integer,intent(in) :: nV(nspin) - integer,intent(in) :: nR(nspin) - integer,intent(in) :: nS(nspin) + 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(nBas) double precision,intent(in) :: cHF(nBas,nBas) diff --git a/src/GF/GGF.f90 b/src/GF/GGF.f90 index c8064ae..cb0d492 100644 --- a/src/GF/GGF.f90 +++ b/src/GF/GGF.f90 @@ -32,11 +32,11 @@ subroutine GGF(dotest,doG0F2,doevGF2,doqsGF2,maxSCF,thresh,max_diis,dophBSE,dopp double precision,intent(in) :: ENuc integer,intent(in) :: nBas - integer,intent(in) :: nC(nspin) - integer,intent(in) :: nO(nspin) - integer,intent(in) :: nV(nspin) - integer,intent(in) :: nR(nspin) - integer,intent(in) :: nS(nspin) + 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(nBas) diff --git a/src/GF/RGF.f90 b/src/GF/RGF.f90 index 490dd49..db4c2af 100644 --- a/src/GF/RGF.f90 +++ b/src/GF/RGF.f90 @@ -39,11 +39,11 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm,maxSCF,thresh double precision,intent(in) :: ENuc integer,intent(in) :: nBas - integer,intent(in) :: nC(nspin) - integer,intent(in) :: nO(nspin) - integer,intent(in) :: nV(nspin) - integer,intent(in) :: nR(nspin) - integer,intent(in) :: nS(nspin) + 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(nBas) diff --git a/src/GT/GTpp_QP_graph.f90 b/src/GT/GTpp_QP_graph.f90 index bc53745..ff36fb0 100644 --- a/src/GT/GTpp_QP_graph.f90 +++ b/src/GT/GTpp_QP_graph.f90 @@ -42,9 +42,6 @@ subroutine GTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s double precision,intent(out) :: eGT(nBas) double precision,intent(out) :: Z(nBas) - SigC = 0d0 - dSigC = 0d0 - ! Run Newton's algorithm to find the root write(*,*)'-----------------------------------------------------' diff --git a/src/GT/RG0T0pp.f90 b/src/GT/RG0T0pp.f90 index b1158e6..43d35d2 100644 --- a/src/GT/RG0T0pp.f90 +++ b/src/GT/RG0T0pp.f90 @@ -39,6 +39,8 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d ! Local variables + logical :: print_T = .false. + integer :: ispin integer :: iblock integer :: nOOs,nOOt @@ -76,7 +78,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d ! TDA for T if(TDA_T) then - write(*,*) 'Tamm-Dancoff approximation for pp T-matrix!' + write(*,*) 'Tamm-Dancoff approximation activated for pp T-matrix!' write(*,*) end if @@ -128,8 +130,8 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d deallocate(Bpp,Cpp,Dpp) - call print_excitation_energies('ppRPA@RHF','2p (alpha-beta)',nVVs,Om1s(:)) - call print_excitation_energies('ppRPA@RHF','2h (alpha-beta)',nOOs,Om2s(:)) + if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-beta)',nVVs,Om1s) + if(print_T) call print_excitation_energies('ppRPA@RHF','2h (alpha-beta)',nOOs,Om2s) !---------------------------------------------- ! alpha-alpha block @@ -151,8 +153,8 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d deallocate(Bpp,Cpp,Dpp) - call print_excitation_energies('ppRPA@RHF','2p (alpha-alpha)',nVVt,Om1t) - call print_excitation_energies('ppRPA@RHF','2h (alpha-beta)',nOOt,Om2t) + if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-alpha)',nVVt,Om1t) + if(print_T) call print_excitation_energies('ppRPA@RHF','2h (alpha-alpha)',nOOt,Om2t) !---------------------------------------------- ! Compute excitation densities @@ -214,32 +216,32 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d ! iblock = 1 iblock = 3 - allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) +! 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,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) - 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(ispin)) - deallocate(Bpp,Cpp,Dpp) +! deallocate(Bpp,Cpp,Dpp) ispin = 2 ! iblock = 2 iblock = 4 - allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) +! 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,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) - 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(ispin)) - deallocate(Bpp,Cpp,Dpp) +! deallocate(Bpp,Cpp,Dpp) - EcRPA(1) = EcRPA(1) - EcRPA(2) - EcRPA(2) = 3d0*EcRPA(2) +! EcRPA(1) = EcRPA(1) - EcRPA(2) +! EcRPA(2) = 3d0*EcRPA(2) call print_RG0T0pp(nBas,nO,eHF,ENuc,ERHF,Sig,Z,eGT,EcGM,EcRPA) @@ -296,10 +298,10 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp correlation energy (singlet) =',EcBSE(1),' au' - write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp correlation energy (triplet) =',EcBSE(2),' au' - write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp correlation energy =',sum(EcBSE),' au' - write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp total energy =',ENuc + ERHF + sum(EcBSE),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp@RHF correlation energy (singlet) = ',EcBSE(1),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp@RHF correlation energy (triplet) = ',EcBSE(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp@RHF correlation energy = ',sum(EcBSE),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) @@ -315,10 +317,10 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp correlation energy (singlet) =',EcBSE(1),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp correlation energy (triplet) =',EcBSE(2),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp correlation energy =',sum(EcBSE),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp total energy =',ENuc + ERHF + sum(EcBSE),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp@RHF correlation energy (singlet) = ',EcBSE(1),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp@RHF correlation energy (triplet) = ',EcBSE(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp@RHF correlation energy = ',sum(EcBSE),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/GT/RGT.f90 b/src/GT/RGT.f90 index e35890c..363910a 100644 --- a/src/GT/RGT.f90 +++ b/src/GT/RGT.f90 @@ -1,7 +1,7 @@ -subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thresh,max_diis,doACFDT, & - exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & - linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, & - ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF) +subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thresh,max_diis, & + doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & + linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, & + ERI_AO,ERI_MO,dipole_int_AO,dipole_int,PHF,cHF,eHF) ! T-matrix module @@ -44,14 +44,14 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxS double precision,intent(in) :: ENuc integer,intent(in) :: nBas - integer,intent(in) :: nC(nspin) - integer,intent(in) :: nO(nspin) - integer,intent(in) :: nV(nspin) - integer,intent(in) :: nR(nspin) - integer,intent(in) :: nS(nspin) + 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(nBas) + double precision,intent(in) :: ERHF + double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: S(nBas,nBas) @@ -60,7 +60,7 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxS double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: X(nBas,nBas) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) @@ -77,13 +77,21 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxS call wall_time(start_GT) call RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & - linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) + linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_GT,' seconds' write(*,*) +! call wall_time(start_GT) +! call ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) +! call wall_time(end_GT) +! +! t_GT = end_GT - start_GT +! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0T0 = ',t_GT,' seconds' +! write(*,*) + end if !------------------------------------------------------------------------ @@ -94,7 +102,7 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxS 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,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) + linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT @@ -110,9 +118,9 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxS 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,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int, & - PHF,cHF,epsHF) + call qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & + eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int, & + PHF,cHF,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT @@ -129,7 +137,7 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxS 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,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) + linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT @@ -146,7 +154,7 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxS 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,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) + singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT @@ -163,8 +171,8 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxS 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,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int, & - PHF,cHF,epsHF) + eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int, & + PHF,cHF,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT diff --git a/src/GT/print_RG0T0pp.f90 b/src/GT/print_RG0T0pp.f90 index 1028c99..de5e911 100644 --- a/src/GT/print_RG0T0pp.f90 +++ b/src/GT/print_RG0T0pp.f90 @@ -36,7 +36,7 @@ subroutine print_RG0T0pp(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) ! Dump results write(*,*)'-------------------------------------------------------------------------------' - write(*,*)' One-shot G0T0pp calculation ' + write(*,*)' G0T0pp@RHF calculation ' write(*,*)'-------------------------------------------------------------------------------' write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & '|','#','|','e_HF (eV)','|','Sig_GTpp (eV)','|','Z','|','e_GTpp (eV)','|' @@ -48,16 +48,16 @@ subroutine print_RG0T0pp(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) enddo write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A60,F15.6,A3)') 'RG0T0pp HOMO energy =',eGT(HOMO)*HaToeV,' eV' - write(*,'(2X,A60,F15.6,A3)') 'RG0T0pp LUMO energy =',eGT(LUMO)*HaToeV,' eV' - write(*,'(2X,A60,F15.6,A3)') 'RG0T0pp HOMO-LUMO gap =',Gap*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@RHF HOMO energy = ',eGT(HOMO)*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@RHF LUMO energy = ',eGT(LUMO)*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@RHF HOMO-LUMO gap = ',Gap*HaToeV,' eV' write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A60,F15.6,A3)') 'ppRPA@RG0T0pp correlation energy (singlet) =',EcRPA(1),' au' - write(*,'(2X,A60,F15.6,A3)') 'ppRPA@RG0T0pp correlation energy (triplet) =',EcRPA(2),' au' - write(*,'(2X,A60,F15.6,A3)') 'ppRPA@RG0T0pp correlation energy =',EcRPA(1) + EcRPA(2),' au' - write(*,'(2X,A60,F15.6,A3)') 'ppRPA@RG0T0pp total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2),' au' - write(*,'(2X,A60,F15.6,A3)') ' GM@RG0T0pp correlation energy =',EcGM,' au' - write(*,'(2X,A60,F15.6,A3)') ' GM@RG0T0pp total energy =',ENuc + ERHF + EcGM,' au' + write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp@RHF correlation energy (singlet) = ',EcRPA(1),' au' + write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp@RHF correlation energy (triplet) = ',EcRPA(2),' au' + write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp@RHF correlation energy = ',sum(EcRPA),' au' + write(*,'(2X,A60,F15.6,A3)') 'ppRPA@G0T0pp@RHF total energy = ',ENuc + ERHF + sum(EcRPA),' au' + write(*,'(2X,A60,F15.6,A3)') ' GM@G0T0pp@RHF correlation energy = ',EcGM,' au' + write(*,'(2X,A60,F15.6,A3)') ' GM@G0T0pp@RHF total energy = ',ENuc + ERHF + EcGM,' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/GT/ufG0T0pp.f90 b/src/GT/ufG0T0pp.f90 index 1936d58..cb60511 100644 --- a/src/GT/ufG0T0pp.f90 +++ b/src/GT/ufG0T0pp.f90 @@ -1,4 +1,4 @@ -subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) +subroutine ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Upfolded G0T0pp equations @@ -9,7 +9,7 @@ subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) logical,intent(in) :: dotest - logical,intent(in) :: TDA_W + logical,intent(in) :: TDA_T integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO @@ -27,23 +27,31 @@ subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) integer :: s integer :: i,j,k,l integer :: a,b,c,d - integer :: jb,kc,ia,ja + integer :: ij,ab integer :: klc,kcd,ija,ijb,iab,jab logical :: dRPA integer :: ispin - double precision :: EcRPA + integer :: iblock + integer :: nOO,nOOs,nOOt + integer :: nVV,nVVs,nVVt + double precision :: EcRPA(nspin) integer :: n2h1p,n2p1h,nH double precision,external :: Kronecker_delta double precision,allocatable :: H(:,:) - double precision,allocatable :: eGW(:) + double precision,allocatable :: eGT(:) double precision,allocatable :: Z(:) - double precision,allocatable :: Aph(:,:) - double precision,allocatable :: Bph(:,:) - double precision,allocatable :: Om(:) - double precision,allocatable :: XpY(:,:) - double precision,allocatable :: XmY(:,:) - double precision,allocatable :: rho(:,:,:) + double precision,allocatable :: Bpp(:,:) + double precision,allocatable :: Cpp(:,:) + double precision,allocatable :: Dpp(:,:) + double precision,allocatable :: Om1s(:),Om1t(:) + double precision,allocatable :: X1s(:,:),X1t(:,:) + double precision,allocatable :: Y1s(:,:),Y1t(:,:) + double precision,allocatable :: rho1s(:,:,:),rho1t(:,:,:) + double precision,allocatable :: Om2s(:),Om2t(:) + double precision,allocatable :: X2s(:,:),X2t(:,:) + double precision,allocatable :: Y2s(:,:),Y2t(:,:) + double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:) logical :: verbose = .true. double precision,parameter :: cutoff1 = 0.01d0 @@ -63,15 +71,25 @@ subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) write(*,*)'******************************************' write(*,*) +! Dimensions of the ppRPA linear reponse matrices + + nOOs = nO*(nO + 1)/2 + nVVs = nV*(nV + 1)/2 + +! nOOs = nO*nO +! nVVs = nV*nV + + nOOt = nO*(nO - 1)/2 + nVVt = nV*(nV - 1)/2 + + nOO = nO*nO + nVV = nV*nV + ! Dimension of the supermatrix - n2h1p = nO*nO*nV - n2p1h = nV*nV*nO + n2h1p = nOO*nV + n2p1h = nVV*nO nH = 1 + n2h1p + n2p1h - -! Memory allocation - - allocate(H(nH,nH),eGW(nH),Z(nH)) ! Initialization @@ -80,6 +98,10 @@ subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) eF = 0.5d0*(eHF(nO+1) + eHF(nO)) +! Memory allocation + + allocate(H(nH,nH),eGT(nH),Z(nH)) + !-------------------------! ! Main loop over orbitals ! !-------------------------! @@ -88,9 +110,9 @@ subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) H(:,:) = 0d0 - if (TDA_W) then + if (TDA_T) then - ! TDA for W + ! TDA for T write(*,*) 'Tamm-Dancoff approximation actived!' write(*,*) @@ -214,7 +236,7 @@ subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) else - ! RPA for W + ! RPA for T write(*,*) 'Tamm-Dancoff approximation deactivated!' write(*,*) @@ -242,8 +264,9 @@ subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! alpha-beta block - ispin = 1 - iblock = 3 + ispin = 1 + iblock = 1 +! iblock = 3 ! Compute linear response @@ -263,7 +286,8 @@ subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! alpha-alpha block ispin = 2 - iblock = 4 + iblock = 2 +! iblock = 4 ! Compute linear response @@ -290,6 +314,8 @@ subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) iblock = 4 call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) + deallocate(Om1s,X1s,Y1s,Om2s,X2s,Y2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t) + call wall_time(start_timing) !---------! @@ -303,11 +329,20 @@ subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) !-------------! ija = 0 - do i=nC+1,nO - do ja=1,nS + do ij=1,nOOs + do a=nO+1,nBas-nR ija = ija + 1 - H(1+ija,1+ija) = eHF(i) - Om(ja) + H(1+ija,1+ija) = - eHF(i) + Om2s(ij) + + end do + end do + + do ij=1,nOOt + do a=nO+1,nBas-nR + ija = ija + 1 + + H(1+ija,1+ija) = - eHF(i) + Om2t(ij) end do end do @@ -317,12 +352,22 @@ subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) !-------------! ija = 0 - do i=nC+1,nO - do ja=1,nS + do ij=1,nOOs + do a=nO+1,nBas-nR ija = ija + 1 - H(1 ,1+ija) = sqrt(2d0)*rho(p,i,ja) - H(1+ija,1 ) = sqrt(2d0)*rho(p,i,ja) + H(1 ,1+ija) = rho2s(p,a,ij) + H(1+ija,1 ) = rho2s(p,a,ij) + + end do + end do + + do ij=1,nOOt + do a=nO+1,nBas-nR + ija = ija + 1 + + H(1 ,1+ija) = rho2t(p,a,ij) + H(1+ija,1 ) = rho2t(p,a,ij) end do end do @@ -332,11 +377,20 @@ subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) !-------------! iab = 0 - do ia=1,nS - do b=nO+1,nBas-nR + do ab=1,nVVs + do i=nC+1,nO iab = iab + 1 - H(1+n2h1p+iab,1+n2h1p+iab) = eHF(b) + Om(ia) + H(1+n2h1p+iab,1+n2h1p+iab) = - eHF(i) + Om1s(ab) + + end do + end do + + do ab=1,nVVt + do i=nC+1,nO + iab = iab + 1 + + H(1+n2h1p+iab,1+n2h1p+iab) = - eHF(i) + Om2s(ab) end do end do @@ -346,19 +400,29 @@ subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) !-------------! iab = 0 - do ia=1,nS - do b=nO+1,nBas-nR + do ab=1,nVVs + do i=nC+1,nO iab = iab + 1 - H(1 ,1+n2h1p+iab) = sqrt(2d0)*rho(p,b,ia) - H(1+n2h1p+iab,1 ) = sqrt(2d0)*rho(p,b,ia) + H(1 ,1+n2h1p+iab) = rho1s(p,i,ab) + H(1+n2h1p+iab,1 ) = rho1s(p,i,ab) + + end do + end do + + do ab=1,nVVt + do i=nC+1,nO + iab = iab + 1 + + H(1 ,1+n2h1p+iab) = rho1t(p,i,ab) + H(1+n2h1p+iab,1 ) = rho1t(p,i,ab) end do end do ! Memory deallocation - deallocate(Om,Aph,Bph,XpY,XmY,rho) + deallocate(rho1s,rho2s,rho1t,rho2t) call wall_time(end_timing) @@ -375,7 +439,7 @@ subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) call wall_time(start_timing) - call diagonalize_matrix(nH,H,eGW) + call diagonalize_matrix(nH,H,eGT) call wall_time(end_timing) @@ -397,17 +461,17 @@ subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) !--------------! write(*,*)'-------------------------------------------' - write(*,'(1X,A32,I3,A8)')'| G0W0 energies (eV) for orbital',p,' |' + write(*,'(1X,A32,I3,A8)')'| G0T0pp energies (eV) for orbital',p,' |' write(*,*)'-------------------------------------------' write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') & '|','#','|','e_QP','|','Z','|' write(*,*)'-------------------------------------------' do s=1,nH - if(eGW(s) < eF .and. eGW(s) > eF - window) then + if(eGT(s) < eF .and. eGT(s) > eF - window) then ! if(Z(s) > cutoff1) then write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & - '|',s,'|',eGW(s)*HaToeV,'|',Z(s),'|' + '|',s,'|',eGT(s)*HaToeV,'|',Z(s),'|' end if end do @@ -416,63 +480,63 @@ subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) if(verbose) then - if(TDA_W) then + if(TDA_T) then ! TDA printing format - do s=1,nH - - if(eGW(s) < eF .and. eGW(s) > eF - window) then - - write(*,*)'-------------------------------------------------------------' - write(*,'(1X,A7,1X,I3,A6,I3,A1,1X,A7,F12.6,A13,F6.4,1X)') & - 'Orbital',p,' and #',s,':','e_QP = ',eGW(s)*HaToeV,' eV and Z = ',Z(s) - write(*,*)'-------------------------------------------------------------' - write(*,'(1X,A20,1X,A20,1X,A15,1X)') & - ' Configuration ',' Coefficient ',' Weight ' - write(*,*)'-------------------------------------------------------------' - - if(p <= nO) & - write(*,'(1X,A7,I3,A16,1X,F15.6,1X,F15.6)') & - ' (',p,') ',H(1,s),H(1,s)**2 - if(p > nO) & - write(*,'(1X,A16,I3,A7,1X,F15.6,1X,F15.6)') & - ' (',p,') ',H(1,s),H(1,s)**2 - - ija = 0 - do i=nC+1,nO - do j=nC+1,nO - do a=nO+1,nBas-nR - ija = ija + 1 - - if(abs(H(1+ija,s)) > cutoff2) & - write(*,'(1X,A3,I3,A1,I3,A6,I3,A7,1X,F15.6,1X,F15.6)') & - ' (',i,',',j,') -> (',a,') ',H(1+ija,s),H(1+ija,s)**2 - - end do - end do - end do - - iab = 0 - do i=nC+1,nO - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR - iab = iab + 1 - - if(abs(H(1+n2h1p+iab,s)) > cutoff2) & - write(*,'(1X,A7,I3,A6,I3,A1,I3,A3,1X,F15.6,1X,F15.6)') & - ' (',i,') -> (',a,',',b,') ',H(1+n2h1p+iab,s),H(1+n2h1p+iab,s)**2 - - end do - end do - end do +! do s=1,nH +! +! if(eGT(s) < eF .and. eGT(s) > eF - window) then +! +! write(*,*)'-------------------------------------------------------------' +! write(*,'(1X,A7,1X,I3,A6,I3,A1,1X,A7,F12.6,A13,F6.4,1X)') & +! 'Orbital',p,' and #',s,':','e_QP = ',eGT(s)*HaToeV,' eV and Z = ',Z(s) +! write(*,*)'-------------------------------------------------------------' +! write(*,'(1X,A20,1X,A20,1X,A15,1X)') & +! ' Configuration ',' Coefficient ',' Weight ' +! write(*,*)'-------------------------------------------------------------' +! +! if(p <= nO) & +! write(*,'(1X,A7,I3,A16,1X,F15.6,1X,F15.6)') & +! ' (',p,') ',H(1,s),H(1,s)**2 +! if(p > nO) & +! write(*,'(1X,A16,I3,A7,1X,F15.6,1X,F15.6)') & +! ' (',p,') ',H(1,s),H(1,s)**2 +! +! ija = 0 +! do i=nC+1,nO +! do j=nC+1,nO +! do a=nO+1,nBas-nR +! ija = ija + 1 +! +! if(abs(H(1+ija,s)) > cutoff2) & +! write(*,'(1X,A3,I3,A1,I3,A6,I3,A7,1X,F15.6,1X,F15.6)') & +! ' (',i,',',j,') -> (',a,') ',H(1+ija,s),H(1+ija,s)**2 +! +! end do +! end do +! end do +! +! iab = 0 +! do i=nC+1,nO +! do a=nO+1,nBas-nR +! do b=nO+1,nBas-nR +! iab = iab + 1 +! +! if(abs(H(1+n2h1p+iab,s)) > cutoff2) & +! write(*,'(1X,A7,I3,A6,I3,A1,I3,A3,1X,F15.6,1X,F15.6)') & +! ' (',i,') -> (',a,',',b,') ',H(1+n2h1p+iab,s),H(1+n2h1p+iab,s)**2 +! +! end do +! end do +! end do - write(*,*)'-------------------------------------------------------------' - write(*,*) +! write(*,*)'-------------------------------------------------------------' +! write(*,*) - end if +! end if - end do +! end do else @@ -480,14 +544,14 @@ subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) do s=1,nH - if(eGW(s) < eF .and. eGW(s) > eF - window) then + if(eGT(s) < eF .and. eGT(s) > eF - window) then write(*,*)'-------------------------------------------------------------' write(*,'(1X,A7,1X,I3,A6,I3,A1,1X,A7,F12.6,A13,F6.4,1X)') & - 'Orbital',p,' and #',s,':','e_QP = ',eGW(s)*HaToeV,' eV and Z = ',Z(s) + 'Orbital',p,' and #',s,':','e_QP = ',eGT(s)*HaToeV,' eV and Z = ',Z(s) write(*,*)'-------------------------------------------------------------' write(*,'(1X,A20,1X,A20,1X,A15,1X)') & - ' Conf. (p,ia) ',' Coefficient ',' Weight ' + ' Conf. (i,ab) or (a,ij) ',' Coefficient ',' Weight ' write(*,*)'-------------------------------------------------------------' if(p <= nO) & @@ -498,25 +562,25 @@ subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ' (',p,') ',H(1,s),H(1,s)**2 ija = 0 - do i=nC+1,nO - do ja=1,nS + do ij=1,nOO + do a=nO+1,nBas-nR ija = ija + 1 if(abs(H(1+ija,s)) > cutoff2) & write(*,'(1X,A7,I3,A1,I3,A12,1X,F15.6,1X,F15.6)') & - ' (',i,',',ja,') ',H(1+ija,s),H(1+ija,s)**2 + ' (',a,',',ij,') ',H(1+ija,s),H(1+ija,s)**2 end do end do iab = 0 - do ia=1,nS - do b=nO+1,nBas-nR + do ab=1,nVV + do i=nC+1,nO iab = iab + 1 if(abs(H(1+n2h1p+iab,s)) > cutoff2) & write(*,'(1X,A7,I3,A1,I3,A12,1X,F15.6,1X,F15.6)') & - ' (',ia,',',b,') ',H(1+n2h1p+iab,s),H(1+n2h1p+iab,s)**2 + ' (',i,',',ab,') ',H(1+n2h1p+iab,s),H(1+n2h1p+iab,s)**2 end do end do @@ -534,4 +598,6 @@ subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) end do + deallocate(H,eGT,Z) + end subroutine diff --git a/src/GW/GW_QP_graph.f90 b/src/GW/GW_QP_graph.f90 index d4e77a8..cf240fd 100644 --- a/src/GW/GW_QP_graph.f90 +++ b/src/GW/GW_QP_graph.f90 @@ -58,7 +58,6 @@ subroutine GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eOld,eGW,Z) dSigC = GW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho) f = w - eHF(p) - SigC df = 1d0/(1d0 - dSigC) - w = w - df*f end do @@ -78,6 +77,7 @@ subroutine GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eOld,eGW,Z) end if end do + write(*,*)'-----------------------------------------------------' write(*,*) diff --git a/src/GW/RGW.f90 b/src/GW/RGW.f90 index c137e74..696039d 100644 --- a/src/GW/RGW.f90 +++ b/src/GW/RGW.f90 @@ -44,11 +44,11 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre double precision,intent(in) :: ENuc integer,intent(in) :: nBas - integer,intent(in) :: nC(nspin) - integer,intent(in) :: nO(nspin) - integer,intent(in) :: nV(nspin) - integer,intent(in) :: nR(nspin) - integer,intent(in) :: nS(nspin) + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS double precision,intent(in) :: ERHF double precision,intent(in) :: eHF(nBas) diff --git a/src/HF/print_RHF.f90 b/src/HF/print_RHF.f90 index 89b93ec..3a06d31 100644 --- a/src/HF/print_RHF.f90 +++ b/src/HF/print_RHF.f90 @@ -27,7 +27,7 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole) double precision :: Gap double precision :: S,S2 - logical :: dump_orb = .true. + logical :: dump_orb = .false. ! HOMO and LUMO diff --git a/src/LR/ppLR.f90 b/src/LR/ppLR.f90 index 9ce816a..2e54490 100644 --- a/src/LR/ppLR.f90 +++ b/src/LR/ppLR.f90 @@ -1,6 +1,6 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) -! Compute the pp channel of the linear response: see Scuseria et al. JCP 139, 104113 (2013) +! Solve the pp-RPA linear eigenvalue problem implicit none include 'parameters.h' @@ -69,6 +69,8 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) M( 1:nVV ,nVV+1:nOO+nVV) = - Bpp(1:nVV,1:nOO) M(nVV+1:nOO+nVV, 1:nVV) = + transpose(Bpp(1:nVV,1:nOO)) +! call matout(nOO,nOO,Dpp) + ! Diagonalize the p-p matrix if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Om,Z) diff --git a/src/MP/RMP.f90 b/src/MP/RMP.f90 index 0a88427..9fba986 100644 --- a/src/MP/RMP.f90 +++ b/src/MP/RMP.f90 @@ -14,10 +14,10 @@ subroutine RMP(dotest,doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) logical,intent(in) :: regularize integer,intent(in) :: nBas - integer,intent(in) :: nC(nspin) - integer,intent(in) :: nO(nspin) - integer,intent(in) :: nV(nspin) - integer,intent(in) :: nR(nspin) + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF double precision,intent(in) :: eHF(nBas) diff --git a/src/RPA/RRPA.f90 b/src/RPA/RRPA.f90 index 281f39d..93bf457 100644 --- a/src/RPA/RRPA.f90 +++ b/src/RPA/RRPA.f90 @@ -21,11 +21,11 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker logical,intent(in) :: singlet logical,intent(in) :: triplet integer,intent(in) :: nBas - integer,intent(in) :: nC(nspin) - integer,intent(in) :: nO(nspin) - integer,intent(in) :: nV(nspin) - integer,intent(in) :: nR(nspin) - integer,intent(in) :: nS(nspin) + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF double precision,intent(in) :: eHF(nBas) diff --git a/src/RPA/ppRRPA.f90 b/src/RPA/ppRRPA.f90 index fc5dd77..9c6a330 100644 --- a/src/RPA/ppRRPA.f90 +++ b/src/RPA/ppRRPA.f90 @@ -1,4 +1,4 @@ -subroutine ppRRPA(dotest,TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e) +subroutine ppRRPA(dotest,TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipole_int,eHF) ! Perform ppRPA calculation @@ -19,8 +19,8 @@ subroutine ppRRPA(dotest,TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,E integer,intent(in) :: nV integer,intent(in) :: nR double precision,intent(in) :: ENuc - double precision,intent(in) :: EHF - double precision,intent(in) :: e(nBas) + double precision,intent(in) :: 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) @@ -71,8 +71,8 @@ subroutine ppRRPA(dotest,TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,E Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)) 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_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eHF,ERI,Cpp) + call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eHF,ERI,Dpp) call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA(ispin)) @@ -103,8 +103,8 @@ subroutine ppRRPA(dotest,TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,E Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)) 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_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eHF,ERI,Cpp) + call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eHF,ERI,Dpp) call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA(ispin)) @@ -121,10 +121,10 @@ subroutine ppRRPA(dotest,TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,E write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10,A3)') 'Tr@ppRRPA correlation energy (singlet) = ',EcRPA(1),'au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@ppRRPA correlation energy (triplet) = ',EcRPA(2),'au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@ppRRPA correlation energy = ',sum(EcRPA),'au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@ppRRPA total energy = ',ENuc + EHF + sum(EcRPA),'au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppRPA@RHF correlation energy (singlet) = ',EcRPA(1),'au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppRPA@RHF correlation energy (triplet) = ',EcRPA(2),'au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppRPA@RHF correlation energy = ',sum(EcRPA),'au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppRPA@RHF total energy = ',ENuc + ERHF + sum(EcRPA),'au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) @@ -137,14 +137,14 @@ subroutine ppRRPA(dotest,TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,E write(*,*) '--------------------------------------------------------' write(*,*) - call ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ERI,e,EcRPA) + call ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ERI,eHF,EcRPA) write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10,A3)') 'AC@ppRRPA correlation energy (singlet) = ',EcRPA(1),' au' - write(*,'(2X,A50,F20.10,A3)') 'AC@ppRRPA correlation energy (triplet) = ',EcRPA(2),' au' - write(*,'(2X,A50,F20.10,A3)') 'AC@ppRRPA correlation energy = ',EcRPA(1) + EcRPA(2),' au' - write(*,'(2X,A50,F20.10,A3)') 'AC@ppRRPA total energy = ',ENuc + EHF + EcRPA(1) + EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA@RHF correlation energy (singlet) = ',EcRPA(1),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA@RHF correlation energy (triplet) = ',EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA@RHF correlation energy = ',EcRPA(1) + EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA@RHF total energy = ',ENuc + ERHF + EcRPA(1) + EcRPA(2),' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/RPA/sort_ppRPA.f90 b/src/RPA/sort_ppRPA.f90 index 9b61219..e7037ad 100644 --- a/src/RPA/sort_ppRPA.f90 +++ b/src/RPA/sort_ppRPA.f90 @@ -41,19 +41,17 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) ! Memory allocation - allocate(M(nOO+nVV,nOO+nVV), & - Z1(nOO+nVV,nVV),Z2(nOO+nVV,nOO), & - order1(nVV),order2(nOO)) + allocate(M(nOO+nVV,nOO+nVV),Z1(nOO+nVV,nVV),Z2(nOO+nVV,nOO),order1(nVV),order2(nOO)) ! Initializatiom - Om1(:) = 0d0 - X1(:,:) = 0d0 - Y1(:,:) = 0d0 + Om1(:) = 0d0 + X1(:,:) = 0d0 + Y1(:,:) = 0d0 - Om2(:) = 0d0 - X2(:,:) = 0d0 - Y2(:,:) = 0d0 + Om2(:) = 0d0 + X2(:,:) = 0d0 + Y2(:,:) = 0d0 ! Compute metric @@ -90,9 +88,8 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) end do - if(minval(Om1(:)) < 0d0 .or. ab /= nVV) call print_warning('You may have instabilities in pp-RPA!!') - if(maxval(Om2(:)) > 0d0 .or. ij /= nOO) call print_warning('You may have instabilities in pp-RPA!!') - + if(minval(Om1) < 0d0 .or. ab /= nVV) call print_warning('You may have instabilities in pp-RPA!!') + if(maxval(Om2) > 0d0 .or. ij /= nOO) call print_warning('You may have instabilities in pp-RPA!!') if(nVV > 0) then @@ -100,8 +97,8 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) order1(ab) = ab end do - call quick_sort(Om1(:),order1(:),nVV) - call set_order(Z1(:,:),order1(:),nOO+nVV,nVV) + call quick_sort(Om1,order1,nVV) + call set_order(Z1,order1,nOO+nVV,nVV) end if @@ -111,8 +108,8 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) order2(ij) = ij end do - call quick_sort(Om2(:),order2(:),nOO) - call set_order(Z2(:,:),order2(:),nOO+nVV,nOO) + call quick_sort(Om2,order2,nOO) + call set_order(Z2,order2,nOO+nVV,nOO) end if @@ -205,21 +202,24 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) if(nVV > 0) call dgemm ('N', 'N', nOO+nVV, nVV, nOO+nVV, 1d0, M, nOO+nVV, Z1, nOO+nVV, 0d0, tmp1, nOO+nVV) if(nVV > 0) call dgemm ('T', 'N', nVV , nVV, nOO+nVV, 1d0, Z1, nOO+nVV, tmp1, nOO+nVV, 0d0, S1, nVV) - !S1 = + matmul(transpose(Z1),matmul(M,Z1)) if(nOO > 0) call dgemm ('N', 'N', nOO+nVV, nOO, nOO+nVV, 1d0, M, nOO+nVV, -1d0*Z2, nOO+nVV, 0d0, tmp2, nOO+nVV) if(nOO > 0) call dgemm ('T', 'N', nOO , nOO, nOO+nVV, 1d0, Z2, nOO+nVV, tmp2, nOO+nVV, 0d0, S2, nOO) - ! S2 = - matmul(transpose(Z2),matmul(M,Z2)) - if(nVV > 0) call orthogonalization_matrix(1,nVV,S1,O1) - if(nOO > 0) call orthogonalization_matrix(1,nOO,S2,O2) +! S1 = + matmul(transpose(Z1),matmul(M,Z1)) +! S2 = - matmul(transpose(Z2),matmul(M,Z2)) + + if(nVV > 0) call orthogonalization_matrix(nVV,S1,O1) + if(nOO > 0) call orthogonalization_matrix(nOO,S2,O2) - !Z1 = matmul(Z1,O1) if(nVV > 0) call dgemm ('N', 'N', nOO+nVV,nVV,nVV, 1d0, Z1, nOO+nVV, O1, nVV,0d0, tmp1, nOO+nVV) Z1 = tmp1 if(nOO > 0) call dgemm ('N', 'N', nOO+nVV,nOO,nOO, 1d0, Z2, nOO+nVV, O2, nOO,0d0, tmp2, nOO+nVV) Z2 = tmp2 +! Z1 = matmul(Z1,O1) +! Z2 = matmul(Z2,O2) + ! Define submatrices X1, Y1, X2, & Y2 X1(1:nVV,1:nVV) = + Z1( 1: nVV,1:nVV) diff --git a/src/utils/wrap_lapack.f90 b/src/utils/wrap_lapack.f90 index c8f3797..3801928 100644 --- a/src/utils/wrap_lapack.f90 +++ b/src/utils/wrap_lapack.f90 @@ -178,6 +178,7 @@ subroutine linear_solve(N,A,b,x,rcond) ! Find optimal size for temporary arrays allocate(work(1)) + allocate(AF(N,N),ipiv(N),iwork(N)) lwork = -1 call dsysvx('N','U',N,1,A,N,AF,N,ipiv,b,N,x,N,rcond,ferr,berr,work,lwork,iwork,info) @@ -185,7 +186,7 @@ subroutine linear_solve(N,A,b,x,rcond) deallocate(work) - allocate(AF(N,N),ipiv(N),work(lwork),iwork(N)) + allocate(work(lwork)) call dsysvx('N','U',N,1,A,N,AF,N,ipiv,b,N,x,N,rcond,ferr,berr,work,lwork,iwork,info)