Bug: Incorrect calculation of partial charge for multi-k-points under both PW and LCAO when symmetry is turned on

[AsTonyshment]

Describe the bug

I am currently writing a tutorial on how to calculate and visualize band-decomposed (partial) charge density. However, when I attempt to calculate and draw the .cube charge file, the results from PW and LCAO basis set appear to be inconsistent. Attached are some benchmark results using 2D graphene with a 2×2×2 k-points sampling.

BAND1 (VASP)

BAND1 (ABACUS, LCAO)

BAND1 (ABACUS, PW)

It appears that under the PW basis set, only the charge of certain specific k-points is taken into account. The observed pattern matches the result of the k-point (0.5, 0.5, 0.5) (in terms of k-points sampling using DIRECT coordinates in output file kpoints):

BAND1 K(0.5, 0.5, 0.5) (VASP)

BAND1 K(0.5, 0.5, 0.5) (ABACUS, LCAO)

Task list for Issue attackers (only for developers)

• Verify the issue is not a duplicate.
• Describe the bug.
• Steps to reproduce.
• Expected behavior.
• Error message.
• Environment details.
• Additional context.
• Assign a priority level (low, medium, high, urgent).
• Assign the issue to a team member.
• Label the issue with relevant tags.
• Identify possible related issues.
• Create a unit test or automated test to reproduce the bug (if applicable).
• Fix the bug.
• Test the fix.
• Update documentation (if necessary).
• Close the issue and inform the reporter (if applicable).

[AsTonyshment]

Furthermore, the current PW basis set does not support outputting partial charge densities separately for different spins and different k-points. This functionality should also be implemented. This will be addressed in a separate Feature Request issue.

[AsTonyshment]

New updates: After further testing, I found that the incorrect result only occurs when symmetry = 1. When I set symmetry = -1, the results are still correct:

BAND1 (ABACUS, PW, turning symmetry off)

This is very wired... Further investigation is needed.

[AsTonyshment]

Now, it turns out that all results of ABACUS with symmetry = 1 are wrong, including LCAO... I suspect this issue is related to the symmetry analysis in ABACUS.

[AsTonyshment]

I have a rough idea of what the problem is. When symmetry = 1, the symmetry reduction is based on lattice symmetry. Therefore, the charge density needs to undergo corresponding symmetry operations to return to its true state in real space. However, ABACUS currently does not perform such operations on the band-decomposed charge density (both PW and LCAO).

For example, in the diagram below, K2, K3, and K4 are all reduced to the M point. However, their charge densities in real space differ by a lattice symmetry operation. The current algorithm directly sums them according to the k-point weights without considering these symmetry operations.