Fix #5021 #5939: warn on ignored x/y magnetization and remove nspin=4 autoset

.ci/slurm/config.ini

@@ -491,7 +491,7 @@ runner = autotest_gpu
 
 [case.080]
 suite = 01_PW
-name = 035_PW_15_SO
+name = scf_spin4
 resource = pw_gpu1
 runner = autotest_gpu
 

docs/advanced/input_files/input-main.md

@@ -1365,6 +1365,8 @@
   - 1: Spin degeneracy
   - 2: Collinear spin polarized.
   - 4: Noncollinear or spin-orbit calculations. Set nspin to 4 explicitly when noncolin or lspinorb is enabled.
+  - Note: With nspin=2 and no initial magnetization in STRU, a moment of 1.0 is autoset for every atom, unless symmetry is 1.
+  - Note: With nspin=4 no moment is ever autoset. The calculation starts from zero magnetic moment and a warning is printed; set 'mag' explicitly in STRU for the magnetic atoms if a magnetic ground state is expected.
 - **Default**: 1
 
 ### gga_grad
@@ -1606,6 +1608,8 @@
   - noncolin=0, lspinorb=1: SOC with z-axis magnetism only (for non-magnetic materials with SOC)
   - noncolin=1, lspinorb=0: Non-collinear magnetism without SOC
   - noncolin=1, lspinorb=1: Both non-collinear magnetism and SOC
+  - Note: When nspin=4 and noncolin=0, only the z component of the initial magnetization in STRU is used; x/y components are ignored and a warning is printed.
+  - Note: When nspin=4 and no initial magnetization is set in STRU, the calculation starts from zero magnetic moment; no automatic magnetization is assigned.
 - **Default**: False
 
 ### soc_lambda

docs/advanced/input_files/stru.md

@@ -281,7 +281,7 @@ For general usage requirements, the APNSv1.0 pseudopotential and orbital set is
       0.0 0.0 0.0 m 0 0 0
       0.5 0.5 0.5 m 1 1 1
       ```
-      For `nspin==2`, we will autoset atomic magmon is `1.0`:
+      For `nspin==2`, the atomic magnetic moment is autoset to `1.0`:
       ```
       Fe
       1.0
@@ -295,21 +295,32 @@ For general usage requirements, the APNSv1.0 pseudopotential and orbital set is
       0.0 0.0 0.0 m 0 0 0
       0.5 0.5 0.5 m 1 1 1
       ```
-      For `nspin==4`, we will autoset atomic magmon as follow:
+      This autoset is skipped as soon as `STRU` specifies a finite moment for
+      any single atom, and also when [`symmetry`](input-main.md#symmetry) is
+      `1`, because an all-zero moment is a legitimate nonmagnetic choice under
+      the full point group.
+
+      For `nspin==4`, **no** moment is autoset: the calculation starts from a
+      zero magnetic moment and a warning is printed. Zero magnetization is a
+      stationary point, so such a run stays nonmagnetic. Set `mag` explicitly
+      on the magnetic atoms if a magnetic ground state is expected, e.g.:
       ```
       Fe
       0.0
       2
-      0.0 0.0 0.0 m 0 0 0 mag 1 1 1
-      0.5 0.5 0.5 m 1 1 1 mag 1 1 1
+      0.0 0.0 0.0 m 0 0 0 mag 0.0 0.0 2.0
+      0.5 0.5 0.5 m 1 1 1 mag 0.0 0.0 2.0
 
       O
       0.0
       2
