fix: restore product-rule factor in LCAO DeltaSpin forces

source/source_lcao/module_operator_lcao/dspin_fs.cpp

@@ -215,12 +215,12 @@ void DeltaSpin<OperatorLCAO<TK, TR>>::cal_force_stress(const bool cal_force,
 #ifdef __MPI
         Parallel_Reduce::reduce_all(force.c, force.nr * force.nc);
 #endif
-        if (this->nspin != 4)
+        // cal_force_IJR differentiates only one overlap in B_mu * B_nu.
+        // The Hermitian partner contributes equally for every spin representation;
+        // this factor is from the product rule, not spin degeneracy.
+        for (int i = 0; i < force.nr * force.nc; i++)
         {
-            for (int i = 0; i < force.nr * force.nc; i++)
-            {
-                force.c[i] *= 2.0;
-            }
+            force.c[i] *= 2.0;
         }
     }
 

tests/03_NAO_multik/scf_deltaspin4/INPUT

@@ -1,32 +1,42 @@
 INPUT_PARAMETERS
 suffix            autotest
 calculation       scf
-basis_type        lcao
 ecutwfc           12
-gamma_only        0
-noncolin          1
-nspin             4
-#nbands           40
-scf_thr           1.0e-5
-scf_nmax          100
-out_chg           0
+scf_thr_type      2
+scf_thr           1e-10
+scf_ene_thr       1e-9
+scf_nmax          10
 smearing_method   gaussian
 smearing_sigma    0.01
 mixing_type       broyden
-mixing_beta       0.4
-ks_solver         scalapack_gvx
-symmetry          0
-cal_force         1
-cal_stress        1
-
-# DeltaSpin parameters
+mixing_beta       0.2
+mixing_ndim       12
+scf_os_stop       0
+scf_os_thr        -0.01
+scf_os_ndim       20
+basis_type        lcao
+kpar              1
+symmetry          -1
+noncolin          1
+nspin             4
+lspinorb          1
+onsite_radius     3.0
 sc_mag_switch     1
-sc_thr            5e-3
-nsc               20
+sc_thr            1e-10
+nsc               100
 nsc_min           2
 alpha_trial       0.01
 sccut             3.0
-sc_scf_thr        1e-2
-
+sc_scf_thr        10
+sc_os_ndim        5
+cal_stress        1
+cal_force         1
+gga_grad          2
+gamma_only        0
+out_mul           0
+sc_drop_thr       0
+test_force        1
+device            cpu
+ks_solver         scalapack_gvx
 pseudo_dir        ../../PP_ORB
 orbital_dir       ../../PP_ORB

tests/03_NAO_multik/scf_deltaspin4/README

@@ -1,47 +1,2 @@
-Test LCAO DeltaSpin with noncollinear spin (nspin=4), Z-only magnetization constraint, multi-k. Force and stress computed.
-
-## Known issue: final-step stress is not reproducible across MPI ranks (np >= 3)
-
-The SCF of this case never reaches scf_thr within scf_nmax=100 steps, so the
-test compares the stress of the *last, non-converged* SCF step. On that
-trajectory the transverse spin-density components (rho_mx, rho_my) start at
-round-off level (~1e-15) and are amplified exponentially (~x1.2 per step),
-because they are zero modes of the constrained problem. After ~100 steps the
-amplified noise perturbs the final stress at the O(1-100) kbar level.
-
-The seed of the noise is the ScaLAPACK solver (scalapack_gvx / pzhegvx):
-its internal global reductions are order-dependent on MPI message arrival
-timing, so for np >= 3 the eigenvectors within degenerate subspaces differ
-run-to-run at the 1e-16 level even for bitwise identical input. This was
-verified with a standalone pzhegvx reproducible test and is a property of
-the ScaLAPACK library, not of ABACUS. With OMP_NUM_THREADS > 1 the OpenMP
-reductions (e.g. grid-side charge sums) provide an additional noise seed
-with the same effect. np=2 with OMP_NUM_THREADS=1 is deterministic because
-a two-term floating-point sum is order-independent.
-
-Measured totalstress deviation vs the reference (generated at np=2):
-- np=2, OMP=1: deterministic, deviation 0.0 (bitwise reproducible).
-- np=4, OMP=1, Debian libscalapack 2.1: deterministic, deviation 0.093,
-  bitwise identical before and after the SpinConstrain refactor
-  (3a96cd744 vs dae8f0466), i.e. the refactor did not change this result.
-- np=4, OMP=1, CI toolchain ScaLAPACK: chaotic, deviation 3.124 observed.
-- np=4, OMP=14 (Autotest.sh default on a 56-core host): chaotic,
-  deviation 166 observed.
-- np=3: chaotic, deviations from 0.3 up to ~180 kbar, both before and
-  after the refactor.
-
-Therefore the per-case `threshold` file in this directory relaxes the
-thresholds for this case only. NOTE that Autotest.sh exits 1 (and fails
-CI) on ANY warning, so the warning thresholds, not only the fatal one,
-must cover the chaotic spread:
-- threshold (etot) = 1.0 eV, force_threshold = 10, stress_threshold = 500:
-  observed chaotic deviations are up to ~0.4 eV, ~5.3 and ~180
-  respectively; anything above these still prints a WARNING for review.
-- fatal_threshold = 1000: far above the chaos, but still catches
-  sign-flip-scale disasters (~1e5 kbar).
-
-Proper fixes (to be implemented later, in order of preference):
-1. Make the SCF converge (or reduce scf_nmax) so the compared state is a
-   converged one, which is insensitive to the noise; regenerate result.ref.
-2. Compare np=2 runs only (deterministic).
-3. Use a deterministic diagonalization setup (e.g. ELPA) for this case.
+Test LCAO DeltaSpin atomic forces with noncollinear spin (nspin=4) and SOC.
+Check the two-overlap product-rule force regression at 12 Ry with at most 10 SCF steps.

