add_rdm_1_to_rdm_2 Subroutine

    subroutine add_rdm_1_to_rdm_2(rdm_1, rdm_2, scale_factor)

        ! Take the RDM elements in the rdm_1 object, and add them to the rdm_2
        ! object. The hash table for rdm_2 will also be updated. This literally
        ! performs the numerical addition of the two RDM objects.
        use SystemData, only: nel, nBasis
        use hash, only: hash_table_lookup, add_hash_table_entry
        use Parallel_neci, only: MPISumAll

        type(rdm_list_t), intent(in) :: rdm_1
        type(rdm_list_t), intent(inout) :: rdm_2
        real(dp), intent(in), optional :: scale_factor

        integer :: ielem
        integer(int_rdm) :: ijkl
        integer :: ij, kl, i, j, k, l ! spin or spatial orbitals
        integer :: ind, hash_val
        real(dp) :: real_sign_old(rdm_2%sign_length), real_sign_new(rdm_2%sign_length)
        real(dp) :: spawn_sign(rdm_2%sign_length)
        real(dp) :: internal_scale_factor(rdm_1%sign_length), rdm_trace(rdm_1%sign_length)
        logical :: tSuccess
        character(*), parameter :: this_routine = 'add_rdm_1_to_rdm_2'

        if (present(scale_factor)) then
            if (tGUGA) then
                call stop_all(this_routine, "check scale factor and trace for GUGA!")
            end if
            ! normalize and rescale the rdm_1 if requested here
            call calc_rdm_trace(rdm_1, rdm_trace)
            call MPISumAll(rdm_trace, internal_scale_factor)
            internal_scale_factor = scale_factor * (nel * (nel - 1)) / (2 * internal_scale_factor)
        else
            internal_scale_factor = 1.0_dp
        end if

        if (rdm_1%sign_length /= rdm_2%sign_length) &
            call stop_all(this_routine, "nrdms mismatch")

        do ielem = 1, rdm_1%nelements
            ! Decode the compressed RDM labels.
            ijkl = rdm_1%elements(0, ielem)
            if (tGUGA) then
                call extract_2_rdm_ind(ijkl, i, j, k, l)
            else
                call calc_separate_rdm_labels(ijkl, ij, kl, i, j, k, l)
            end if

            ! Extract the spawned sign.
            call extract_sign_rdm(rdm_1%elements(:, ielem), spawn_sign)

            ! Search to see if this RDM element is already in the RDM 2.
            ! If it, tSuccess will be true and ind will hold the position of the
            ! element in rdm.
            ASSERT(all([i, j, k, l] > 0) .and. all([i, j, k, l] <= nBasis))

            call hash_table_lookup((/i, j, k, l/), (/ijkl/), 0, rdm_2%hash_table, rdm_2%elements, ind, hash_val, tSuccess)

            if (tSuccess) then
                ! Extract the existing sign.
                call extract_sign_rdm(rdm_2%elements(:, ind), real_sign_old)

                ! Update the total sign.
                real_sign_new = real_sign_old + spawn_sign * internal_scale_factor
                ! Encode the new sign.
                call encode_sign_rdm(rdm_2%elements(:, ind), real_sign_new)
            else
                ! If we don't have enough memory in rdm_2, try increasing its
                ! size to be 1.5 times bigger.
                if (rdm_2%nelements + 1 > rdm_2%max_nelements) then
                    call try_rdm_list_realloc(rdm_2, int(1.5 * rdm_2%max_nelements), .false.)
                end if

                ! Update the rdm array, and its hash table, and the number of
                ! RDM elements.
                rdm_2%nelements = rdm_2%nelements + 1
                rdm_2%elements(0, rdm_2%nelements) = ijkl
                call encode_sign_rdm(rdm_2%elements(:, rdm_2%nelements), spawn_sign)
                call add_hash_table_entry(rdm_2%hash_table, rdm_2%nelements, hash_val)
            end if
        end do

    end subroutine add_rdm_1_to_rdm_2