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

docs/advanced/input_files/input-main.md

@@ -1604,6 +1604,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/parameters.yaml

@@ -996,6 +996,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_helper.cpp

@@ -102,19 +102,14 @@ void autoset_magnetization(UnitCell& ucell, int nspin,
     {
         if(nspin==4)
         {
-            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);
-                }
-            }
+            // nspin=4 with all-zero magnetization: do NOT autoset,
+            // start from zero and let the user decide explicitly.
+            std::string msg = "nspin=4 but no initial magnetization is set in STRU.\n"
+                "  All atoms start from zero magnetic moment.\n"
+                "  If a magnetic ground state is expected, set 'mag' explicitly "
+                "in STRU for the magnetic atoms.";
+            std::cout << " Warning: " << msg << std::endl;
+            ModuleBase::GlobalFunc::OUT(ofs_running, "Warning", msg);
         }
         else if(nspin==2)
         {
@@ -307,6 +302,25 @@ 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
+            // triggering the warning 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.
+            if(std::abs(atom.m_loc_[ia].x) > 1e-5 || std::abs(atom.m_loc_[ia].y) > 1e-5)
+            {
+                std::stringstream ss;
+                ss << "Atom " << ia+1 << " of type " << atom.label
+                   << ": STRU gives non-zero x/y magnetization ("
+                   << atom.m_loc_[ia].x << ", " << atom.m_loc_[ia].y << ", "
+                   << atom.m_loc_[ia].z << "), but nspin=4 with noncolin=0 "
+                   << "is a collinear calculation; only the z component is used.\n"
+                   << "  x/y components are IGNORED. Set 'noncolin 1' to use the full vector.";
+                std::cout << " Warning: " << ss.str() << std::endl;
+                ModuleBase::GlobalFunc::OUT(ofs_running, "Warning", ss.str());
+            }
             atom.m_loc_[ia].x = 0;
             atom.m_loc_[ia].y = 0;
         }

source/source_cell/test/CMakeLists.txt

