contract_molcas_2_rdm_index Function

    pure function contract_molcas_2_rdm_index(p, q, r, s) result(pqrs)
        ! function using the molcas rdm index convention
        integer, intent(in) :: p, q, r, s
        integer :: pqrs

        integer :: pq, rs, ij, kl

        pq = contract_molcas_1_rdm_index(p, q)
        rs = contract_molcas_1_rdm_index(r, s)

        ij = max(pq, rs)
        kl = min(pq, rs)

        pqrs = kl + ij * (ij - 1) / 2

    end function contract_molcas_2_rdm_index