diff --git a/input/methods b/input/methods index 91ecc0d..6afd5a5 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 T + F F T F # 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 dc7e3fb..0ebb0db 100644 --- a/input/options +++ b/input/options @@ -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 - T F F + T T F # BSE: phBSE phBSE2 ppBSE dBSE dTDA evDyn F F F T F F diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90 index fdff910..1ddfbe2 100644 --- a/src/GW/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -229,7 +229,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD end if - call phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF,eGW,EcAC) + call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF,eGW,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GW/SRG_qsGW.f90 b/src/GW/SRG_qsGW.f90 index 80104e7..a574ffc 100644 --- a/src/GW/SRG_qsGW.f90 +++ b/src/GW/SRG_qsGW.f90 @@ -342,7 +342,7 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE end if - call phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC) + call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GW/evGW.f90 b/src/GW/evGW.f90 index 1547f0f..658c255 100644 --- a/src/GW/evGW.f90 +++ b/src/GW/evGW.f90 @@ -290,7 +290,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, end if - call phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC) + call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GW/qsGW.f90 b/src/GW/qsGW.f90 index 98abeb3..6faa700 100644 --- a/src/GW/qsGW.f90 +++ b/src/GW/qsGW.f90 @@ -331,7 +331,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, end if - call phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC) + call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/LR/phLR.f90 b/src/LR/phLR.f90 index 4b0e48e..9e24cf9 100644 --- a/src/LR/phLR.f90 +++ b/src/LR/phLR.f90 @@ -1,4 +1,4 @@ -subroutine phLR(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Ec,Omega,XpY,XmY) +subroutine phLR(TDA,nS,A,B,EcRPA,Om,XpY,XmY) ! Compute linear response @@ -7,27 +7,16 @@ subroutine phLR(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Ec,Omega,XpY ! Input variables - logical,intent(in) :: dRPA logical,intent(in) :: TDA - double precision,intent(in) :: eta - integer,intent(in) :: ispin - integer,intent(in) :: nBas - integer,intent(in) :: nC - integer,intent(in) :: nO - integer,intent(in) :: nV - integer,intent(in) :: nR integer,intent(in) :: nS - double precision,intent(in) :: lambda - double precision,intent(in) :: e(nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: A(nS,nS) + double precision,intent(in) :: B(nS,nS) ! Local variables integer :: i,j,k double precision :: trace_matrix - double precision,allocatable :: A(:,:) - double precision,allocatable :: B(:,:) double precision,allocatable :: ApB(:,:) double precision,allocatable :: AmB(:,:) double precision,allocatable :: AmBSq(:,:) @@ -37,33 +26,26 @@ subroutine phLR(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Ec,Omega,XpY ! Output variables - double precision,intent(out) :: Ec - double precision,intent(out) :: Omega(nS) + double precision,intent(out) :: EcRPA + double precision,intent(out) :: Om(nS) double precision,intent(out) :: XpY(nS,nS) double precision,intent(out) :: XmY(nS,nS) ! Memory allocation - allocate(A(nS,nS),B(nS,nS),ApB(nS,nS),AmB(nS,nS),AmBSq(nS,nS),AmBIv(nS,nS),Z(nS,nS),tmp(nS,nS)) - -! Build A and B matrices - - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A) + allocate(ApB(nS,nS),AmB(nS,nS),AmBSq(nS,nS),AmBIv(nS,nS),Z(nS,nS),tmp(nS,nS)) ! Tamm-Dancoff approximation if(TDA) then - B(:,:) = 0d0 XpY(:,:) = A(:,:) - call diagonalize_matrix(nS,XpY,Omega) + call diagonalize_matrix(nS,XpY,Om) XpY(:,:) = transpose(XpY(:,:)) XmY(:,:) = XpY(:,:) else - call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B) - ! Build A + B and A - B matrices ApB = A + B @@ -71,42 +53,34 @@ subroutine phLR(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Ec,Omega,XpY ! Diagonalize linear response matrix - call diagonalize_matrix(nS,AmB,Omega) + call diagonalize_matrix(nS,AmB,Om) - if(minval(Omega) < 0d0) & + if(minval(Om) < 0d0) & call print_warning('You may have instabilities in linear response: A-B is not positive definite!!') -! do ia=1,nS -! if(Omega(ia) < 0d0) Omega(ia) = 0d0 -! end do - - call ADAt(nS,AmB,1d0*dsqrt(Omega),AmBSq) - call ADAt(nS,AmB,1d0/dsqrt(Omega),AmBIv) + call ADAt(nS,AmB,1d0*dsqrt(Om),AmBSq) + call ADAt(nS,AmB,1d0/dsqrt(Om),AmBIv) call dgemm('N','N',nS,nS,nS,1d0,ApB,size(ApB,1),AmBSq,size(AmBSq,1),0d0,tmp,size(tmp,1)) call dgemm('N','N',nS,nS,nS,1d0,AmBSq,size(AmBSq,1),tmp,size(tmp,1),0d0,Z,size(Z,1)) - call diagonalize_matrix(nS,Z,Omega) + call diagonalize_matrix(nS,Z,Om) - if(minval(Omega) < 0d0) & + if(minval(Om) < 0d0) & call print_warning('You may have instabilities in linear response: negative excitations!!') - ! do ia=1,nS - ! if(Omega(ia) < 0d0) Omega(ia) = 0d0 - ! end do - - Omega = dsqrt(Omega) + Om = sqrt(Om) call dgemm('T','N',nS,nS,nS,1d0,Z,size(Z,1),AmBSq,size(AmBSq,1),0d0,XpY,size(XpY,1)) - call DA(nS,1d0/dsqrt(Omega),XpY) + call DA(nS,1d0/dsqrt(Om),XpY) call dgemm('T','N',nS,nS,nS,1d0,Z,size(Z,1),AmBIv,size(AmBIv,1),0d0,XmY,size(XmY,1)) - call DA(nS,1d0*dsqrt(Omega),XmY) + call DA(nS,1d0*dsqrt(Om),XmY) end if ! Compute the RPA correlation energy - Ec = 0.5d0*(sum(Omega) - trace_matrix(nS,A)) + EcRPA = 0.5d0*(sum(Om) - trace_matrix(nS,A)) end subroutine diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index c0d6578..02ce218 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -581,7 +581,7 @@ program QuAcK else - call phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF) + call phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF) end if call cpu_time(end_RPA) @@ -606,7 +606,7 @@ program QuAcK else - call phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF) + call phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF) end if call cpu_time(end_RPA) @@ -624,7 +624,7 @@ program QuAcK if(docrRPA) then call cpu_time(start_RPA) - call crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF) + call crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF) call cpu_time(end_RPA) t_RPA = end_RPA - start_RPA diff --git a/src/RPA/crACFDT.f90 b/src/RPA/crACFDT.f90 index c87c4c0..9e82467 100644 --- a/src/RPA/crACFDT.f90 +++ b/src/RPA/crACFDT.f90 @@ -1,4 +1,4 @@ -subroutine crACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,e,EcAC) +subroutine crACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) ! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem ! for the crossed-ring contribution @@ -9,41 +9,31 @@ subroutine crACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta, ! Input variables - logical,intent(in) :: doXBS logical,intent(in) :: exchange_kernel logical,intent(in) :: dRPA - logical,intent(in) :: TDA_W logical,intent(in) :: TDA - logical,intent(in) :: BSE logical,intent(in) :: singlet logical,intent(in) :: triplet - double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nS - double precision,intent(in) :: eW(nBas) double precision,intent(in) :: e(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables integer :: ispin - integer :: isp_W integer :: iAC double precision :: lambda double precision,allocatable :: Ec(:,:) double precision :: EcRPA - double precision,allocatable :: KA(:,:) - double precision,allocatable :: KB(:,:) - double precision,allocatable :: OmRPA(:) - double precision,allocatable :: XpY_RPA(:,:) - double precision,allocatable :: XmY_RPA(:,:) - double precision,allocatable :: rho_RPA(:,:,:) + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) double precision,allocatable :: Om(:) double precision,allocatable :: XpY(:,:) @@ -56,8 +46,7 @@ subroutine crACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta, ! Memory allocation allocate(Ec(nAC,nspin)) - allocate(KA(nS,nS),KB(nS,nS),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) - allocate(Om(nS),XpY(nS,nS),XmY(nS,nS)) + allocate(Aph(nS,nS),Bph(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS)) ! Antisymmetrized kernel version @@ -72,17 +61,6 @@ subroutine crACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta, EcAC(:) = 0d0 Ec(:,:) = 0d0 -! Compute (singlet) RPA screening - - isp_W = 1 - EcRPA = 0d0 - - call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - - call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA) - call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB) - ! Singlet manifold if(singlet) then @@ -102,18 +80,10 @@ subroutine crACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta, lambda = -rAC(iAC) - if(doXBS) then + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph) + if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph) - call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) -! call print_excitation('W^lambda: ',isp_W,nS,OmRPA) - - call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA) - call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB) - - end if - - call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,KA,KB,EcAC(ispin),Om,XpY,XmY) + call phLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY) call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin)) @@ -151,17 +121,10 @@ subroutine crACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta, lambda = -rAC(iAC) - if(doXBS) then + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph) + if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph) - call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - - call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA) - call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB) - - end if - - call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,KA,KB,EcAC(ispin),Om,XpY,XmY) + call phLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY) call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin)) diff --git a/src/RPA/crRPA.f90 b/src/RPA/crRPA.f90 index 197136d..97d72b2 100644 --- a/src/RPA/crRPA.f90 +++ b/src/RPA/crRPA.f90 @@ -1,5 +1,4 @@ -subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & - ERI,dipole_int,eHF) +subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e) ! Crossed-ring channel of the random phase approximation @@ -13,7 +12,6 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n logical,intent(in) :: doACFDT logical,intent(in) :: exchange_kernel logical,intent(in) :: singlet - double precision,intent(in) :: eta logical,intent(in) :: triplet integer,intent(in) :: nBas integer,intent(in) :: nC @@ -22,19 +20,22 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n integer,intent(in) :: nR integer,intent(in) :: nS double precision,intent(in) :: ENuc - double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas) + 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) ! Local variables integer :: ispin - double precision,allocatable :: Om(:,:) - double precision,allocatable :: XpY(:,:,:) - double precision,allocatable :: XmY(:,:,:) + logical :: dRPA + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) + double precision,allocatable :: Om(:) + double precision,allocatable :: XpY(:,:) + double precision,allocatable :: XmY(:,:) - double precision :: EcRPAx(nspin) + double precision :: EcTr(nspin) double precision :: EcAC(nspin) ! Hello world @@ -54,12 +55,15 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n ! Initialization - EcRPAx(:) = 0d0 + dRPA = .false. + + EcTr(:) = 0d0 EcAC(:) = 0d0 ! Memory allocation - allocate(Om(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) + allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS)) + if(.not.TDA) allocate(Bph(nS,nS)) ! Singlet manifold @@ -67,9 +71,12 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n ispin = 1 - call phLR(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,EcRPAx(ispin),Om(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('crRPA@HF ',ispin,nS,Om(:,ispin)) - call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,e,ERI,Aph) + if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,ERI,Bph) + + call phLR(TDA,nS,Aph,Bph,EcTr(ispin),Om,XpY,XmY) + call print_excitation('crRPA@HF ',ispin,nS,Om) + call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) endif @@ -79,25 +86,28 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n ispin = 2 - call phLR(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,EcRPAx(ispin),Om(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('crRPA@HF ',ispin,nS,Om(:,ispin)) - call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,e,ERI,Aph) + if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,ERI,Bph) + + call phLR(TDA,nS,Aph,Bph,EcTr(ispin),Om,XpY,XmY) + call print_excitation('crRPA@HF ',ispin,nS,Om) + call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) endif -! if(exchange_kernel) then + if(exchange_kernel) then - EcRPAx(1) = 0.5d0*EcRPAx(1) - EcRPAx(2) = 1.5d0*EcRPAx(2) + EcTr(1) = 0.5d0*EcTr(1) + EcTr(2) = 1.5d0*EcTr(2) -! end if + end if write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (singlet) =',EcRPAx(1) - write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (triplet) =',EcRPAx(2) - write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy =',EcRPAx(1) + EcRPAx(2) - write(*,'(2X,A50,F20.10)') 'Tr@crRPA total energy =',ENuc + ERHF + EcRPAx(1) + EcRPAx(2) + write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (singlet) =',EcTr(1) + write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (triplet) =',EcTr(2) + write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy =',EcTr(1) + EcTr(2) + write(*,'(2X,A50,F20.10)') 'Tr@crRPA total energy =',ENuc + EHF + EcTr(1) + EcTr(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) @@ -110,15 +120,14 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n write(*,*) '-------------------------------------------------------' write(*,*) - call crACFDT(exchange_kernel,.false.,.false.,.false.,TDA,.false.,singlet,triplet,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC) + call crACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy (singlet) =',EcAC(1) write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy (triplet) =',EcAC(2) write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy =',EcAC(1) + EcAC(2) - write(*,'(2X,A50,F20.10)') 'AC@crRPA total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@crRPA total energy =',ENuc + EHF + EcAC(1) + EcAC(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/RPA/phACFDT.f90 b/src/RPA/phACFDT.f90 index a9d86d1..2cb024c 100644 --- a/src/RPA/phACFDT.f90 +++ b/src/RPA/phACFDT.f90 @@ -1,4 +1,4 @@ -subroutine phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,e,EcAC) +subroutine phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) ! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem @@ -8,41 +8,31 @@ subroutine phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta, ! Input variables - logical,intent(in) :: doXBS logical,intent(in) :: exchange_kernel logical,intent(in) :: dRPA - logical,intent(in) :: TDA_W logical,intent(in) :: TDA - logical,intent(in) :: BSE logical,intent(in) :: singlet logical,intent(in) :: triplet - double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nS - double precision,intent(in) :: eW(nBas) double precision,intent(in) :: e(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables integer :: ispin - integer :: isp_W integer :: iAC double precision :: lambda double precision,allocatable :: Ec(:,:) double precision :: EcRPA - double precision,allocatable :: KA(:,:) - double precision,allocatable :: KB(:,:) - double precision,allocatable :: OmRPA(:) - double precision,allocatable :: XpY_RPA(:,:) - double precision,allocatable :: XmY_RPA(:,:) - double precision,allocatable :: rho_RPA(:,:,:) + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) double precision,allocatable :: Om(:) double precision,allocatable :: XpY(:,:) @@ -55,8 +45,7 @@ subroutine phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta, ! Memory allocation allocate(Ec(nAC,nspin)) - allocate(KA(nS,nS),KB(nS,nS),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) - allocate(Om(nS),XpY(nS,nS),XmY(nS,nS)) + allocate(Aph(nS,nS),Bph(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS)) ! Antisymmetrized kernel version @@ -71,17 +60,6 @@ subroutine phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta, EcAC(:) = 0d0 Ec(:,:) = 0d0 -! Compute (singlet) RPA screening - - isp_W = 1 - EcRPA = 0d0 - - call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - - call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA) - call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB) - ! Singlet manifold if(singlet) then @@ -101,17 +79,10 @@ subroutine phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta, lambda = rAC(iAC) - if(doXBS) then + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph) + if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph) - call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - - call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA) - call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB) - - end if - - call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,KA,KB,EcAC(ispin),Om,XpY,XmY) + call phLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY) call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin)) @@ -149,17 +120,10 @@ subroutine phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta, lambda = rAC(iAC) - if(doXBS) then + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph) + if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph) - call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - - call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA) - call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB) - - end if - - call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,KA,KB,EcAC(ispin),Om,XpY,XmY) + call phLR(TDA,nS,Aph,Bph,EcAc(ispin),Om,XpY,XmY) call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin)) diff --git a/src/RPA/phRPA.f90 b/src/RPA/phRPA.f90 index d4a9766..cd78d14 100644 --- a/src/RPA/phRPA.f90 +++ b/src/RPA/phRPA.f90 @@ -1,4 +1,4 @@ -subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) +subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e) ! Perform a direct random phase approximation calculation @@ -13,7 +13,6 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n logical,intent(in) :: exchange_kernel logical,intent(in) :: singlet logical,intent(in) :: triplet - double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO @@ -21,19 +20,22 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n integer,intent(in) :: nR integer,intent(in) :: nS double precision,intent(in) :: ENuc - double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas) + 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) ! Local variables integer :: ispin + logical :: dRPA + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) double precision,allocatable :: Om(:) double precision,allocatable :: XpY(:,:) double precision,allocatable :: XmY(:,:) - double precision :: EcRPA(nspin) + double precision :: EcTr(nspin) double precision :: EcAC(nspin) ! Hello world @@ -53,12 +55,15 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n ! Initialization - EcRPA(:) = 0d0 + dRPA = .true. + + EcTr(:) = 0d0 EcAC(:) = 0d0 ! Memory allocation - allocate(Om(nS),XpY(nS,nS),XmY(nS,nS)) + allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS)) + if(.not.TDA) allocate(Bph(nS,nS)) ! Singlet manifold @@ -66,8 +71,11 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n ispin = 1 - call phLR(ispin,.true.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,EcRPA(ispin),Om,XpY,XmY) - call print_excitation('RPA@HF ',ispin,nS,Om) + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph) + if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + + call phLR(TDA,nS,Aph,Bph,EcTr(ispin),Om,XpY,XmY) + call print_excitation('phRPA@HF ',ispin,nS,Om) call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) endif @@ -78,60 +86,53 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n ispin = 2 - call phLR(ispin,.true.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,EcRPA(ispin),Om,XpY,XmY) - call print_excitation('RPA@HF ',ispin,nS,Om) + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph) + if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + + call phLR(TDA,nS,Aph,Bph,EcTr(ispin),Om,XpY,XmY) + call print_excitation('phRPA@HF ',ispin,nS,Om) call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) endif -! if(exchange_kernel) then + if(exchange_kernel) then -! EcRPA(1) = 0.5d0*EcRPA(1) -! EcRPA(2) = 1.5d0*EcRPA(2) + EcTr(1) = 0.5d0*EcTr(1) + EcTr(2) = 1.5d0*EcTr(2) -! end if + end if write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy (singlet) =',EcRPA(1) - write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy (triplet) =',EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy =',EcRPA(1) + EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@phRPA total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy (singlet) =',EcTr(1) + write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy (triplet) =',EcTr(2) + write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy =',EcTr(1) + EcTr(2) + write(*,'(2X,A50,F20.10)') 'Tr@phRPA total energy =',ENuc + EHF + EcTr(1) + EcTr(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) - deallocate(Om,XpY,XmY) + deallocate(Om,XpY,XmY,Aph,Bph) ! Compute the correlation energy via the adiabatic connection -! Switch off ACFDT for RPA as the trace formula is equivalent if(doACFDT) then - write(*,*) '------------------------------------------------------' - write(*,*) 'Adiabatic connection version of RPA correlation energy' - write(*,*) '------------------------------------------------------' + write(*,*) '--------------------------------------------------------' + write(*,*) 'Adiabatic connection version of phRPA correlation energy' + write(*,*) '--------------------------------------------------------' write(*,*) - call phACFDT(exchange_kernel,.false.,.true.,.false.,TDA,.false.,singlet,triplet,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC) - - if(exchange_kernel) then - - EcAC(1) = 0.5d0*EcAC(1) - EcAC(2) = 1.5d0*EcAC(2) - - end if + call phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A50,F20.10)') 'AC@phRPA correlation energy (singlet) =',EcAC(1) write(*,'(2X,A50,F20.10)') 'AC@phRPA correlation energy (triplet) =',EcAC(2) write(*,'(2X,A50,F20.10)') 'AC@phRPA correlation energy =',EcAC(1) + EcAC(2) - write(*,'(2X,A50,F20.10)') 'AC@phRPA total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@phRPA total energy =',ENuc + EHF + EcAC(1) + EcAC(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) - end if end subroutine diff --git a/src/RPA/phRPAx.f90 b/src/RPA/phRPAx.f90 index 6fea5b7..24b694b 100644 --- a/src/RPA/phRPAx.f90 +++ b/src/RPA/phRPAx.f90 @@ -1,4 +1,4 @@ -subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) +subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e) ! Perform random phase approximation calculation with exchange (aka TDHF) @@ -13,7 +13,6 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV, logical,intent(in) :: exchange_kernel logical,intent(in) :: singlet logical,intent(in) :: triplet - double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO @@ -21,19 +20,22 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,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) + 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) ! Local variables integer :: ispin + logical :: dRPA + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) double precision,allocatable :: Om(:) double precision,allocatable :: XpY(:,:) double precision,allocatable :: XmY(:,:) - double precision :: EcRPAx(nspin) + double precision :: EcTr(nspin) double precision :: EcAC(nspin) ! Hello world @@ -54,12 +56,15 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV, ! Initialization - EcRPAx(:) = 0d0 - EcAC(:) = 0d0 + dRPA = .false. + + EcTr(:) = 0d0 + EcAC(:) = 0d0 ! Memory allocation - allocate(Om(nS),XpY(nS,nS),XmY(nS,nS)) + allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS)) + if(.not.TDA) allocate(Bph(nS,nS)) ! Singlet manifold @@ -67,8 +72,11 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV, ispin = 1 - call phLR(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,EcRPAx(ispin),Om,XpY,XmY) - call print_excitation('RPAx@HF ',ispin,nS,Om) + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph) + if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + + call phLR(TDA,nS,Aph,Bph,EcTr(ispin),Om,XpY,XmY) + call print_excitation('phRPAx@HF ',ispin,nS,Om) call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) endif @@ -79,50 +87,52 @@ subroutine phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV, ispin = 2 - call phLR(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,EcRPAx(ispin),Om,XpY,XmY) - call print_excitation('RPAx@HF ',ispin,nS,Om) + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph) + if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + + call phLR(TDA,nS,Aph,Bph,EcTr(ispin),Om,XpY,XmY) + call print_excitation('phRPAx@HF ',ispin,nS,Om) call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) endif -! if(exchange_kernel) then + if(exchange_kernel) then - EcRPAx(1) = 0.5d0*EcRPAx(1) - EcRPAx(2) = 1.5d0*EcRPAx(2) + EcTr(1) = 0.5d0*EcTr(1) + EcTr(2) = 1.5d0*EcTr(2) -! end if + end if write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@RPAx correlation energy (singlet) =',EcRPAx(1) - write(*,'(2X,A50,F20.10)') 'Tr@RPAx correlation energy (triplet) =',EcRPAx(2) - write(*,'(2X,A50,F20.10)') 'Tr@RPAx correlation energy =',EcRPAx(1) + EcRPAx(2) - write(*,'(2X,A50,F20.10)') 'Tr@RPAx total energy =',ENuc + ERHF + EcRPAx(1) + EcRPAx(2) + write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy (singlet) =',EcTr(1) + write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy (triplet) =',EcTr(2) + write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy =',EcTr(1) + EcTr(2) + write(*,'(2X,A50,F20.10)') 'Tr@phRPAx total energy =',ENuc + EHF + EcTr(1) + EcTr(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) ! deallocate memory - deallocate(Om,XpY,XmY) + deallocate(Om,XpY,XmY,Aph,Bph) ! Compute the correlation energy via the adiabatic connection if(doACFDT) then - write(*,*) '-------------------------------------------------------' - write(*,*) 'Adiabatic connection version of RPAx correlation energy' - write(*,*) '-------------------------------------------------------' + write(*,*) '---------------------------------------------------------' + write(*,*) 'Adiabatic connection version of phRPAx correlation energy' + write(*,*) '---------------------------------------------------------' write(*,*) - call phACFDT(exchange_kernel,.false.,.false.,.false.,TDA,.false.,singlet,triplet,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC) + call phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy (singlet) =',EcAC(1) - write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy (triplet) =',EcAC(2) - write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy =',EcAC(1) + EcAC(2) - write(*,'(2X,A50,F20.10)') 'AC@RPAx total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@phRPAx correlation energy (singlet) =',EcAC(1) + write(*,'(2X,A50,F20.10)') 'AC@phRPAx correlation energy (triplet) =',EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@phRPAx correlation energy =',EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@phRPAx total energy =',ENuc + EHF + EcAC(1) + EcAC(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*)