LCAO non-collinear (nspin=4, no SOC) SCF magnetization collapses in v3.11: σy sign mismatch between H assembly (clx_j) and LCAO density convention

[fukatou]

Describe the bug

ABACUS v3.11.0-beta7 (085b91e) fails to converge in LCAO non-collinear (nspin=4, noncolin=1, no SOC) spin-polarized SCF calculations: the density residual DRHO stalls at ~0.004–0.015 while the total magnetization collapses (e.g. AMAG 12.5→0.77 within the first 3 iterations with genelpa; 12.0→3.5 at iteration 9 with scalapack_gvx). This is a systematic numerical error introduced by the v3.11 gint refactor: the identical input converges in v3.10.1/v3.10.0 in 24 iterations (AFM, ±3.27 μB per Cr, moments along y).

Root cause (source-level, verified by a build-level experiment): in source/source_lcao/module_gint/gint_common.cpp (merge_hr_part_to_hR), the 2×2 spin-block assembly was changed to the textbook Pauli convention (clx_j = {0,-1,+1,0}, i.e. H↑↓ = Bx − i·By, plus conjugate for the lower triangle), but the LCAO density-matrix convention was not updated in sync: func_xyz_to_updown in source/source_estate/module_dm/density_matrix.cpp still defines the y-magnetization density as rho_y = −2·Im(ρ↑↓) (opposite to the PW path +2·Im(ρ↑↓) in elecstate_op.cpp). The σy feedback loop in the LCAO path therefore has the wrong sign (my ∝ −my), driving the y-magnetization to zero. Only the y-component is affected (x/z channels use identical conventions in both versions), which matches the observed collapse pattern.

Expected behavior

The SCF calculation should converge to the same physical state as v3.10.x: AFM state with moments ±3.27 μB along y, total energy −11710.1785 eV (24 iterations), instead of collapsing the magnetization and stalling.

To Reproduce

1. Clone the source code from https://github.com/deepmodeling/abacus-develop.git (v3.11.0-beta7, commit 085b91e)
2. Configure CMake with options -DCMAKE_CXX_COMPILER=mpiicpc -DCMAKE_C_COMPILER=mpiicc -DCMAKE_Fortran_COMPILER=mpiifort -DENABLE_LCAO=ON -DENABLE_MPI=ON (oneAPI 2022.1)
3. Build ABACUS with cmake --build build -j16
4. Run ABACUS with mpirun -n 48 ./abacus using the attached case (bi-layer CrSBr, SG15 NC pseudopotentials, SR numerical atomic orbitals, 12 atoms) with the following key parameters:
   • nspin=4 noncolin=1 lspinorb=0, basis_type=lcao, ecutwfc=100
   • vdw_method=d3_bj, scf_thr=1e-8, scf_nmax=300
   • mixing_type=pulay, mixing_beta=0.2, mixing_beta_mag=0.2, mixing_ndim=12
   • ks_solver=scalapack_gvx (also reproduced with genelpa)
   • Gaussian smearing, σ = 1 meV; k-grid 15×11×1 (Gamma centered)
   • STRU initial moments mag 0.0 5.0 0.0 / mag 0.0 -5.0 0.0 (AFM along y)

Observed iteration log (v3.11, unpatched, scalapack_gvx, run stopped manually at 150 iters — never converged):

GV1   DRHO=3.6e-01  AMAG=14.4   E=-11725.51
GV2   DRHO=2.8e-01  AMAG=13.0   E=-11740.05
...
GV8   DRHO=1.8e-01  AMAG=12.0   E=-11717.24
GV9   DRHO=5.7e-02  AMAG=3.5    E=-11703.90   <- magnetization collapses
GV40  DRHO=2.7e-03  AMAG=0.03   E=-11703.63
GV80  DRHO=1.2e-03  AMAG=0.07   E=-11703.63
GV120 DRHO=5.6e-03  AMAG=0.25   E=-11703.63
GV150 DRHO=4.8e-03  AMAG=0.22   E=-11703.63   <- stopped, DRHO stalled ~5e-3

