From ce10bbaf5631b5397c1ca98bcc056aea1978a27f Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 23 Sep 2020 09:46:44 +0200 Subject: [PATCH] URPA --- input/methods | 4 +- src/QuAcK/QuAcK.f90 | 24 ++++- src/QuAcK/RPAx.f90 | 13 ++- src/QuAcK/URPAx.f90 | 152 ++++++++++++++++++++++++++++++++ src/QuAcK/UdRPA.f90 | 152 ++++++++++++++++++++++++++++++++ src/QuAcK/{RPA.f90 => dRPA.f90} | 4 +- 6 files changed, 334 insertions(+), 15 deletions(-) create mode 100644 src/QuAcK/URPAx.f90 create mode 100644 src/QuAcK/UdRPA.f90 rename src/QuAcK/{RPA.f90 => dRPA.f90} (97%) diff --git a/input/methods b/input/methods index ba17592..022df90 100644 --- a/input/methods +++ b/input/methods @@ -9,11 +9,11 @@ # CIS CID CISD F F F # RPA RPAx ppRPA - F F F + T F F # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0 evGW qsGW - T F F + F F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index d65bfcc..3bc4577 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -634,8 +634,16 @@ program QuAcK if(doRPA) then call cpu_time(start_RPA) - call RPA(doACFDT,exchange_kernel,singlet,triplet,0d0, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + if(unrestricted) then + + call UdRPA(doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & + ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,eHF) + + else + + call dRPA(doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + + end if call cpu_time(end_RPA) t_RPA = end_RPA - start_RPA @@ -651,8 +659,16 @@ program QuAcK if(doRPAx) then call cpu_time(start_RPAx) - call RPAx(doACFDT,exchange_kernel,singlet,triplet,0d0, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + if(unrestricted) then + + call URPAx(doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & + ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,eHF) + + else + + call RPAx(doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + + end if call cpu_time(end_RPAx) t_RPAx = end_RPAx - start_RPAx diff --git a/src/QuAcK/RPAx.f90 b/src/QuAcK/RPAx.f90 index 8afd56d..06c298a 100644 --- a/src/QuAcK/RPAx.f90 +++ b/src/QuAcK/RPAx.f90 @@ -1,5 +1,4 @@ -subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e) +subroutine RPAx(doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e) ! Perform random phase approximation calculation with exchange (aka TDHF) @@ -11,9 +10,9 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & logical,intent(in) :: doACFDT logical,intent(in) :: exchange_kernel - logical,intent(in) :: singlet_manifold + logical,intent(in) :: singlet double precision,intent(in) :: eta - logical,intent(in) :: triplet_manifold + logical,intent(in) :: triplet integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO @@ -55,7 +54,7 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & ! Singlet manifold - if(singlet_manifold) then + if(singlet) then ispin = 1 @@ -69,7 +68,7 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & ! Triplet manifold - if(triplet_manifold) then + if(triplet) then ispin = 2 @@ -105,7 +104,7 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & write(*,*) '-------------------------------------------------------' write(*,*) - call ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,singlet_manifold,triplet_manifold,eta, & + call ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,singlet,triplet,eta, & nBas,nC,nO,nV,nR,nS,ERI,e,e,Omega,XpY,XmY,rho,EcAC) if(exchange_kernel) then diff --git a/src/QuAcK/URPAx.f90 b/src/QuAcK/URPAx.f90 new file mode 100644 index 0000000..4ae48cc --- /dev/null +++ b/src/QuAcK/URPAx.f90 @@ -0,0 +1,152 @@ +subroutine URPAx(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,e) + +! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism + + implicit none + include 'parameters.h' + include 'quadrature.h' + +! Input variables + + double precision,intent(in) :: eta + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: spin_conserved + logical,intent(in) :: spin_flip + 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) :: ENuc + double precision,intent(in) :: EUHF + 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) + double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: ispin + + integer :: nS_aa,nS_bb,nS_sc + double precision,allocatable :: Omega_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 :: XpY_sf(:,:) + double precision,allocatable :: XmY_sf(:,:) + + double precision :: rho_sc,rho_sf + double precision :: EcRPAx(nspin) + double precision :: EcAC(nspin) + +! Hello world + + write(*,*) + write(*,*)'*********************************************************************' + write(*,*)'| Unrestricted random phase approximation calculation with exchange |' + write(*,*)'*********************************************************************' + write(*,*) + +! Initialization + + EcRPAx(:) = 0d0 + EcAC(:) = 0d0 + +! Spin-conserved transitions + + if(spin_conserved) then + + ispin = 1 + + ! Memory allocation + + nS_aa = nS(1) + nS_bb = nS(2) + nS_sc = nS_aa + nS_bb + + allocate(Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) + + call unrestricted_linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_sc,EcRPAx(ispin),Omega_sc,XpY_sc,XmY_sc) + call print_excitation('URPAx ',5,nS_sc,Omega_sc) +! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + + + endif + +! Spin-flip transitions + + if(spin_flip) then + + ispin = 2 + + ! Memory allocation + + 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(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) + + call unrestricted_linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,1d0,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_sf,EcRPAx(ispin),Omega_sf,XpY_sf,XmY_sf) + call print_excitation('URPAx ',6,nS_sf,Omega_sf) +! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + + + endif + +! if(exchange_kernel) then + +! EcRPAx(1) = 0.5d0*EcRPAx(1) +! EcRPAx(2) = 1.5d0*EcRPAx(2) + +! 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(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +! Compute the correlation energy via the adiabatic connection + +! if(doACFDT) then + +! write(*,*) '-------------------------------------------------------' +! write(*,*) 'Adiabatic connection version of RPAx correlation energy' +! write(*,*) '-------------------------------------------------------' +! write(*,*) + +! call ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,singlet,triplet,eta, & +! nBas,nC,nO,nV,nR,nS,ERI,e,e,Omega,XpY,XmY,rho,EcAC) + +! if(exchange_kernel) then + +! EcAC(1) = 0.5d0*EcAC(1) +! EcAC(2) = 1.5d0*EcAC(2) + +! end if + +! 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 + EUHF + EcAC(1) + EcAC(2) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! end if + +end subroutine URPAx diff --git a/src/QuAcK/UdRPA.f90 b/src/QuAcK/UdRPA.f90 new file mode 100644 index 0000000..da3d019 --- /dev/null +++ b/src/QuAcK/UdRPA.f90 @@ -0,0 +1,152 @@ +subroutine UdRPA(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,e) + +! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism + + implicit none + include 'parameters.h' + include 'quadrature.h' + +! Input variables + + double precision,intent(in) :: eta + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: spin_conserved + logical,intent(in) :: spin_flip + 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) :: ENuc + double precision,intent(in) :: EUHF + 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) + double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: ispin + + integer :: nS_aa,nS_bb,nS_sc + double precision,allocatable :: Omega_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 :: XpY_sf(:,:) + double precision,allocatable :: XmY_sf(:,:) + + double precision :: rho_sc,rho_sf + double precision :: EcRPA(nspin) + double precision :: EcAC(nspin) + +! Hello world + + write(*,*) + write(*,*)'**************************************************************' + write(*,*)'| Unrestricted direct random phase approximation calculation |' + write(*,*)'**************************************************************' + write(*,*) + +! Initialization + + EcRPA(:) = 0d0 + EcAC(:) = 0d0 + +! Spin-conserved transitions + + if(spin_conserved) then + + ispin = 1 + + ! Memory allocation + + nS_aa = nS(1) + nS_bb = nS(2) + nS_sc = nS_aa + nS_bb + + allocate(Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) + + call unrestricted_linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) + call print_excitation('URPA ',5,nS_sc,Omega_sc) +! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + + + endif + +! Spin-flip transitions + + if(spin_flip) then + + ispin = 2 + + ! Memory allocation + + 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(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) + + call unrestricted_linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,1d0,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_sf,EcRPA(ispin),Omega_sf,XpY_sf,XmY_sf) + call print_excitation('URPA ',6,nS_sf,Omega_sf) +! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + + + endif + +! if(exchange_kernel) then + +! EcRPA(1) = 0.5d0*EcRPA(1) +! EcRPA(2) = 1.5d0*EcRPA(2) + +! end if + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@URPA correlation energy (spin-conserved) =',EcRPA(1) + write(*,'(2X,A50,F20.10)') 'Tr@URPA correlation energy (spin-flip) =',EcRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@URPA correlation energy =',EcRPA(1) + EcRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@URPA total energy =',ENuc + EUHF + EcRPA(1) + EcRPA(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +! Compute the correlation energy via the adiabatic connection + +! if(doACFDT) then + +! write(*,*) '-------------------------------------------------------' +! write(*,*) 'Adiabatic connection version of RPA correlation energy' +! write(*,*) '-------------------------------------------------------' +! write(*,*) + +! call ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,singlet,triplet,eta, & +! nBas,nC,nO,nV,nR,nS,ERI,e,e,Omega,XpY,XmY,rho,EcAC) + +! if(exchange_kernel) then + +! EcAC(1) = 0.5d0*EcAC(1) +! EcAC(2) = 1.5d0*EcAC(2) + +! end if + +! 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 + EUHF + EcAC(1) + EcAC(2) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! end if + +end subroutine UdRPA diff --git a/src/QuAcK/RPA.f90 b/src/QuAcK/dRPA.f90 similarity index 97% rename from src/QuAcK/RPA.f90 rename to src/QuAcK/dRPA.f90 index 32d9587..b807497 100644 --- a/src/QuAcK/RPA.f90 +++ b/src/QuAcK/dRPA.f90 @@ -1,4 +1,4 @@ -subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & +subroutine dRPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e) ! Perform a direct random phase approximation calculation @@ -125,4 +125,4 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & end if -end subroutine RPA +end subroutine dRPA