Band-decomposed (partial) charge densities using PW basis set appear unnormalized

[AsTonyshment]

Describe the bug

While performing band-decomposed (partial) charge density calculations using the PW basis set, I've observed that the charge density integrals are not normalized and are incrementally increasing across bands (using a simple script written by Python attached below). This trend suggests that there might be an issue with array reinitialization in the subsequent computational steps:

$ python3 grid_integral.py
电荷密度积分：band1.cube 1.9999983547907862
电荷密度积分：band2.cube 3.999996709581772
电荷密度积分：band3.cube 5.999995064372504
电荷密度积分：band4.cube 7.999993419163174
电荷密度积分：band5.cube 9.999991773953775
电荷密度积分：band6.cube 11.999990128744749
电荷密度积分：band7.cube 13.999988483535311
电荷密度积分：band8.cube 15.999986838326754

Expected behavior

• For LCAO (one might have to add additional INPUT parameters such as calculation=get_pchg), the integral values should remain consistent and close to 2. Here's one of my results:

电荷密度积分：BAND1_SPIN0_CHG.cube 1.9999983010347646
电荷密度积分：BAND2_SPIN0_CHG.cube 1.9999983636754213
电荷密度积分：BAND3_SPIN0_CHG.cube 1.9999982713382998
电荷密度积分：BAND4_SPIN0_CHG.cube 1.9999982579638917

• For PW, note the increasing trend in the integral values across bands, which is unexpected and indicates a possible bug.

To Reproduce

Environment Setup

• Input Files: Provided in the attachments, including INPUT, STRU, and a Python script grid_integral.py.
• Script Functionality: grid_integral.py is designed to compute the total space integral of charge density from .cube files. A result close to 2 is expected when nspin=1.

Steps to Reproduce

1. Run ABACUS using the provided INPUT and STRU files.
2. Execute the grid_integral.py script to calculate the charge density integrals.

Input files.zip

Environment

No response

Additional Context

I suspect that this bug was introduced in PR #3501. I plan to fix it along with PR #3934.

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).

[WHUweiqingzhou]

@AsTonyshment
As far as I know, only LCAO calculations use renormalization of rho. There is no renormalization of rho in PW calculations.

[AsTonyshment]

@AsTonyshment As far as I know, only LCAO calculations use renormalization of rho. There is no renormalization of rho in PW calculations.

What I try to express is that, in the code source/module_esolver/esolver_ks_pw.cpp, the initialization of the band-decomposed charge density array appears to be placed before the main loop over the band index, starting at line 1030:

    for(int i = 0; i < this->pw_rho->nxyz; i++)
    {
        rho_band[i] = 0.0;
    }

    for(int i = 0; i < INPUT.band_print_num; i++)
    {
        int ib = INPUT.bands_to_print[i];
        // further processing

This setup results in the charge density of the next band being directly added to the cumulative result of the previous loop.