With genelpa the magnetization does not simply collapse to zero but oscillates (AMAG 0.77 → 3.1 → 13.5 over 132 iterations) with DRHO stalled at ~0.1 — also never converging.

Case files (attached inline, self-contained):

INPUT:

INPUT_PARAMETERS
calculation          scf
pseudo_dir           <pseudo_dir>
orbital_dir          <orbital_dir>
basis_type           lcao
ecutwfc              100
scf_nmax             300
scf_thr              1.0e-8
nspin                4
noncolin             1
lspinorb             0
smearing_method      gaussian
smearing_sigma       0.00007350
mixing_type          pulay
mixing_beta          0.2
mixing_beta_mag      0.2
mixing_ndim          12
ks_solver            scalapack_gvx
out_freq_elec        10
vdw_method           d3_bj
out_mul              1

STRU:

ATOMIC_SPECIES
Cr1  51.996  Cr_ONCV_PBE-1.2.upf
Cr2  51.996  Cr_ONCV_PBE-1.2.upf
S    32.065  S_ONCV_PBE-1.2.upf
Br   79.904  Br_ONCV_PBE-1.2.upf
NUMERICAL_ORBITAL
Cr_gga_8au_100Ry_4s2p2d1f.orb
Cr_gga_8au_100Ry_4s2p2d1f.orb
S_gga_7au_100Ry_2s2p1d.orb
Br_gga_7au_100Ry_2s2p1d.orb
LATTICE_CONSTANT
1.8897261246
LATTICE_VECTORS
3.4949676754 0.0000000000 0.0000000000
0.0000000000 4.7138923664 0.0000000000
0.0000000000 0.0000000000 27.4286647980
ATOMIC_POSITIONS
Direct
Cr1
0.0
2
0.5000000000 0.0000000000 0.3201256800  mag 0.0 5.0 0.0  0 0 0
0.0000000000 0.5000000000 0.3926176240  mag 0.0 5.0 0.0  0 0 0
Cr2
0.0
2
0.5000000000 0.0000000000 0.6073823760  mag 0.0 -5.0 0.0  0 0 0
0.0000000000 0.5000000000 0.6798743200  mag 0.0 -5.0 0.0  0 0 0
S
0.0
4
0.5000000000 0.5000000000 0.3343522027  0 0 0
0.0000000000 0.0000000000 0.3783911012  0 0 0
0.5000000000 0.5000000000 0.6216088988  0 0 0
0.0000000000 0.0000000000 0.6656477973  0 0 0
Br
0.0
4
0.0000000000 0.0000000000 0.2552074646  0 0 0
0.5000000000 0.5000000000 0.4575358393  0 0 0
0.0000000000 0.0000000000 0.5424641607  0 0 0
0.5000000000 0.5000000000 0.7447925354  0 0 0

KPT:

K_POINTS
0
Gamma
15 11 1 0 0 0

Pseudopotentials: SG15 NC v1.2 (Cr_ONCV_PBE-1.2.upf, S_ONCV_PBE-1.2.upf, Br_ONCV_PBE-1.2.upf); numerical atomic orbitals (SR): Cr_gga_8au_100Ry_4s2p2d1f.orb, S_gga_7au_100Ry_2s2p1d.orb, Br_gga_7au_100Ry_2s2p1d.orb. Replace <pseudo_dir>/<orbital_dir> in INPUT with the corresponding local directories.

Environment

• OS: RHEL-like cluster OS (Linux)
• Compiler: Intel oneAPI 2022.1 (mpiicpc/mpiicc/mpiifort)
• Dependencies: Intel MKL 2022.1.0, Intel MPI 2021.9.0
• Machine: 2× Intel Xeon (48 cores), also reproduced on a different cluster (same toolchain)

Additional Context

Source-level evidence:

v3.10.1 source/module_hamilt_lcao/module_gint/gint_k_pvpr.cpp (transfer_pvpR, NSPIN=4):

