calc_hermitian_errors Subroutine

    subroutine calc_hermitian_errors(rdm, rdm_recv, spawn, rdm_norm, max_error_herm_all, sum_error_herm_all)

        ! Calculate the hermiticity errors, i.e.
        !
        ! \Gamma_{ij,kl} - \Gamma_{kl,ij}^*,
        !
        ! which should be equal to zero for an exact 2-RDM.

        ! Specifically we return the largest such error, and the sum of all
        ! such errors, to max_error_herm_all and sum_error_herm_all.

        ! Use the rdm_recv and spawn objects as space for performing a
        ! communication of a 2-RDM. The hermiticity error 2-RDM which is
        ! communicated is defined as
        !
        ! \Gamma^{error}_{ij,kl} = \Gamma_{ij,kl} - \Gamma_{kl,ij}^*,
        !
        ! for ij < kl.

        use hash, only: clear_hash_table
        use Parallel_neci, only: MPIAllReduce, nProcessors
        use MPI_wrapper, only: MPI_SUM, MPI_MAX
        use rdm_data_utils, only: add_to_rdm_spawn_t, communicate_rdm_spawn_t_wrapper, annihilate_rdm_list

        type(rdm_list_t), intent(in) :: rdm
        type(rdm_list_t), intent(inout) :: rdm_recv
        type(rdm_spawn_t), intent(inout) :: spawn
        real(dp), intent(in) :: rdm_norm(rdm%sign_length)
        real(dp), intent(out) :: max_error_herm_all(rdm%sign_length)
        real(dp), intent(out) :: sum_error_herm_all(rdm%sign_length)
        character(*), parameter :: this_routine = "calc_hermitian_errors"

        integer(int_rdm) :: ijkl
        integer :: ielem
        integer :: ij, kl, i, j, k, l ! spin orbitals
        integer(int_rdm) :: ij_, kl_
        real(dp) :: rdm_sign(rdm%sign_length)
        real(dp) :: max_error_herm(rdm%sign_length), sum_error_herm(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

        ! Clear the spawn object before we use it.
        spawn%free_slots = spawn%init_free_slots(0:nProcessors - 1)
        call clear_hash_table(spawn%rdm_send%hash_table)

        ! Loop over all RDM elements.
        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 RDM element sign
            call extract_sign_rdm(rdm%elements(:, ielem), rdm_sign)

            if (tGUGA) then
                call extract_2_rdm_ind(ijkl, i, j, k, l, ij_, kl_)
                call add_to_rdm_spawn_t(spawn, i, j, k, l, rdm_sign, .true., nearly_full)
                call add_to_rdm_spawn_t(spawn, k, l, i, j, -rdm_sign, .true., nearly_full)
            else
                ! Obtain spin orbital labels
                call calc_separate_rdm_labels(ijkl, ij, kl, i, j, k, l)
                ! If in the lower half of the RDM, reflect to the upper half and
                ! include with a minus sign.
                if (ij > kl) then
                    call add_to_rdm_spawn_t(spawn, k, l, i, j, -rdm_sign, .false., nearly_full)
                else if (ij < kl) then
                    call add_to_rdm_spawn_t(spawn, i, j, k, l, rdm_sign, .false., nearly_full)
                end if
            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)

        max_error_herm = 0.0_dp
        sum_error_herm = 0.0_dp

        ! Find the largest error and sum of errrors on this processor.
        do ielem = 1, rdm_recv%nelements
            call extract_sign_rdm(rdm_recv%elements(:, ielem), rdm_sign)
            rdm_sign = abs(rdm_sign)
            max_error_herm = max(max_error_herm, rdm_sign)
            sum_error_herm = sum_error_herm + rdm_sign
        end do

        ! The input 2-RDM wasn't normalised, so need to normalise these.
        max_error_herm = max_error_herm / rdm_norm
        sum_error_herm = sum_error_herm / rdm_norm

        ! Find the largest error and sum of errors across all processors.
        call MPIAllReduce(max_error_herm, MPI_MAX, max_error_herm_all)
        call MPIAllReduce(sum_error_herm, MPI_SUM, sum_error_herm_all)

    end subroutine calc_hermitian_errors