diff --git a/input/methods b/input/methods index a73f4b4..2305a6a 100644 --- a/input/methods +++ b/input/methods @@ -13,7 +13,7 @@ # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW - F F T F F F + F F F F F F # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh - F F F F F F + F F F F T F # * unrestricted version available diff --git a/src/RPA/unrestricted_ACFDT.f90 b/src/GW/GW_UACFDT.f90 similarity index 85% rename from src/RPA/unrestricted_ACFDT.f90 rename to src/GW/GW_UACFDT.f90 index 80093d2..fe4b1c8 100644 --- a/src/RPA/unrestricted_ACFDT.f90 +++ b/src/GW/GW_UACFDT.f90 @@ -1,5 +1,5 @@ -subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin_flip,eta, & - nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eW,e,EcAC) +subroutine GW_UACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin_flip,eta, & + nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eW,e,EcAC) ! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem @@ -46,12 +46,12 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons double precision,allocatable :: rho_RPA(:,:,:,:) integer :: nS_aa,nS_bb,nS_sc - double precision,allocatable :: Omega_sc(:) + double precision,allocatable :: Om_sc(:) double precision,allocatable :: XpY_sc(:,:) double precision,allocatable :: XmY_sc(:,:) integer :: nS_ab,nS_ba,nS_sf - double precision,allocatable :: Omega_sf(:) + double precision,allocatable :: Om_sf(:) double precision,allocatable :: XpY_sf(:,:) double precision,allocatable :: XmY_sf(:,:) @@ -103,7 +103,7 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons ispin = 1 - allocate(Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) + allocate(Om_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) write(*,*) '------------------------' write(*,*) 'Spin-conserved manifold ' @@ -127,10 +127,10 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons end if call phULR(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),Omega_sc,XpY_sc,XmY_sc) + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),Om_sc,XpY_sc,XmY_sc) - call unrestricted_ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc, & - ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sc,XmY_sc,Ec(iAC,ispin)) + call UACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc, & + ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sc,XmY_sc,Ec(iAC,ispin)) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) @@ -145,7 +145,7 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons write(*,*) '-----------------------------------------------------------------------------------' write(*,*) - deallocate(Omega_sc,XpY_sc,XmY_sc) + deallocate(Om_sc,XpY_sc,XmY_sc) end if @@ -157,7 +157,7 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons ! Memory allocation - allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) + allocate(Om_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) write(*,*) '--------------------' write(*,*) ' Spin-flip manifold ' @@ -181,10 +181,10 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons end if call phULR(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,lambda,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),Omega_sf,XpY_sf,XmY_sf) + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),Om_sf,XpY_sf,XmY_sf) - call unrestricted_ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf, & - ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sf,XmY_sf,Ec(iAC,ispin)) + call UACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf, & + ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sf,XmY_sf,Ec(iAC,ispin)) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) @@ -199,7 +199,7 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons write(*,*) '-----------------------------------------------------------------------------------' write(*,*) - deallocate(Omega_sf,XpY_sf,XmY_sf) + deallocate(Om_sf,XpY_sf,XmY_sf) end if diff --git a/src/GW/UG0W0.f90 b/src/GW/UG0W0.f90 index 8598672..73eea6f 100644 --- a/src/GW/UG0W0.f90 +++ b/src/GW/UG0W0.f90 @@ -223,8 +223,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons end if - call unrestricted_ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,spin_conserved,spin_flip,eta, & - nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGW,EcAC) + call GW_UACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,spin_conserved,spin_flip,eta, & + nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGW,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GW/evUGW.f90 b/src/GW/evUGW.f90 index 93a0fc6..a459c18 100644 --- a/src/GW/evUGW.f90 +++ b/src/GW/evUGW.f90 @@ -271,8 +271,8 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, end if - call unrestricted_ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,spin_conserved,spin_flip, & - eta,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eGW,eGW,EcAC) + call GW_UACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,spin_conserved,spin_flip, & + eta,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eGW,eGW,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GW/qsUGW.f90 b/src/GW/qsUGW.f90 index 148f890..883949d 100644 --- a/src/GW/qsUGW.f90 +++ b/src/GW/qsUGW.f90 @@ -401,8 +401,8 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W, end if - call unrestricted_ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,spin_conserved,spin_flip, & - eta,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eGW,eGW,EcAC) + call GW_UACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,spin_conserved,spin_flip, & + eta,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eGW,eGW,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/RPA/UACFDT.f90 b/src/RPA/UACFDT.f90 new file mode 100644 index 0000000..dacfed8 --- /dev/null +++ b/src/RPA/UACFDT.f90 @@ -0,0 +1,206 @@ +subroutine UACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin_flip,eta, & + nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eW,e,EcAC) + +! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem + + implicit none + include 'parameters.h' + include 'quadrature.h' + +! Input variables + + logical,intent(in) :: doXBS + logical,intent(in) :: dRPA + logical,intent(in) :: TDA_W + logical,intent(in) :: TDA + logical,intent(in) :: BSE + logical,intent(in) :: exchange_kernel + logical,intent(in) :: spin_conserved + logical,intent(in) :: spin_flip + + double precision,intent(in) :: eta + 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) + double precision,intent(in) :: eW(nBas,nspin) + double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_bbbb(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 :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:,:) + + integer :: nS_aa,nS_bb,nS_sc + double precision,allocatable :: Om_sc(:) + double precision,allocatable :: XpY_sc(:,:) + double precision,allocatable :: XmY_sc(:,:) + + integer :: nS_ab,nS_ba,nS_sf + double precision,allocatable :: Om_sf(:) + double precision,allocatable :: XpY_sf(:,:) + double precision,allocatable :: XmY_sf(:,:) + +! Output variables + + double precision,intent(out) :: EcAC(nspin) + +! Memory allocation + + allocate(Ec(nAC,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 + + ! Memory allocation + + nS_aa = nS(1) + nS_bb = nS(2) + nS_sc = nS_aa + nS_bb + + nS_ab = (nO(1) - nC(1))*(nV(2) - nR(2)) + nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1)) + nS_sf = nS_ab + nS_ba + + allocate(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin)) + + call phULR(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,eW, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + +! Spin-conserved manifold + + if(spin_conserved) then + + ispin = 1 + + allocate(Om_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) + + write(*,*) '------------------------' + write(*,*) 'Spin-conserved manifold ' + 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 phULR(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + + end if + + call phULR(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),Om_sc,XpY_sc,XmY_sc) + + call UACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc, & + ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sc,XmY_sc,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(*,*) + + deallocate(Om_sc,XpY_sc,XmY_sc) + + end if + +! spin-flip manifold + + if(spin_flip) then + + ispin = 2 + + ! Memory allocation + + allocate(Om_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) + + write(*,*) '--------------------' + write(*,*) ' Spin-flip manifold ' + 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 phULR(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + + end if + + call phULR(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,lambda,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),Om_sf,XpY_sf,XmY_sf) + + call UACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf, & + ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sf,XmY_sf,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(*,*) + + deallocate(Om_sf,XpY_sf,XmY_sf) + + end if + +end subroutine diff --git a/src/RPA/unrestricted_ACFDT_correlation_energy.f90 b/src/RPA/UACFDT_correlation_energy.f90 similarity index 95% rename from src/RPA/unrestricted_ACFDT_correlation_energy.f90 rename to src/RPA/UACFDT_correlation_energy.f90 index 3b1e8d9..b788960 100644 --- a/src/RPA/unrestricted_ACFDT_correlation_energy.f90 +++ b/src/RPA/UACFDT_correlation_energy.f90 @@ -1,5 +1,5 @@ -subroutine unrestricted_ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt, & - ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,XmY,EcAC) +subroutine UACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt, & + ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,XmY,EcAC) ! Compute the correlation energy via the adiabatic connection formula @@ -230,4 +230,4 @@ subroutine unrestricted_ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,n + trace_matrix(nSt,matmul(X,matmul(Ap,transpose(X))) + matmul(Y,matmul(Ap,transpose(Y)))) & - trace_matrix(nSt,Ap) -end subroutine unrestricted_ACFDT_correlation_energy +end subroutine diff --git a/src/RPA/URPA.f90 b/src/RPA/URPA.f90 index 64874dd..eb8b81d 100644 --- a/src/RPA/URPA.f90 +++ b/src/RPA/URPA.f90 @@ -48,7 +48,6 @@ subroutine URPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC double precision :: rho_sc,rho_sf double precision :: EcRPA(nspin) - double precision :: EcAC(nspin) ! Hello world @@ -68,7 +67,7 @@ subroutine URPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC ! Initialization EcRPA(:) = 0d0 - EcAC(:) = 0d0 + EcRPA(:) = 0d0 ! Spin-conserved transitions @@ -143,25 +142,25 @@ subroutine URPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC write(*,*) '---------------------------------------------------------' write(*,*) - call unrestricted_ACFDT(exchange_kernel,.false.,.true.,.false.,TDA,.false.,spin_conserved,spin_flip,eta, & - nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,e,e,EcAC) + call UACFDT(exchange_kernel,.false.,.true.,.false.,TDA,.false.,spin_conserved,spin_flip,eta, & + nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,e,e,EcRPA) if(exchange_kernel) then - EcAC(1) = 0.5d0*EcAC(1) - EcAC(2) = 1.5d0*EcAC(2) + EcRPA(1) = 0.5d0*EcRPA(1) + EcRPA(2) = 1.5d0*EcRPA(2) end if write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'AC@URPA correlation energy (spin-conserved) =',EcAC(1) - write(*,'(2X,A50,F20.10)') 'AC@URPA correlation energy (spin-flip) =',EcAC(2) - write(*,'(2X,A50,F20.10)') 'AC@URPA correlation energy =',EcAC(1) + EcAC(2) - write(*,'(2X,A50,F20.10)') 'AC@URPA total energy =',ENuc + EUHF + EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@URPA correlation energy (spin-conserved) =',EcRPA(1) + write(*,'(2X,A50,F20.10)') 'AC@URPA correlation energy (spin-flip) =',EcRPA(2) + write(*,'(2X,A50,F20.10)') 'AC@URPA correlation energy =',EcRPA(1) + EcRPA(2) + write(*,'(2X,A50,F20.10)') 'AC@URPA total energy =',ENuc + EUHF + EcRPA(1) + EcRPA(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) end if -end subroutine URPA +end subroutine diff --git a/src/RPA/URPAx.f90 b/src/RPA/URPAx.f90 index 5ef5b16..eb1b081 100644 --- a/src/RPA/URPAx.f90 +++ b/src/RPA/URPAx.f90 @@ -47,8 +47,7 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n double precision,allocatable :: XmY_sf(:,:) double precision :: rho_sc,rho_sf - double precision :: EcRPAx(nspin) - double precision :: EcAC(nspin) + double precision :: EcRPA(nspin) ! Hello world @@ -68,8 +67,8 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n ! Initialization - EcRPAx(:) = 0d0 - EcAC(:) = 0d0 + EcRPA(:) = 0d0 + EcRPA(:) = 0d0 ! Spin-conserved transitions @@ -86,7 +85,7 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n allocate(Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) call phULR(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,Omega_sc,rho_sc,EcRPAx(ispin),Omega_sc,XpY_sc,XmY_sc) + ERI_aaaa,ERI_aabb,ERI_bbbb,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) call print_excitation('URPAx ',5,nS_sc,Omega_sc) call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & c,S,Omega_sc,XpY_sc,XmY_sc) @@ -110,7 +109,7 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) call phULR(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sf,1d0,e, & - ERI_aaaa,ERI_aabb,ERI_bbbb,Omega_sf,rho_sf,EcRPAx(ispin),Omega_sf,XpY_sf,XmY_sf) + ERI_aaaa,ERI_aabb,ERI_bbbb,Omega_sf,rho_sf,EcRPA(ispin),Omega_sf,XpY_sf,XmY_sf) call print_excitation('URPAx ',6,nS_sf,Omega_sf) call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & c,S,Omega_sf,XpY_sf,XmY_sf) @@ -121,21 +120,21 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n if(exchange_kernel) then - EcRPAx(1) = 0.5d0*EcRPAx(1) - EcRPAx(2) = 0.5d0*EcRPAx(2) + EcRPA(1) = 0.5d0*EcRPA(1) + EcRPA(2) = 0.5d0*EcRPA(2) else - EcRPAx(2) = 0d0 + EcRPA(2) = 0d0 end if write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@URPAx correlation energy (spin-conserved) =',EcRPAx(1) - write(*,'(2X,A50,F20.10)') 'Tr@URPAx correlation energy (spin-flip) =',EcRPAx(2) - write(*,'(2X,A50,F20.10)') 'Tr@URPAx correlation energy =',EcRPAx(1) + EcRPAx(2) - write(*,'(2X,A50,F20.10)') 'Tr@URPAx total energy =',ENuc + EUHF + EcRPAx(1) + EcRPAx(2) + write(*,'(2X,A50,F20.10)') 'Tr@URPAx correlation energy (spin-conserved) =',EcRPA(1) + write(*,'(2X,A50,F20.10)') 'Tr@URPAx correlation energy (spin-flip) =',EcRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@URPAx correlation energy =',EcRPA(1) + EcRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@URPAx total energy =',ENuc + EUHF + EcRPA(1) + EcRPA(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) @@ -148,18 +147,18 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n write(*,*) '----------------------------------------------------------' write(*,*) - call unrestricted_ACFDT(exchange_kernel,.false.,.false.,.false.,TDA,.false.,spin_conserved,spin_flip,eta, & - nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,e,e,EcAC) + call UACFDT(exchange_kernel,.false.,.false.,.false.,TDA,.false.,spin_conserved,spin_flip,eta, & + nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,e,e,EcRPA) write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'AC@URPAx correlation energy (spin-conserved) =',EcAC(1) - write(*,'(2X,A50,F20.10)') 'AC@URPAx correlation energy (spin-flip) =',EcAC(2) - write(*,'(2X,A50,F20.10)') 'AC@URPAx correlation energy =',EcAC(1) + EcAC(2) - write(*,'(2X,A50,F20.10)') 'AC@URPAx total energy =',ENuc + EUHF + EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@URPAx correlation energy (spin-conserved) =',EcRPA(1) + write(*,'(2X,A50,F20.10)') 'AC@URPAx correlation energy (spin-flip) =',EcRPA(2) + write(*,'(2X,A50,F20.10)') 'AC@URPAx correlation energy =',EcRPA(1) + EcRPA(2) + write(*,'(2X,A50,F20.10)') 'AC@URPAx total energy =',ENuc + EUHF + EcRPA(1) + EcRPA(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) end if -end subroutine URPAx +end subroutine