Bug report: Is the orbital symmetrization in `enforce_symmetry` correct?

[xfzhang-phys]

Dear all,

I may be misunderstanding the implementation, but it appears that the _symmetrize_orbital function in symmetries.py may pair upper- and lower-triangular orbital indices incorrectly.

As a simple check, I tested the _split_orbital_triangle function with the following script:

import numpy as np


def _split_orbital_triangle(norb_):
    norb = norb_
    upper_triangle_slice = np.triu_indices(norb, k=1)
    lower_triangle_slice = np.tril_indices(norb, k=-1)
    combined_slice = upper_triangle_slice + lower_triangle_slice

    for upper_1, upper_2, lower_1, lower_2 in zip(*combined_slice):
        upper_triangle = [(upper_1, upper_2)]
        lower_triangle = [(lower_1, lower_2)]
        yield upper_triangle, lower_triangle


for upper_triangle, lower_triangle in _split_orbital_triangle(5):
    print(upper_triangle)
    print(lower_triangle)
    print()

The output begins as follows:

[(np.int64(0), np.int64(1))]
[(np.int64(1), np.int64(0))]

[(np.int64(0), np.int64(2))]
[(np.int64(2), np.int64(0))]

[(np.int64(0), np.int64(3))]
[(np.int64(2), np.int64(1))]

[(np.int64(0), np.int64(4))]
[(np.int64(3), np.int64(0))]

[(np.int64(1), np.int64(2))]
[(np.int64(3), np.int64(1))]

[(np.int64(1), np.int64(3))]
[(np.int64(3), np.int64(2))]

[(np.int64(1), np.int64(4))]
[(np.int64(4), np.int64(0))]

[(np.int64(2), np.int64(3))]
[(np.int64(4), np.int64(1))]

[(np.int64(2), np.int64(4))]
...

For example, the upper-triangular index (0, 3) is paired with the lower-triangular index (2, 1), whereas I would have expected it to be paired with (3, 0).

Could you please confirm whether this behavior is intentional? If so, I would appreciate some clarification on how these index pairs are intended to be used in the orbital symmetrization. Otherwise, this may lead to incorrect symmetrization of the gap function.

best,
xf