Feature: fully utilized magnetic space group (MSG) enabling SOC (#7692)

* SOC symmetry: magmom filter, charge&DM symmetrization

* Hexx(R) rotation

* Feat: magnetic (Shubnikov) group

* remove TRS_first_

* Fix rhog_symmetry: anti-unitary ops (used for reducing k-points) should be included

* Perf: antiunitary ops are also used to reduce Hexx(R)

* do not autoset magmom at symmetry=1

* Fix SOC bug: θ=π branch of so3_to_su2, and M⊗U convention

* fix the phase factor: compatible with the formula

* Fix: enable mixing_restart for EXX (instead of segfault)

* print magnetic space group

* doc and testing

* refactor: parameterize symmetry in module_cell

* refactor: remove default values of new parameters

* fix parameter description out-of-sync

* fix: if symmetry is not set, default -1 at lspinorb=1 (as before)

* Fix: explicitly set mag 1 for nspin-2+symmetry test cases (was autoset to 1 before)

* fix a UT format

* refine the warning lines skipping autoset magmom

* update results of nspin4+sym1 cases (wrong before spin-mixed charge density symmetrization is implemented)

Author: maki49

docs/advanced/input_files/input-main.md

@@ -614,7 +614,7 @@
 - **Description**: Takes value 1, 0 or -1.
   - -1: No symmetry will be considered. It is recommended to set -1 for non-colinear + soc calculations, where time reversal symmetry is broken sometimes.
   - 0: Only time reversal symmetry would be considered in symmetry operations, which implied k point and -k point would be treated as a single k point with twice the weight.
-  - 1: Symmetry analysis will be performed to determine the type of Bravais lattice and associated symmetry operations. (point groups, space groups, primitive cells, and irreducible k-points)
+  - 1: Symmetry analysis will be performed to determine the type of Bravais lattice and associated symmetry operations (point groups, space groups, primitive cells, and irreducible k-points). For a magnetic system, the symmetry of the initial magnetic structure will be analyzed and preserved.
 
   > Note: When symmetry is enabled (value 1), k-points are reduced to the irreducible Brillouin zone (IBZ). For explicit k-point lists with custom weights (see KPT file), the custom weights are preserved during symmetry reduction. For Monkhorst-Pack grids, uniform weights are used.
 - **Default**: default

docs/parameters.yaml

@@ -63,7 +63,7 @@ parameters:
       Takes value 1, 0 or -1.
       * -1: No symmetry will be considered. It is recommended to set -1 for non-colinear + soc calculations, where time reversal symmetry is broken sometimes.
       * 0: Only time reversal symmetry would be considered in symmetry operations, which implied k point and -k point would be treated as a single k point with twice the weight.
-      * 1: Symmetry analysis will be performed to determine the type of Bravais lattice and associated symmetry operations. (point groups, space groups, primitive cells, and irreducible k-points)
+      * 1: Symmetry analysis will be performed to determine the type of Bravais lattice and associated symmetry operations (point groups, space groups, primitive cells, and irreducible k-points). For a magnetic system, the symmetry of the initial magnetic structure will be analyzed and preserved.
 
       [NOTE] When symmetry is enabled (value 1), k-points are reduced to the irreducible Brillouin zone (IBZ). For explicit k-point lists with custom weights (see KPT file), the custom weights are preserved during symmetry reduction. For Monkhorst-Pack grids, uniform weights are used.
     default_value: default

python/pyabacus/src/ModuleDriver/py_driver.cpp

@@ -428,7 +428,8 @@ CalculationResult PyDriver::run(
         impl_->ucell_->setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running, PARAM.inp.symmetry_prec, PARAM.inp.dfthalf_type, PARAM.inp.pseudo_dir, PARAM.inp.nspin,
             PARAM.inp.basis_type, PARAM.inp.orbital_dir, PARAM.inp.init_wfc,
             PARAM.inp.onsite_radius, PARAM.globalv.deepks_setorb, PARAM.inp.rpa,
-            PARAM.inp.fixed_atoms, PARAM.inp.noncolin, PARAM.inp.calculation, PARAM.inp.esolver_type);
+            PARAM.inp.fixed_atoms, PARAM.inp.noncolin, PARAM.inp.calculation, PARAM.inp.esolver_type,
+            std::stoi(PARAM.inp.symmetry));
 
         // Check atomic structure
         unitcell::check_atomic_stru(*impl_->ucell_, PARAM.inp.min_dist_coef);

