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