Negative beta rho in DFT calculations

[BariumOxide13716]

Describe the bug

This issue is related to #6086 .

In module_elecstate/module_charge/charge_mixing_rho.cpp, the mixed alpah and beta rhos are computed as:
chr->rho[0][ir] = 0.5 * (rho_mag[ir] + rho_mag[ir+nrxx]);
chr->rho[1][ir] = 0.5 * (rho_mag[ir] - rho_mag[ir+nrxx]);
,where rho_mag[ir] is the mixed spin-summed charge, and rho_mag[ir+nrxx] is the mixed spin-diifference charge.

This leads to an error in the first SCF calculation that chr->rho[1][ir] maybe negative, when, according to the manual, the default mixing coefficients for spin-difference charge mixing_beta_mag is 4 times the value of the spin-summed charge mixing_beta.

This further leads to an undefined behavior in DFT calculations, which can be simplied as
(chr->rho[0][ir] * sgn[0][ir] + chr->rho[1][ir] * sgn[1][ir]) * e_xc
(cf module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp,
in the line etxc += ModuleBase::e2 * exc[ir] * rho[ir*nspin+is] * sgn[ir*nspin+is])

,where sgn is zero if the corresponding chr->rho is less than a threshold.

In a scenerio when chr->rho[1][ir] is negative, although sgn[1][ir] is zero and excludes the negative beta charge to be multiplied in the calculation, e_xc is still calculated with a negative chr->rho[1][ir], and multipled with the alpha charge if sgn[0][ir] is 1. This further leads to some disagreements between two calculations.

To identify the negative beta charge, I add a printout before
chr->rho[0][ir] = 0.5 * (rho_mag[ir] + rho_mag[ir+nrxx]);
chr->rho[1][ir] = 0.5 * (rho_mag[ir] - rho_mag[ir+nrxx]);
namely if (rho_mag[ir] < rho_mag[ir+nrxx]) std::cout << "m greater than rho at " << ir << std::endl;

In a calculation using R2SCAN together with libxc, when I set mixing_beta to 0.4 without setting mixing_beta_mag, I obtain different intermediate SCF energies as follows, with the left running on two processors and right on 1:

Note that at some grid points, m is greater than rho. And also note that the SCF energies are different in two calculations.

However, when I set mixing_beta_mag the same as mixing_beta, (actually there are a few values between 0.4 and 1.6 that give the same results as shown below) the negative beta charge is gone and the results are much more consistent between two calculations.

Note that the information of "m is greater than rho" is gone.

So should we add some warnings to avoid producing negative beta charges on grid points?

Expected behavior

No response

To Reproduce

To reproduce, set different values for mixing_beta_mag in the INPUT file. test_218.zip

Environment

No response

Additional Context

No response

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