std::vector<int> clx_i = {1, 0, 0, -1};
std::vector<int> clx_j = {0, 1, -1, 0};     // H↑↓ = Bx + i·By
...
lower_mat->get_value(icol, irow) = upper_mat->get_value(irow, icol);   // no conj

v3.11 source/source_lcao/module_gint/gint_common.cpp (merge_hr_part_to_hR):

// Pauli-to-spinor conversion: H = V_0*I + B_x*sigma_x + B_y*sigma_y + B_z*sigma_z
// sigma_y = [[0,-i],[i,0]], so H_{up,down} = B_x - i*B_y ...
std::vector<int> clx_i = {1, 0, 0, -1};
std::vector<int> clx_j = {0, -1, 1, 0};     // H↑↓ = Bx − i·By (textbook Pauli)
...
lower_mat->get_value(icol, irow) = std::conj(upper_mat->get_value(irow, icol));  // conj

Note: the two versions' H↓↑(−R) are actually identical (Bx + i·By) because the added conj compensates the clx_j flip; the whole difference reduces to the sign of the σy (B_y) coefficient in H↑↓(R).

The density/potential conventions are identical in both versions (not updated in sync):

• LCAO DM→rho (density_matrix.cpp, func_xyz_to_updown): rho_y = -Im(rho_updown) + Im(rho_downup) = -2·Im(rho_updown)
• PW psi→rho (elecstate_op.cpp): rho_y = +2·Im(psi_up* psi_dn)
• XC magnetic potential (xc_pot.cpp): v(ipol) = e2 * vs * rho[ipol] / amag (B ∝ rho, no extra sign)

Self-consistency analysis of the σy channel (LCAO path, rho[2] = −my):

	v3.10.1 (H↑↓ = Bx + iBy)	v3.11 (H↑↓ = Bx − iBy)
σy coefficient in H	−B_y ∝ +my	+B_y ∝ −my
ground state my ∝	+my → stable	−my → my = 0, collapse

This is consistent with all observations: z-channel (rho[3], σz, clx_i) and the PW path are unaffected in both versions (v3.11 PW+NC with z-initialized moments converges to a z-AFM state, ±2.61 μB/Cr); SOC (lspinorb=1) calculations are unaffected (different operator path).

Verification (build-level experiment):

Patching only line 179 of gint_common.cpp:

std::vector<int> clx_j = {0, -1, 1, 0};  ->  std::vector<int> clx_j = {0, 1, -1, 0};

(keeping the conj in the lower-triangle fill), recompiling and rerunning the identical input:

Quantity	v3.10.0 reference	v3.11 + patch	Match
Total energy (eV)	−11710.1785	−11710.178525	✅ (2.5e-5)
AMAG (μB)	13.5	13.5 (stable throughout)	✅
Cr1–4 moments (μB)	(0, +3.2733, 0) / (0, +3.2711, 0) / (0, −3.2710, 0) / (0, −3.2733, 0)	identical, digit by digit	✅
Iterations to converge	24	23 (DRHO 1.25e-8)	✅

Suggested fixes:

• Option A (recommended — unify conventions): in func_xyz_to_updown, change the y-component to match the PW path:

  // rho_y = +Im(rho_updown) - Im(rho_downup)  (= +2·Im(rho_updown))
  target_DMR_mat[icol + step_trace[2]] = tmp[1].imag() - tmp[2].imag();

  keeping the standard-Pauli H assembly. This unifies the LCAO and PW conventions and also fixes the latent PW y-channel issue in v3.10.x.
• Option B (minimal, verified by this report): restore clx_j = {0, 1, -1, 0} in merge_hr_part_to_hR (keeping the conj), i.e. keep the v3.10.1 assembly, which is the correct companion of the LCAO density convention.

Task list for Issue attackers (only for developers)

• Reproduce the bug with the provided input case
• Confirm the σy sign mismatch (clx_j vs func_xyz_to_updown)
• Apply Option A or B and verify against the v3.10 reference values above