tests/03_NAO_multik/scf_deltaspin4/STRU

@@ -1,22 +1,21 @@
 ATOMIC_SPECIES
-Fe 1.000 Fe.upf
+Fe 55.845 Fe.upf
 
 NUMERICAL_ORBITAL
 Fe_gga_6au_100Ry_4s2p2d1f.orb
 
 LATTICE_CONSTANT
-8.190
+1.8897261254578281
 
 LATTICE_VECTORS
-1.00    0.50    0.50
-0.50    1.00    0.50
-0.50    0.50    1.00
+2.86 0 0
+0 2.86 0
+0 0 2.86
 
 ATOMIC_POSITIONS
-Direct
-
+Cartesian
 Fe
 0.0
 2
-0.00    0.00    0.00    magmom  0.0 0.0  2.0    sc 0 0 1
-0.51    0.51    0.51    magmom  0.0 0.0  -2.0   sc 0 0 1
+0.00000000000000 0.00000000000000 0.00000000000000 1 1 1 magmom 0.00000000000000 0.00000000000000 3.00000000000000 sc 1 1 1
+1.57300000000000 1.28700000000000 1.43000000000000 1 1 1 magmom 2.59807621135332 0.00000000000000 1.50000000000000 sc 1 1 1

tests/03_NAO_multik/scf_deltaspin4/result.ref

@@ -1,4 +1,5 @@
-etotref -6782.770693326635
-etotperatomref -3391.385346663318
-totalforceref 59.195649
-totalstressref 2526.777101
+etotref -6807.4629648003083275
+etotperatomref -3403.7314824002
+totalforceref 93.685868
+totalstressref 30117.121565
+totaltimeref 27.24

tests/03_NAO_multik/scf_deltaspin4/threshold

@@ -1,20 +1 @@
-# Per-case thresholds for this test only (see README for rationale).
-#
-# The final-step stress of this non-converged 100-step SCF trajectory is
-# chaotic across MPI ranks (np >= 3) and across OpenMP thread counts,
-# due to ScaLAPACK pzhegvx run-to-run non-determinism amplified by the
-# SCF. Measured totalstress deviations: 3.124 (CI, np=4, OMP=1), 166
-# (local np=4, OMP=14), up to ~180 (local np=3). etot (~0.4 eV) and
-# totalforce (up to ~5.3 eV/Angstrom) are less affected but not clean.
-#
-# NOTE: Autotest.sh exits 1 (CI failure) on ANY warning, so the warning
-# thresholds must also cover the chaotic spread. fatal_threshold is kept
-# far above the chaos but still catches sign-flip-scale disasters
-# (~1e5 kbar).
-#
-# default: threshold 1e-7, force_threshold 1e-4, stress_threshold 1e-3,
-#          fatal_threshold 1
-threshold         1.0
-force_threshold   10.0
-stress_threshold  500.0
-fatal_threshold   1000.0
+force_threshold 0.000001