-      0.0 0.0 0.0 m 0 0 0 mag 1 1 1
-      0.5 0.5 0.5 m 1 1 1 mag 1 1 1
+      0.0 0.0 0.0 m 0 0 0
+      0.5 0.5 0.5 m 1 1 1
       ```
-      However, this autoset will not be vaild once `STRU` specalize a finite magnetic for any single atom.
+      With [`noncolin`](input-main.md#noncolin) `0`, only the z component of
+      `mag` is used; the x/y components are discarded and the number of
+      affected atoms is reported per atom type. A moment given purely along
+      x/y therefore also leaves a zero starting moment.
 
   - `lambda`: Lagrange multiplier vector for spin constraint method. Can specify one value (z-component) or three values for x, y, z components (e.g., `lambda 0.5` or `lambda 0.1 0.2 0.3`). Values are in eV and will be converted to Rydberg internally. Used with spin-constrained DFT (enable with `sc_mag_switch` in INPUT file).
 

docs/parameters.yaml

@@ -709,6 +709,8 @@ parameters:
       * 1: Spin degeneracy
       * 2: Collinear spin polarized.
       * 4: Noncollinear or spin-orbit calculations. Set nspin to 4 explicitly when noncolin or lspinorb is enabled.
+      * Note: With nspin=2 and no initial magnetization in STRU, a moment of 1.0 is autoset for every atom, unless symmetry is 1.
+      * Note: With nspin=4 no moment is ever autoset. The calculation starts from zero magnetic moment and a warning is printed; set 'mag' explicitly in STRU for the magnetic atoms if a magnetic ground state is expected.
     default_value: "1"
     unit: ""
     availability: ""
@@ -996,6 +998,8 @@ parameters:
         * noncolin=0, lspinorb=1: SOC with z-axis magnetism only (for non-magnetic materials with SOC)
         * noncolin=1, lspinorb=0: Non-collinear magnetism without SOC
         * noncolin=1, lspinorb=1: Both non-collinear magnetism and SOC
+      * Note: When nspin=4 and noncolin=0, only the z component of the initial magnetization in STRU is used; x/y components are ignored and a warning is printed.
+      * Note: When nspin=4 and no initial magnetization is set in STRU, the calculation starts from zero magnetic moment; no automatic magnetization is assigned.
     default_value: "False"
     unit: ""
     availability: ""

source/source_cell/cal_ux.cpp

@@ -2,6 +2,19 @@
 
 namespace unitcell {
 
+// Compute the global spin quantization axis ux_ for non-collinear (nspin=4)
+// calculations. ux_ is a single axis shared by the whole system (NOT one axis
+// per atom); it is consumed by XC_Functional::noncolin_rho to split the spin
+// density matrix into spin-up/down scalars for the LSDA/GGA kernel.
+//
+// The axis is derived from the *initial* moments m_loc_ read from STRU
+// (direction of the first atom with |m| > threshold), not from the converged
+// density. Consequently the STRU initial moments silently select the XC
+// projection axis: with all-zero initial moments lsign_ stays false and the
+// XC projection becomes direction-blind (neg = +1 everywhere in
+// noncolin_rho), which can follow a different SCF path than a run started
+// with an explicit non-zero moment even when the starting density (e.g. from
+// init_chg = wfc/file) is identical.
 void cal_ux(UnitCell& ucell, const int nspin) {
 
     if (nspin != 4)

source/source_cell/read_atoms.cpp

@@ -70,6 +70,11 @@ bool unitcell::read_atom_positions(UnitCell& ucell,
             int na = ucell.atoms[it].na;
             ucell.nat += na;
 
+            // Number of atoms of this type whose x/y moment was discarded by
+            // the collinear projection (nspin=4 with noncolin=0). Reported
+            // once per type after the atom loop instead of once per atom.
+            int n_xy_discarded = 0;
+
             if (na > 0)
             {
                 unitcell::allocate_atom_properties(ucell.atoms[it], na);
@@ -103,9 +108,14 @@ bool unitcell::read_atom_positions(UnitCell& ucell,
                     }
 
                     // Process magnetization
-                    unitcell::process_magnetization(ucell.atoms[it], it, ia, nspin,
-                                        input_vec_mag, input_angle_mag, ofs_running,
-                                        noncolin);
+                    const bool xy_discarded
+                        = unitcell::process_magnetization(ucell.atoms[it], it, ia, nspin,
+                                            input_vec_mag, input_angle_mag, ofs_running,
+                                            noncolin);
+                    if (xy_discarded)
+                    {
+                        ++n_xy_discarded;
+                    }
 
                     // Transform coordinates
                     unitcell::transform_atom_coordinates(ucell.atoms[it], ia, Coordinate,
@@ -117,6 +127,9 @@ bool unitcell::read_atom_positions(UnitCell& ucell,
                     ucell.atoms[it].dis[ia].set(0, 0, 0);
                 }//endj
             }    // end na
+
+            unitcell::warn_xy_magnetization_ignored(ucell.atoms[it].label,
+                                        n_xy_discarded, na, ofs_running);
             // reset some useless parameters
             if (set_element_mag_zero)
             {
@@ -130,15 +143,29 @@ bool unitcell::read_atom_positions(UnitCell& ucell,
         // so do not override it with an autoset seed. Warn instead.
         if (symmetry == 1)
         {
-            ofs_running << "\n WARNING: initial magmom is all zero and symmetry=1; "
-                        << "autoset magnetism is SKIPPED to preserve the symmetry of the initial (nonmagnetic) structure.\n"
-                        << "          If spontaneous magnetism is expected, set magmom explicitly "
-                        << "in STRU, or use symmetry = 0 or -1." << std::endl;
+            // Only report the skipped autoset when there is in fact
+            // nothing set: this message used to be printed even for an
+            // explicitly magnetized STRU.
+            if (unitcell::is_magnetization_all_zero(ucell))
+            {
+                ofs_running << "\n WARNING: initial magmom is all zero and symmetry=1; "
+                            << "autoset magnetism is SKIPPED to preserve the symmetry of the initial (nonmagnetic) structure.\n"
+                            << "          If spontaneous magnetism is expected, set magmom explicitly "
+                            << "in STRU, or use symmetry = 0 or -1." << std::endl;
+            }
         }
         else
         {
             unitcell::autoset_magnetization(ucell, nspin, ofs_running);
         }
+
+        // nspin=4 never autosets a moment, so the zero-moment start has to
+        // be reported for every symmetry setting, not only on the autoset
+        // path that symmetry=1 skips.
+        if (nspin == 4)
+        {
+            unitcell::warn_zero_magnetization_nspin4(ucell, ofs_running);
+        }
     }   // end scan_begin
 
     // Final validation and output

source/source_cell/read_atoms_helper.cpp

@@ -17,6 +17,9 @@ namespace {
     constexpr char LOWER_A = 'a';
     constexpr char LOWER_Z = 'z';
     constexpr char MINUS_SIGN = '-';
+
+    // Below this value a magnetic moment component is treated as zero.
+    constexpr double mag_threshold = 1e-5;
 }
 
 namespace unitcell {
@@ -79,56 +82,97 @@ void set_atom_movement_flags(Atom& atom, int ia,
     }
 }
 
-void autoset_magnetization(UnitCell& ucell, int nspin,
-                           std::ofstream& ofs_running)
+bool is_magnetization_all_zero(const UnitCell& ucell)
 {
-    const int ntype = ucell.ntype;
-
-    // Check if any atom has non-zero magnetization
-    int autoset_mag = 1;
-    for (int it = 0; it < ntype; it++)
+    for (int it = 0; it < ucell.ntype; it++)
     {
         for (int ia = 0; ia < ucell.atoms[it].na; ia++)
         {
-            if(std::abs(ucell.atoms[it].mag[ia]) > 1e-5)
+            if (std::abs(ucell.atoms[it].mag[ia]) > mag_threshold)
             {
-                autoset_mag = 0;
-                break;
+                return false;
             }
         }
     }
+    return true;
+}
+
+void autoset_magnetization(UnitCell& ucell, int nspin,
+                           std::ofstream& ofs_running)
+{
+    // nspin=4 never autosets a moment. An all-zero moment in STRU is a
+    // legitimate nonmagnetic start; warn_zero_magnetization_nspin4() reports
+    // it instead, independently of the symmetry switch that gates this call.
+    if (nspin != 2)
+    {
+        return;
+    }
 
-    if (autoset_mag)
+    if (!is_magnetization_all_zero(ucell))
+    {
+        return;
+    }
+
+    const int ntype = ucell.ntype;
+    for (int it = 0; it < ntype; it++)
     {
-        if(nspin==4)
+        for (int ia = 0; ia < ucell.atoms[it].na; ia++)
         {
-            for (int it = 0; it < ntype; it++)
-            {
-                for (int ia = 0; ia < ucell.atoms[it].na; ia++)
-                {
-                    ucell.atoms[it].m_loc_[ia].x = 1.0;
-                    ucell.atoms[it].m_loc_[ia].y = 1.0;
-                    ucell.atoms[it].m_loc_[ia].z = 1.0;
-                    ucell.atoms[it].mag[ia] = sqrt(pow(ucell.atoms[it].m_loc_[ia].x,2)
-                            +pow(ucell.atoms[it].m_loc_[ia].y,2)
-                            +pow(ucell.atoms[it].m_loc_[ia].z,2));
-                    ModuleBase::GlobalFunc::OUT(ofs_running,"Autoset magnetism for this atom", 1.0, 1.0, 1.0);
-                }
-            }
+            ucell.atoms[it].mag[ia] = 1.0;
+            ucell.atoms[it].m_loc_[ia].z = ucell.atoms[it].mag[ia];
+            ModuleBase::GlobalFunc::OUT(ofs_running,"Autoset magnetism for this atom", 1.0);
         }
-        else if(nspin==2)
+    }
+}
+
+void warn_zero_magnetization_nspin4(const UnitCell& ucell,
+                                    std::ofstream& ofs_running)
+{
+    const double threshold_sq = mag_threshold * mag_threshold;
+    for (int it = 0; it < ucell.ntype; it++)
+    {
+        for (int ia = 0; ia < ucell.atoms[it].na; ia++)
         {
-            for (int it = 0; it < ntype; it++)
+            const ModuleBase::Vector3<double>& m = ucell.atoms[it].m_loc_[ia];
+            const double norm_sq = m.x * m.x + m.y * m.y + m.z * m.z;
+            if (norm_sq > threshold_sq)
             {
-                for (int ia = 0; ia < ucell.atoms[it].na; ia++)
-                {
-                    ucell.atoms[it].mag[ia] = 1.0;
-                    ucell.atoms[it].m_loc_[ia].z = ucell.atoms[it].mag[ia];
-                    ModuleBase::GlobalFunc::OUT(ofs_running,"Autoset magnetism for this atom", 1.0);
-                }
+                return;
             }
         }
     }
+
+    const std::string summary = "nspin=4 but no initial magnetization is set in STRU; "
+        "all atoms start from zero magnetic moment.";
+    std::cout << " Warning: " << summary << std::endl;
+    ModuleBase::WARNING("read_atom_positions", summary);
+    ofs_running << "\n WARNING: " << summary << std::endl;
+    ofs_running << "          If a magnetic ground state is expected, set 'mag' explicitly "
+                << "in STRU for the magnetic atoms." << std::endl;
+    ofs_running << "          With noncolin=0 the x/y components of 'mag' are discarded, so a "
+                << "moment given purely along x/y also leaves a zero starting moment."
+                << std::endl;
+}
+
+void warn_xy_magnetization_ignored(const std::string& label,
+                                   int n_discarded,
+                                   int na,
+                                   std::ofstream& ofs_running)
+{
+    if (n_discarded <= 0)
+    {
+        return;
+    }
+
+    std::stringstream ss;
+    ss << "atom type " << label << ": " << n_discarded << " of " << na
+       << " atoms give non-zero x/y magnetization in STRU, but nspin=4 with noncolin=0 "
+       << "is a collinear calculation; only the z component is used. "
+       << "x/y components are IGNORED. Set 'noncolin 1' to use the full vector.";
+    const std::string summary = ss.str();
+    std::cout << " Warning: " << summary << std::endl;
+    ModuleBase::WARNING("read_atom_positions", summary);
+    ofs_running << "\n WARNING: " << summary << std::endl;
 }
 
 bool finalize_atom_positions(UnitCell& ucell,
@@ -265,12 +309,14 @@ void transform_atom_coordinates(Atom& atom, int ia,
     }
 }
 
-void process_magnetization(Atom& atom, int it, int ia,
+bool process_magnetization(Atom& atom, int it, int ia,
                           int nspin, bool input_vec_mag,
                           bool input_angle_mag,
                           std::ofstream& ofs_running,
                           const bool noncolin)
 {
+    bool xy_discarded = false;
+
     // Recalculate mag and m_loc_ from read in angle1, angle2 and mag or mx, my, mz
     if(input_angle_mag)
     {
@@ -307,6 +353,23 @@ void process_magnetization(Atom& atom, int it, int ia,
         if(!noncolin)
         {
             // collinear case with nspin = 4, only z component is used
+            // Note: a scalar "mag <value>" in STRU is stored purely in
+            // m_loc_.z (x/y stay zero), so it passes the check below without
+            // flagging a discard and still defines a non-zero z moment.
+            // cal_ux() then turns that into the global XC quantization axis
+            // ux_ = (0,0,1) with lsign_ = true. Any positive scalar suffices;
+            // only the *direction* enters ux_, the magnitude (e.g. 1.7320508
+            // vs 1.0) is irrelevant.
+            //
+            // The discard is only reported back to the caller, which
+            // aggregates it per atom type: emitting a message here would
+            // print one warning per atom for a cell where every atom tilts
+            // the moment.
+            if(std::abs(atom.m_loc_[ia].x) > mag_threshold
+                || std::abs(atom.m_loc_[ia].y) > mag_threshold)
+            {
+                xy_discarded = true;
+            }
             atom.m_loc_[ia].x = 0;
             atom.m_loc_[ia].y = 0;
         }
@@ -347,6 +410,8 @@ void process_magnetization(Atom& atom, int it, int ia,
             ModuleBase::GlobalFunc::OUT(ofs_running, ss.str(),atom.mag[ia]);
         }
     }
+
+    return xy_discarded;
 }
 
 bool parse_atom_properties(std::ifstream& ifpos,