create_spinfree_2rdm Subroutine

public subroutine create_spinfree_2rdm(rdm, nrdms_standard, spawn, rdm_recv)

Arguments

Type IntentOptional Attributes Name
type(rdm_list_t), intent(in) :: rdm
integer, intent(in) :: nrdms_standard
type(rdm_spawn_t), intent(inout) :: spawn
type(rdm_list_t), intent(inout) :: rdm_recv

Contents

Source Code


Source Code

    subroutine create_spinfree_2rdm(rdm, nrdms_standard, spawn, rdm_recv)

        ! Take an standard (spinned) 2-RDM, stored in rdm, and output the
        ! spinfree version of it to the spawn%rdm_recv object.

        ! The input RDM has elements equal to:
        !
        ! \Gamma_{ij,kl} = < a^+_i a^+_j a_l a_k >
        !
        ! where i, j, k and l are spin orbital labels, and the output spinfree
        ! RDM has elements equal to:
        !
        ! \Gamma^{spinfree}_{pq,rs} = \sum_{x,y} < a^+_{p,x} a^+_{q,y} a_{s,y} a_{r,x} >
        !
        ! where p, q, r, s are spatial orbital labels, and x and y are spin
        ! labels (alpha or beta) which are summed over.

        ! Thus, all terms with spin signature aaaa, abab, bbbb or baba are
        ! summed together. Terms with spin signature abba or baab have their
        ! final two spin orbital labels swapped (introducing a minus sign), so
        ! that they give a contribution to the resulting spinfree RDM element.

        use rdm_data_utils, only: annihilate_rdm_list
        use SystemData, only: nbasis
        use UMatCache, only: spatial

        type(rdm_list_t), intent(in) :: rdm
        integer, intent(in) :: nrdms_standard
        type(rdm_spawn_t), intent(inout) :: spawn
        type(rdm_list_t), intent(inout) :: rdm_recv

        integer(int_rdm) :: ijkl
        integer :: ielem
        integer :: ij, kl, i, j, k, l, k_orig, l_orig ! spin orbitals
        integer :: pq, rs, p, q, r, s ! spatial orbitals
        real(dp) :: rdm_sign(rdm%sign_length)
        logical :: nearly_full, finished, all_finished

        ! If we're about to fill up the spawn list, perform a communication.
        nearly_full = .false.
        ! Have we finished adding RDM elements to the spawned list?
        finished = .false.
        rdm_recv%nelements = 0

        do ielem = 1, rdm%nelements
            ! If the spawned list is nearly full, perform a communication.
            if (nearly_full) then
                call communicate_rdm_spawn_t_wrapper(spawn, rdm_recv, finished, all_finished)
                nearly_full = .false.
            end if

            ijkl = rdm%elements(0, ielem)
            ! Obtain spin orbital labels and the RDM element.
            call calc_separate_rdm_labels(ijkl, ij, kl, i, j, k, l)
            call extract_sign_rdm(rdm%elements(:, ielem), rdm_sign)

            ! Store the original labels, before we possibly swap them.
            k_orig = k; l_orig = l;
            ! If this term is abba or baab then we can make it abab or baba by
            ! swapping the last two indices, which introduces a minus sign.
            ! It will then contribute to a spinfree 2-RDM element.
            if (.not. same_spin(i, k)) then
                l = k_orig
                k = l_orig
                rdm_sign = -rdm_sign
            end if

            ! Get the spatial orbital labels from the spin orbital ones.
            p = spatial(i); q = spatial(j);
            r = spatial(k); s = spatial(l);
            ! The 'combined' labels.
            pq = (p - 1) * nbasis + q
            rs = (r - 1) * nbasis + s

            ! If the RDM is not symmetrised then the same term will be added
            ! from both above below the diagonal, so in this case we want a
            ! factor of a half to average and not double count. However,
            ! we do not apply hermiticty to transition RDMs, as they are not
            ! hermitian. So only apply this to non-transition RDMs.
            if (pq /= rs) rdm_sign(1:nrdms_standard) = rdm_sign(1:nrdms_standard) * 0.5_dp

            ! Due to the fact that RDM elements are only stored with p < q and
            ! r < s, the following terms are only stored with baba spin, never
            ! with abab. Double this term to make up for it.
            if (p == q .and. r == s) rdm_sign = 2.0_dp * rdm_sign

            ! Add all spinfree 2-RDM elements corresponding to these labels.
            call add_rdm_elements(p, q, r, s, nrdms_standard, rdm_sign, spawn, nearly_full)

            ! If this is an aaaa or bbbb term then *minus* this RDM element will
            ! be equal to the equivalent RDM element with the last two labels
            ! swapped. So, add this contribution into that RDM element. We
            ! don't have to do this, but doing so applies some extra averaging.
            ! Want to apply all the averaging possible over equivalent elements.
            if (same_spin(i, j)) then
                ! Re-extract sign in case it has been modified.
                call extract_sign_rdm(rdm%elements(:, ielem), rdm_sign)

                ! Swap the spatial labels.
                r = spatial(l_orig); s = spatial(k_orig);
                rs = (r - 1) * nbasis + s

                ! Half the sign of non-transition RDMs, where we apply
                ! hermiticity to average above and below the diagonal, to
                ! prevent double counting.
                if (pq /= rs) rdm_sign(1:nrdms_standard) = rdm_sign(1:nrdms_standard) * 0.5_dp
                rdm_sign = -rdm_sign

                call add_rdm_elements(p, q, r, s, nrdms_standard, rdm_sign, spawn, nearly_full)
            end if

        end do

        finished = .true.
        ! Keep performing communications until all RDM spawnings on every
        ! processor have been communicated.
        do
            call communicate_rdm_spawn_t_wrapper(spawn, rdm_recv, finished, all_finished)
            if (all_finished) exit
        end do
        call annihilate_rdm_list(rdm_recv)

    contains

        subroutine add_rdm_elements(p, q, r, s, nrdms_standard, rdm_sign, spawn, nearly_full)

            ! Add in the single contribution rdm_sign to the following elements
            ! of the spinfree 2-RDM:
            !
            ! \Gamma^{spinfree}_{pq,rs} = \sum_{x,y} < a^+_{p,x} a^+_{q,y} a_{s,y} a_{r,x} >
            ! \Gamma^{spinfree}_{qp,sr} = \sum_{x,y} < a^+_{q,x} a^+_{p,y} a_{r,y} a_{s,x} >
            !
            ! And for non-transition RDMs, which are hermitian, also add in
            ! the following elements:
            !
            ! \Gamma^{spinfree}_{rs,pq} = \sum_{x,y} < a^+_{r,x} a^+_{s,y} a_{q,y} a_{p,x} >
            ! \Gamma^{spinfree}_{sr,qp} = \sum_{x,y} < a^+_{s,x} a^+_{r,y} a_{p,y} a_{q,x} >
            !
            ! where x and y are spin labels which are summed over in the final
            ! result.
            !
            ! For a *REAL* spinfree 2-RDM, all of these elements are rigorously
            ! equal, so it is appropriate that we add all contributions in
            ! together like this (the last two aren't equal to the first two for
            ! transition RDMs, so don't add these for transition RDMs).
            !
            ! The if-statements in here prevent adding to the same RDM element
            ! twice.

            integer, intent(in) :: p, q, r, s ! spatial orbitals
            integer, intent(in) :: nrdms_standard
            real(dp), intent(in) :: rdm_sign(:)
            type(rdm_spawn_t), intent(inout) :: spawn
            logical, intent(inout) :: nearly_full

            real(dp) :: rdm_sign_temp(size(rdm_sign))

            rdm_sign_temp = 0.0_dp

            ! RDM element \Gamma_{pq,rs}.
            call add_to_rdm_spawn_t(spawn, p, q, r, s, rdm_sign, .true., nearly_full)

            ! RDM element \Gamma_{qp,sr}.
            if (.not. (p == q .and. r == s)) then
                call add_to_rdm_spawn_t(spawn, q, p, s, r, rdm_sign, .true., nearly_full)
            end if

            if (pq /= rs) then
                ! We only want to apply hermiticity symmetry to non-transition
                ! RDMs, since transition RDMs are not hermitian.
                rdm_sign_temp(1:nrdms_standard) = rdm_sign(1:nrdms_standard)

                ! RDM element \Gamma_{rs,pq}.
                call add_to_rdm_spawn_t(spawn, r, s, p, q, rdm_sign_temp, .true., nearly_full)

                ! RDM element \Gamma_{sr,qp}.
                if (.not. (p == q .and. r == s)) then
                    call add_to_rdm_spawn_t(spawn, s, r, q, p, rdm_sign_temp, .true., nearly_full)
                end if
            end if

        end subroutine add_rdm_elements

    end subroutine create_spinfree_2rdm