diff --git a/src/RPA/ACFDT_Tmatrix.f90 b/src/RPA/ACFDT_Tmatrix.f90 index c51dc1d..5953979 100644 --- a/src/RPA/ACFDT_Tmatrix.f90 +++ b/src/RPA/ACFDT_Tmatrix.f90 @@ -140,7 +140,7 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple call linear_response_BSE(ispin,.false.,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TAt+TAs,TBt+TBs, & EcAC(ispin),Om(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) + call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) @@ -203,7 +203,7 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple call linear_response_BSE(ispin,.false.,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TAt-TAs,TBt-TBs, & EcAC(ispin),Om(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) + call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) diff --git a/src/RPA/ACFDT_cr.f90 b/src/RPA/ACFDT_cr.f90 deleted file mode 100644 index 02f847a..0000000 --- a/src/RPA/ACFDT_cr.f90 +++ /dev/null @@ -1,181 +0,0 @@ -subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,e,EcAC) - -! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem -! for the crossed-ring contribution - - implicit none - include 'parameters.h' - include 'quadrature.h' - -! 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,nC,nO,nV,nR,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 :: WA(:,:) - double precision,allocatable :: WB(:,:) - double precision,allocatable :: OmRPA(:) - double precision,allocatable :: XpY_RPA(:,:) - double precision,allocatable :: XmY_RPA(:,:) - double precision,allocatable :: rho_RPA(:,:,:) - - double precision,allocatable :: Omega(:,:) - double precision,allocatable :: XpY(:,:,:) - double precision,allocatable :: XmY(:,:,:) - -! Output variables - - double precision,intent(out) :: EcAC(nspin) - -! Memory allocation - - allocate(Ec(nAC,nspin)) - allocate(WA(nS,nS),WB(nS,nS),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) - allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) - -! Antisymmetrized kernel version - - if(exchange_kernel) then - - write(*,*) - write(*,*) '*** Exchange kernel version ***' - write(*,*) - - end if - - 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,WA) - call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WB) - -! Singlet manifold - - if(singlet) then - - ispin = 1 - - write(*,*) '--------------' - write(*,*) 'Singlet states' - write(*,*) '--------------' - write(*,*) - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)' - write(*,*) '-----------------------------------------------------------------------------------' - - do iAC=1,nAC - - lambda = -rAC(iAC) - - if(doXBS) then - - 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,WA) - call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB) - - end if - - call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,WA,WB, & - EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - - call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, & - ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) - - write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) - - end do - - EcAC(ispin) = -0.5d0*dot_product(wAC,Ec(:,ispin)) - - if(exchange_kernel) EcAC(ispin) = 0.5d0*EcAC(ispin) - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin) - write(*,*) '-----------------------------------------------------------------------------------' - write(*,*) - - end if - -! Triplet manifold - - if(triplet) then - - ispin = 2 - - write(*,*) '--------------' - write(*,*) 'Triplet states' - write(*,*) '--------------' - write(*,*) - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)' - write(*,*) '-----------------------------------------------------------------------------------' - - do iAC=1,nAC - - lambda = -rAC(iAC) - - if(doXBS) then - - 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,WA) - call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB) - - end if - - call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,WA,WB, & - EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - - call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) - - write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) - - end do - - EcAC(ispin) = -0.5d0*dot_product(wAC,Ec(:,ispin)) - - if(exchange_kernel) EcAC(ispin) = 1.5d0*EcAC(ispin) - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin) - write(*,*) '-----------------------------------------------------------------------------------' - write(*,*) - - end if - -end subroutine ACFDT_cr diff --git a/src/RPA/crRPA.f90 b/src/RPA/crRPA.f90 index 67f1910..197136d 100644 --- a/src/RPA/crRPA.f90 +++ b/src/RPA/crRPA.f90 @@ -30,7 +30,7 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n ! Local variables integer :: ispin - double precision,allocatable :: Omega(:,:) + double precision,allocatable :: Om(:,:) double precision,allocatable :: XpY(:,:,:) double precision,allocatable :: XmY(:,:,:) @@ -59,7 +59,7 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n ! Memory allocation - allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) + allocate(Om(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) ! Singlet manifold @@ -67,9 +67,9 @@ 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),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('crRPA@HF ',ispin,nS,Omega(:,ispin)) - call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + 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)) endif @@ -79,9 +79,9 @@ 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),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('crRPA@HF ',ispin,nS,Omega(:,ispin)) - call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + 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)) endif @@ -110,8 +110,8 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n write(*,*) '-------------------------------------------------------' write(*,*) - call ACFDT_cr(exchange_kernel,.false.,.false.,.false.,TDA,.false.,singlet,triplet,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC) + call crACFDT(exchange_kernel,.false.,.false.,.false.,TDA,.false.,singlet,triplet,eta, & + nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -124,4 +124,4 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n end if -end subroutine crRPA +end subroutine diff --git a/src/RPA/phACFDT.f90 b/src/RPA/phACFDT.f90 index 6a43f89..a9d86d1 100644 --- a/src/RPA/phACFDT.f90 +++ b/src/RPA/phACFDT.f90 @@ -18,7 +18,12 @@ subroutine phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta, logical,intent(in) :: triplet double precision,intent(in) :: eta - integer,intent(in) :: nBas,nC,nO,nV,nR,nS + 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) @@ -32,16 +37,16 @@ subroutine phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta, double precision,allocatable :: Ec(:,:) double precision :: EcRPA - double precision,allocatable :: WA(:,:) - double precision,allocatable :: WB(:,:) + 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 :: Omega(:,:) - double precision,allocatable :: XpY(:,:,:) - double precision,allocatable :: XmY(:,:,:) + double precision,allocatable :: Om(:) + double precision,allocatable :: XpY(:,:) + double precision,allocatable :: XmY(:,:) ! Output variables @@ -50,8 +55,8 @@ subroutine phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta, ! Memory allocation allocate(Ec(nAC,nspin)) - allocate(WA(nS,nS),WB(nS,nS),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) - allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,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)) ! Antisymmetrized kernel version @@ -74,8 +79,8 @@ subroutine phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta, 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,WA) - call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WB) + 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 @@ -100,18 +105,15 @@ subroutine phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta, 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,WA) - call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB) + 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,WA,WB, & - EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + 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 ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, & - ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) + call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin)) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) @@ -152,15 +154,14 @@ subroutine phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta, 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,WA) - call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB) + 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,WA,WB, & - EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + 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 ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) + call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin)) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) diff --git a/src/RPA/ACFDT_correlation_energy.f90 b/src/RPA/phACFDT_correlation_energy.f90 similarity index 77% rename from src/RPA/ACFDT_correlation_energy.f90 rename to src/RPA/phACFDT_correlation_energy.f90 index e320976..48cfaa5 100644 --- a/src/RPA/ACFDT_correlation_energy.f90 +++ b/src/RPA/phACFDT_correlation_energy.f90 @@ -1,4 +1,4 @@ -subroutine ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,EcAC) +subroutine phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,EcAC) ! Compute the correlation energy via the adiabatic connection formula @@ -9,7 +9,12 @@ subroutine ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ER integer,intent(in) :: ispin logical,intent(in) :: exchange_kernel - integer,intent(in) :: nBas,nC,nO,nV,nR,nS + 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) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XmY(nS,nS) @@ -56,11 +61,9 @@ subroutine ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ER do b=nO+1,nBas-nR jb = jb + 1 - Ap(ia,jb) = (1d0 + delta_spin)*ERI(i,b,a,j) & - - delta_Kx*ERI(i,b,j,a) + Ap(ia,jb) = (1d0 + delta_spin)*ERI(i,b,a,j) - delta_Kx*ERI(i,b,j,a) - Bp(ia,jb) = (1d0 + delta_spin)*ERI(i,j,a,b) & - - delta_Kx*ERI(i,j,b,a) + Bp(ia,jb) = (1d0 + delta_spin)*ERI(i,j,a,b) - delta_Kx*ERI(i,j,b,a) enddo enddo @@ -76,5 +79,4 @@ subroutine ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ER + trace_matrix(nS,matmul(X,matmul(Ap,transpose(X))) + matmul(Y,matmul(Ap,transpose(Y)))) & - trace_matrix(nS,Ap) -end subroutine ACFDT_correlation_energy - +end subroutine diff --git a/src/RPA/phRPA.f90 b/src/RPA/phRPA.f90 index 585a423..d4a9766 100644 --- a/src/RPA/phRPA.f90 +++ b/src/RPA/phRPA.f90 @@ -93,10 +93,10 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@RPA correlation energy (singlet) =',EcRPA(1) - write(*,'(2X,A50,F20.10)') 'Tr@RPA correlation energy (triplet) =',EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@RPA correlation energy =',EcRPA(1) + EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@RPA total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2) + 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(*,*)'-------------------------------------------------------------------------------' write(*,*) @@ -124,10 +124,10 @@ subroutine phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy (singlet) =',EcAC(1) - write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy (triplet) =',EcAC(2) - write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy =',EcAC(1) + EcAC(2) - write(*,'(2X,A50,F20.10)') 'AC@RPA total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) + 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(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/RPA/ppACFDT_correlation_energy.f90 b/src/RPA/ppACFDT_correlation_energy.f90 index 117a257..18a9185 100644 --- a/src/RPA/ppACFDT_correlation_energy.f90 +++ b/src/RPA/ppACFDT_correlation_energy.f90 @@ -8,7 +8,12 @@ subroutine ppACFDT_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOO,nVV,X1,Y ! Input variables integer,intent(in) :: ispin - integer,intent(in) :: nBas,nC,nO,nV,nR,nS + 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) :: ERI(nBas,nBas,nBas,nBas) integer,intent(in) :: nOO