source/source_cell/k_vector_utils.cpp

@@ -539,7 +539,23 @@ void kvec_ibz_kpoint(K_Vectors& kv,
             kgmatrix[i] = symm.kgmatrix[i];
         }
 
-        if (!include_inv)
+        if (symm.magnetic_nspin4)
+        {
+            // (nspin=4, magnetic) Time reversal Theta reverses the magnetization, so Theta alone is
+            // NOT a symmetry and the blanket "-k is always equivalent" doubling below is invalid.
+            // Only the antiunitary elements Theta*g with g in the moment-reversing coset belong to
+            // the Shubnikov group; append exactly those, keeping the index convention
+            // j + nrotk  <->  Theta * gmatrix_anti[j]  (decoded the same way in restore_dm).
+            // (nspin=2 is unaffected: there the antiunitary operation is plain conjugation K, which
+            //  does not touch the spin, so D_s(-k)=D_s^*(k) holds even for a ferromagnet and the
+            //  generic branch below stays correct.)
+            for (int j = 0; j < symm.nrotk_anti; ++j)
+            {
+                kgmatrix[j + symm.nrotk] = inv * symm.kgmatrix_anti[j];
+            }
+            nrotkm = symm.nrotk + symm.nrotk_anti;
+        }
+        else if (!include_inv)
         {
             for (int i = 0; i < symm.nrotk; ++i)
             {

source/source_cell/module_symmetry/CMakeLists.txt

@@ -12,6 +12,7 @@ add_library(
     symm_pricell.cpp
     symm_rho.cpp
     symmetry.cpp
+    symmetry_rotation_spin.cpp
 )
 
 if(ENABLE_COVERAGE)

source/source_cell/module_symmetry/symm_analysis.cpp

@@ -290,6 +290,13 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms,
 
     this->set_atom_map(atoms); // find the atom mapping according to the symmetry operations
 
+    // (nspin=4 / SOC) restrict to the unitary magnetic subgroup: drop operations that reverse
+    // the magnetization (pseudovector), so they are not applied in k-reduction / density symmetrization.
+    if (nspin == 4)
+    {
+        this->analyze_magnetic_group_nspin4(atoms, st, latvec1);
+    }
+
     // Do this here for debug
     if (calculation == "relax")
     {

source/source_cell/module_symmetry/symm_magnetic.cpp

@@ -1,7 +1,11 @@
 #include "symmetry.h"
 using namespace ModuleSymmetry;
 
+#include "symmetry_rotation_spin.h"
+#include "source_io/module_parameter/parameter.h"
+
 #include <set>
+#include <vector>
 
 void Symmetry::analyze_magnetic_group(const Atom* atoms, const Statistics& st, int& nrot_out, int& nrotk_out)
 {
@@ -74,6 +78,116 @@ void Symmetry::analyze_magnetic_group(const Atom* atoms, const Statistics& st, i
 
 }
 
+void Symmetry::analyze_magnetic_group_nspin4(const Atom* atoms, const Statistics& st, const ModuleBase::Matrix3& latvec)
+{
+    // Restrict the space group to the unitary magnetic subgroup (nspin=4 / SOC):
+    // operation g survives if it preserves the magnetic configuration as a pseudovector,
+    // i.e. W(g) m_i = m_{g(i)} for every atom, with W(g) = spin_so3(gmatc). 
+    // Operations that reverse the moment (only symmetries together with time reversal) are dropped, 
+    // so they are no longer applied in k-reduction or density symmetrization. 
+    // Non-magnetic (m_i=0) keeps every operation.
+    const ModuleBase::Matrix3 ilatvec = latvec.Inverse();
+    std::vector<int> keep;
+    keep.reserve(this->nrotk);
+    int nrot_new = 0;
+
+    // Is the configuration actually magnetic? For m_i = 0 every operation both "preserves" and
+    // "reverses" the moment, so the antiunitary coset is meaningless there: 
+    // Theta (TRS) itself is a symmetry and the grey group is handled by the usual -k shortcut in the k-reduction).
+    bool has_moment = false;
+    for (int iat = 0; iat < this->nat && !has_moment; ++iat)
+    {
+        const ModuleBase::Vector3<double>& m = atoms[st.iat2it[iat]].m_loc_[st.iat2ia[iat]];
+        if (!this->equal(m.x, 0.0) || !this->equal(m.y, 0.0) || !this->equal(m.z, 0.0)) { has_moment = true; }
+    }
+    std::vector<int> anti;   // operations that REVERSE the moment: Theta*g is a symmetry
+    anti.reserve(this->nrotk);
+    for (int isym = 0; isym < this->nrotk; ++isym)
+    {
+        const ModuleBase::Matrix3 gmatc = ilatvec * this->gmatrix[isym] * latvec;
+        const ModuleBase::Matrix3 W = ModuleSymmetry::SpinRotation::spin_so3(gmatc);
+        bool ok = true;
+        for (int iat = 0; iat < this->nat && ok; ++iat)
+        {
+            const ModuleBase::Vector3<double>& m = atoms[st.iat2it[iat]].m_loc_[st.iat2ia[iat]];
+            // pseudovector-rotated moment W*m (column-vector convention: m'^i = W_ij m^j)
+            const ModuleBase::Vector3<double>& mrot = W * m;
+            const int jat = this->get_rotated_atom(isym, iat);
+            const ModuleBase::Vector3<double>& mj = atoms[st.iat2it[jat]].m_loc_[st.iat2ia[jat]];
+            if (!this->equal(mrot.x, mj.x) || !this->equal(mrot.y, mj.y) || !this->equal(mrot.z, mj.z)) { ok = false; }
+        }
+        if (ok)
+        {
+            keep.push_back(isym);
+            if (isym < this->nrot) { ++nrot_new; }   // pure point-group rotations are the first nrot ops
+        }
+        else if (has_moment)
+        {
+            // g does not preserve m; check whether it exactly REVERSES it, i.e.
+            // W(g) m_i = -m_{g(i)} for every atom. Then g alone is not a symmetry but the
+            // antiunitary element Theta*g is, and it belongs to the Shubnikov group.
+            bool anti_ok = true;
+            for (int iat = 0; iat < this->nat && anti_ok; ++iat)
+            {
+                const ModuleBase::Vector3<double>& m = atoms[st.iat2it[iat]].m_loc_[st.iat2ia[iat]];
+                const ModuleBase::Vector3<double>& mrot = W * m;
+                const int jat = this->get_rotated_atom(isym, iat);
+                const ModuleBase::Vector3<double>& mj = atoms[st.iat2it[jat]].m_loc_[st.iat2ia[jat]];
+                if (!this->equal(mrot.x, -mj.x) || !this->equal(mrot.y, -mj.y) || !this->equal(mrot.z, -mj.z)) { anti_ok = false; }
+            }
+            if (anti_ok) { anti.push_back(isym); }
+        }
+    }
+
+    // Capture the antiunitary coset BEFORE the unitary arrays are compacted in place below
+    // (the compaction overwrites gmatrix/kgmatrix/gtrans/isym_rotiat_ and would lose them).
+    this->magnetic_nspin4 = has_moment;
+    this->nrotk_anti = static_cast<int>(anti.size());
+    if (this->nrotk_anti > 0)
+    {
+        this->isym_rotiat_anti_.resize(this->nrotk_anti);
+        for (int j = 0; j < this->nrotk_anti; ++j)
+        {
+            const int isym = anti[j];
+            this->gmatrix_anti[j] = this->gmatrix[isym];
+            this->kgmatrix_anti[j] = this->kgmatrix[isym];
+            this->gtrans_anti[j] = this->gtrans[isym];
+            this->isym_rotiat_anti_[j] = this->isym_rotiat_[isym];
+        }
+        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,
+            "MAGNETIC ANTIUNITARY OPERATIONS (Theta*g)", this->nrotk_anti);
+    }
+
+    const int nrotk_new = static_cast<int>(keep.size());
+    if (nrotk_new == this->nrotk) { return; }   // nothing removed (non-magnetic or fully-preserving group)
+
+    // compact the operation arrays in ascending order (keeps the rotations-first layout).
+    for (int i = 0; i < nrotk_new; ++i)
+    {
+        const int isym = keep[i];
+        if (i != isym)
+        {
+            this->gmatrix[i] = this->gmatrix[isym];
+            this->kgmatrix[i] = this->kgmatrix[isym];
+            this->gtrans[i] = this->gtrans[isym];
+            this->isym_rotiat_[i] = this->isym_rotiat_[isym];
+        }
+    }
+    this->isym_rotiat_.resize(nrotk_new);
+    this->nrot = nrot_new;
+    this->nrotk = nrotk_new;
+
+    // refresh the point-/space-group labels for the reduced (unitary magnetic) group
+    this->pointgroup(this->nrot, this->pgnumber, this->pgname, this->gmatrix, GlobalV::ofs_running, nullptr);
+    this->pointgroup(this->nrotk, this->spgnumber, this->spgname, this->gmatrix, GlobalV::ofs_running, nullptr);
+    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "MAGNETIC POINT GROUP (unitary, nspin=4)", this->pgname);
+    // space-group-consistent name of the unitary magnetic group (from nrotk, matching "POINT GROUP IN
+    // SPACE GROUP"); pgname above is the pure-point-group-block name, which under-detects for hexagonal
+    // (e.g. Co prints S_6 there but is C_6h here).
+    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "MAGNETIC POINT GROUP IN SPACE GROUP", this->spgname);
+    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "MAGNETIC SPACE GROUP OPERATIONS", this->nrotk);
+}
+
 bool Symmetry::magmom_same_check(const Atom* atoms)const
 {
     ModuleBase::TITLE("Symmetry", "magmom_same_check");
@@ -96,3 +210,32 @@ bool Symmetry::magmom_same_check(const Atom* atoms)const
     return pricell_loop;
 }
 