@@ -173,7 +173,43 @@ add_test(NAME MODULE_CELL_parallel_kpoints_test
 AddTest(
   TARGET MODULE_CELL_unitcell_test
   LIBS base device cell_info symmetry
-  SOURCES unitcell_test.cpp ../cal_ux.cpp
+  SOURCES unitcell_test.cpp
+)
+
+AddTest(
+  TARGET MODULE_CELL_cal_ux_test
+  LIBS base device cell_info
+  SOURCES test_cal_ux.cpp ../cal_ux.cpp
+)
+
+AddTest(
+  TARGET MODULE_CELL_read_orb_test
+  LIBS base device cell_info
+  SOURCES test_read_orb.cpp
+)
+
+AddTest(
+  TARGET MODULE_CELL_update_cell_test
+  LIBS base device cell_info
+  SOURCES test_update_cell.cpp
+)
+
+AddTest(
+  TARGET MODULE_CELL_cell_tools_test
+  LIBS base device cell_info
+  SOURCES test_cell_tools.cpp
+)
+
+AddTest(
+  TARGET MODULE_CELL_read_atom_species_test
+  LIBS base device cell_info
+  SOURCES test_read_atom_species.cpp
+)
+
+AddTest(
+  TARGET MODULE_CELL_read_atoms_test
+  LIBS base device cell_info
+  SOURCES test_read_atoms.cpp
 )
 
 AddTest(

source/source_cell/test/read_atoms_hlp_test.cpp

@@ -774,3 +774,207 @@ int main(int argc, char **argv)
     ::testing::InitGoogleTest(&argc, argv);
     return RUN_ALL_TESTS();
 }
+
+// ============================================================
+// Tests for issue #5021: nspin=4 with noncolin=0 should warn
+// when x/y components are non-zero
+// ============================================================
+
+TEST_F(ReadAtomsHelperTest, ProcessMagnetizationNspin4CollinearWarnsOnXY)
+{
+    Atom atom;
+    atom.label = "Fe";
+    atom.mag.resize(1);
+    atom.m_loc_.resize(1);
+    atom.angle1.resize(1);
+    atom.angle2.resize(1);
+
+    atom.m_loc_[0].set(1.0, 0.5, 0.3);
+    atom.mag[0] = sqrt(1.0 + 0.25 + 0.09);
+
+    const int nspin = 4;
+    const bool input_vec_mag = true;
+    const bool input_angle_mag = false;
+    const bool noncolin = false;
+
+    // Capture stdout
+    testing::internal::CaptureStdout();
+    unitcell::process_magnetization(atom, 0, 0, nspin, input_vec_mag, input_angle_mag, ofs_running, noncolin);
+    std::string stdout_output = testing::internal::GetCapturedStdout();
+
+    // x/y should be zeroed
+    EXPECT_DOUBLE_EQ(atom.m_loc_[0].x, 0.0);
+    EXPECT_DOUBLE_EQ(atom.m_loc_[0].y, 0.0);
+    EXPECT_DOUBLE_EQ(atom.m_loc_[0].z, 0.3);
+
+    // Warning should appear on stdout
+    EXPECT_NE(stdout_output.find("Warning"), std::string::npos);
+    EXPECT_NE(stdout_output.find("x/y components are IGNORED"), std::string::npos);
+
+    // Warning should appear in running log
+    ofs_running.flush();
+    std::ifstream ifs("test_running.log");
+    std::stringstream buffer;
+    buffer << ifs.rdbuf();
+    std::string log_content = buffer.str();
+    ifs.close();
+    EXPECT_NE(log_content.find("Warning"), std::string::npos);
+    EXPECT_NE(log_content.find("x/y components are IGNORED"), std::string::npos);
+}
+
+TEST_F(ReadAtomsHelperTest, ProcessMagnetizationNspin4CollinearNoWarnOnZeroXY)
+{
+    Atom atom;
+    atom.label = "Fe";
+    atom.mag.resize(1);
+    atom.m_loc_.resize(1);
+    atom.angle1.resize(1);
+    atom.angle2.resize(1);
+
+    atom.m_loc_[0].set(0.0, 0.0, 0.3);
+    atom.mag[0] = 0.3;
+
+    const int nspin = 4;
+    const bool input_vec_mag = true;
+    const bool input_angle_mag = false;
+    const bool noncolin = false;
+
+    testing::internal::CaptureStdout();
+    unitcell::process_magnetization(atom, 0, 0, nspin, input_vec_mag, input_angle_mag, ofs_running, noncolin);
+    std::string stdout_output = testing::internal::GetCapturedStdout();
+
+    EXPECT_DOUBLE_EQ(atom.m_loc_[0].x, 0.0);
+    EXPECT_DOUBLE_EQ(atom.m_loc_[0].y, 0.0);
+    EXPECT_DOUBLE_EQ(atom.m_loc_[0].z, 0.3);
+
+    // No warning should appear
+    EXPECT_EQ(stdout_output.find("Warning"), std::string::npos);
+}
+
+TEST_F(ReadAtomsHelperTest, ProcessMagnetizationNspin4NoncolinNoWarn)
+{
+    Atom atom;
+    atom.label = "Fe";
+    atom.mag.resize(1);
+    atom.m_loc_.resize(1);
+    atom.angle1.resize(1);
+    atom.angle2.resize(1);
+
+    atom.m_loc_[0].set(1.0, 0.5, 0.3);
+    atom.mag[0] = sqrt(1.0 + 0.25 + 0.09);
+
+    const int nspin = 4;
+    const bool input_vec_mag = true;
+    const bool input_angle_mag = false;
+    const bool noncolin = true;
+
+    testing::internal::CaptureStdout();
+    unitcell::process_magnetization(atom, 0, 0, nspin, input_vec_mag, input_angle_mag, ofs_running, noncolin);
+    std::string stdout_output = testing::internal::GetCapturedStdout();
+
+    // All components should be preserved
+    EXPECT_DOUBLE_EQ(atom.m_loc_[0].x, 1.0);
+    EXPECT_DOUBLE_EQ(atom.m_loc_[0].y, 0.5);
+    EXPECT_DOUBLE_EQ(atom.m_loc_[0].z, 0.3);
+
+    // No warning should appear
+    EXPECT_EQ(stdout_output.find("Warning"), std::string::npos);
+}
+
+// ============================================================
+// Tests for issue #5939: nspin=4 autoset should not force (1,1,1)
+// ============================================================
+
+TEST_F(ReadAtomsHelperTest, AutosetMagnetizationNspin4NoAutoset)
+{
+    UnitCell ucell;
+    ucell.ntype = 1;
+    ucell.atoms = new Atom[1];
+    ucell.atoms[0].label = "Fe";
+    ucell.atoms[0].na = 1;
+    ucell.atoms[0].mag.resize(1, 0.0);
+    ucell.atoms[0].m_loc_.resize(1);
+    ucell.atoms[0].m_loc_[0].set(0.0, 0.0, 0.0);
+
+    const int nspin = 4;
+
+    testing::internal::CaptureStdout();
+    unitcell::autoset_magnetization(ucell, nspin, ofs_running);
+    std::string stdout_output = testing::internal::GetCapturedStdout();
+
+    // Magnetization should remain zero (NOT autoset to (1,1,1))
+    EXPECT_DOUBLE_EQ(ucell.atoms[0].m_loc_[0].x, 0.0);
+    EXPECT_DOUBLE_EQ(ucell.atoms[0].m_loc_[0].y, 0.0);
+    EXPECT_DOUBLE_EQ(ucell.atoms[0].m_loc_[0].z, 0.0);
+    EXPECT_DOUBLE_EQ(ucell.atoms[0].mag[0], 0.0);
+
+    // Warning should appear on stdout
+    EXPECT_NE(stdout_output.find("Warning"), std::string::npos);
+    EXPECT_NE(stdout_output.find("no initial magnetization is set in STRU"), std::string::npos);
+
+    // Warning should appear in running log
+    ofs_running.flush();
+    std::ifstream ifs("test_running.log");
+    std::stringstream buffer;
+    buffer << ifs.rdbuf();
+    std::string log_content = buffer.str();
+    ifs.close();
+    EXPECT_NE(log_content.find("Warning"), std::string::npos);
+    EXPECT_NE(log_content.find("no initial magnetization is set in STRU"), std::string::npos);
+
+    delete[] ucell.atoms;
+}
+
+TEST_F(ReadAtomsHelperTest, AutosetMagnetizationNspin4NoWarnWhenMagSet)
+{
+    UnitCell ucell;
+    ucell.ntype = 1;
+    ucell.atoms = new Atom[1];
+    ucell.atoms[0].label = "Fe";
+    ucell.atoms[0].na = 1;
+    ucell.atoms[0].mag.resize(1, 1.5);
+    ucell.atoms[0].m_loc_.resize(1);
+    ucell.atoms[0].m_loc_[0].set(0.0, 0.0, 1.5);
+
+    const int nspin = 4;
+
+    testing::internal::CaptureStdout();
+    unitcell::autoset_magnetization(ucell, nspin, ofs_running);
+    std::string stdout_output = testing::internal::GetCapturedStdout();
+
+    // No autoset should happen
+    EXPECT_DOUBLE_EQ(ucell.atoms[0].mag[0], 1.5);
+    EXPECT_DOUBLE_EQ(ucell.atoms[0].m_loc_[0].z, 1.5);
+
+    // No warning should appear
+    EXPECT_EQ(stdout_output.find("Warning"), std::string::npos);
+
+    delete[] ucell.atoms;
+}
+
+TEST_F(ReadAtomsHelperTest, AutosetMagnetizationNspin2StillAutoset)
+{
+    UnitCell ucell;
+    ucell.ntype = 1;
+    ucell.atoms = new Atom[1];
+    ucell.atoms[0].label = "Fe";
+    ucell.atoms[0].na = 1;
+    ucell.atoms[0].mag.resize(1, 0.0);
+    ucell.atoms[0].m_loc_.resize(1);
+    ucell.atoms[0].m_loc_[0].set(0.0, 0.0, 0.0);
+
+    const int nspin = 2;
+
+    testing::internal::CaptureStdout();
+    unitcell::autoset_magnetization(ucell, nspin, ofs_running);
+    std::string stdout_output = testing::internal::GetCapturedStdout();
+
+    // nspin=2 should still autoset to 1.0 (VASP-compatible behavior)
+    EXPECT_DOUBLE_EQ(ucell.atoms[0].mag[0], 1.0);
+    EXPECT_DOUBLE_EQ(ucell.atoms[0].m_loc_[0].z, 1.0);
+
+    // No warning for nspin=2 autoset
+    EXPECT_EQ(stdout_output.find("Warning"), std::string::npos);
+
+    delete[] ucell.atoms;
+}