A significant difference in energy is observed during the relaxation process of H2 by different versions

[ieiue]

Describe the bug

When running a relax calculation of H2, there are significant different energies between different versions(INPUT, STRU, KPT are the same, keep sigma=0.015, gauss method).

1. v3.4.0: ETOT=-27.19665397014 eV
   
2. v3.5.4:ETOT=-31.741 eV
   

Expected behavior

Why the energy differences between the different versions are so obvious？More than 4 eV.

To Reproduce

H2_pw.zip

Environment

v3.4.0 & v3.5.4

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

[ieiue]

3. v3.6.4: ETOT=-31.741 eV
   

[Qianruipku]

The results are changed after PR #3226 commit fd7654

The difference is from parameter "mixing_beta_mag"

[Qianruipku]

Now, we can restate this issue:
For the latest version of the abacus, the results are significantly different when the mixing_beta_mag parameter is set to 0.4 and 1.6.
Example:
INPUT:

INPUT_PARAMETERS
suffix 		H2_10A_te_3.4
pseudo_dir 		./
ecutwfc 		100 
scf_thr 		1e-8
basis_type 		pw
calculation 	    scf
mixing_type     	broyden 
mixing_beta     	 0.4
mixing_beta_mag     *0.4 or 1.6*
mixing_gg0		0
scf_nmax        100
cal_stress      	True
nspin          	2
smearing_method	gauss
smearing_sigma		0.015
symmetry		1

For 1.6:

final etot is -31.739723622 eV
 EFERMI = -3.8375104048 eV

 STATE ENERGY(eV) AND OCCUPATIONS    NSPIN == 2
SPIN UP : 
 1/1 kpoint (Cartesian) = 0.0000 0.0000 0.0000 (14226 pws)
       1       -10.3568        1.00000
       2      -0.216514        0.00000
       3       0.389821        0.00000
       4       0.917017        0.00000
       5       0.917017        0.00000
       6        1.30126        0.00000
       7        1.32317        0.00000
       8        1.35560        0.00000
       9        1.89282        0.00000
      10        2.26524        0.00000
      11        2.26524        0.00000
      12        2.65818        0.00000
      13        2.65821        0.00000

SPIN DOWN : 
 1/1 kpoint (Cartesian) = 0.0000 0.0000 0.0000 (14226 pws)
       1       -10.3556        1.00000
       2      -0.192317        0.00000
       3       0.400011        0.00000
       4       0.930567        0.00000
       5       0.930567        0.00000
       6        1.32976        0.00000
       7        1.35412        0.00000
       8        1.37641        0.00000
       9        1.91771        0.00000
      10        2.28104        0.00000
      11        2.28104        0.00000
      12        2.68180        0.00000
      13        2.68179        0.00000

For 0.4:

final etot is -21.562935575 eV
 EFERMI = -1.750123649 eV

 STATE ENERGY(eV) AND OCCUPATIONS    NSPIN == 2
SPIN UP : 
 1/1 kpoint (Cartesian) = 0.0000 0.0000 0.0000 (14226 pws)
       1       -17.0406        1.00000
       2       -2.40799       0.999997
       3       -1.09225    2.57331e-06
       4      -0.176180        0.00000
       5      -0.176180        0.00000
       6       0.540090        0.00000
       7       0.788138        0.00000
       8       0.807214        0.00000
       9       0.876941        0.00000
      10        1.54734        0.00000
      11        1.54735        0.00000
      12        1.76000        0.00000
      13        1.76000        0.00000

SPIN DOWN : 
 1/1 kpoint (Cartesian) = 0.0000 0.0000 0.0000 (14226 pws)
       1       0.100745        0.00000
       2        1.48992        0.00000
       3        1.50678        0.00000
       4        1.57374        0.00000
       5        1.57374        0.00000
       6        1.66962        0.00000
       7        1.97928        0.00000
       8        2.91686        0.00000
       9        2.98261        0.00000
      10        2.98261        0.00000
      11        3.05105        0.00000
      12        3.07624        0.00000
      13        3.15614        0.00000

[Qianruipku]

@WHUweiqingzhou Could you have a look?

[Qianruipku]

When mixing_beta_mag is equal to mixing_beta, the result will have a non-zeror total magnetization; if they are different, the result will tend towards a total magnetization of zero. Maybe the latter is right, since H2 should have a zero total magnetization and the latter has lower energy. The former is only a local miminum.

For QE, mixing_beta_mag is equal to mixing_beta, so they get a "total magnetization = 2.00 Bohr mag/cell" result.

[kirk0830]

When mixing_beta_mag is equal to mixing_beta, the result will have a non-zeror total magnetization; if they are different, the result will tend towards a total magnetization of zero. Maybe the latter is right, since H2 should have a zero total magnetization and the latter has lower energy. The former is only a local miminum.

For QE, mixing_beta_mag is equal to mixing_beta, so they get a "total magnetization = 2.00 Bohr mag/cell" result.

this question is much more complicated and also depends on bond length, initial guess, mixing, etc., see https://mattermodeling.stackexchange.com/questions/11753/unrestricted-hartree-fock-density-matrices-initialization
I also copy one figure here:

[WHUweiqingzhou]

@ieiue, I agree with @Qianruipku. This change originates from PR #3226, where we update the mixing method for magnetic calculation. I think the states with -27 eV is not the ground one, and new mixing can find a true ground state.

[WHUweiqingzhou]

I close this issue due to no further discussion, feel free to reopen it.