+int Symmetry::density_sym_ops(std::vector<ModuleBase::Matrix3>& kgmat,
+                              std::vector<ModuleBase::Vector3<double>>& gtr,
+                              std::vector<double>& trs_inv) const
+{
+    // The density must be symmetrized with the SAME group that was used to fold the k-points
+    // (see KVectorUtils::ibz_kpoint): otherwise the density accumulated over the IBZ is not
+    // restored to the full BZ result. For nspin=4 with a non-zero moment that group is the
+    // Shubnikov group H + Theta*A, so the antiunitary elements' spatial parts are appended here.
+    // Theta leaves the charge invariant and reverses the magnetization, which is what `trs_inv`
+    // encodes; the spatial bookkeeping (orbit grouping, phases) is identical for both kinds.
+    const int nu = this->nrotk;
+    const int na = (this->magnetic_nspin4 ? this->nrotk_anti : 0);
+    kgmat.resize(nu + na);
+    gtr.resize(nu + na);
+    trs_inv.assign(nu + na, 1.0);
+    for (int i = 0; i < nu; ++i)
+    {
+        kgmat[i] = this->kgmatrix[i];
+        gtr[i] = this->gtrans[i];
+    }
+    for (int j = 0; j < na; ++j)
+    {
+        kgmat[nu + j] = this->kgmatrix_anti[j];
+        gtr[nu + j] = this->gtrans_anti[j];
+        trs_inv[nu + j] = -1.0;
+    }
+    return nu + na;
+}
+