feat: SSG for nspin=4: Rm(i)=m(gi) fitting

Author: maki49

source/Makefile.Objects

@@ -523,6 +523,7 @@ OBJS_SYMMETRY=symm_other.o\
     symm_hermite.o\
     symm_lattice.o\
     symm_magnetic.o\
+    symm_spin_group.o\
     symm_pricell.o\
     symm_rho.o\
     symmetry.o\

source/source_cell/module_symmetry/CMakeLists.txt

@@ -9,6 +9,7 @@ add_library(
     symm_hermite.cpp
     symm_lattice.cpp
     symm_magnetic.cpp
+    symm_spin_group.cpp
     symm_pricell.cpp
     symm_rho.cpp
     symmetry.cpp

source/source_cell/module_symmetry/symm_analysis.cpp

@@ -304,7 +304,17 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms,
     // the magnetization (pseudovector), so they are not applied in k-reduction / density symmetrization.
     if (nspin == 4)
     {
-        this->analyze_magnetic_group_nspin4(atoms, st, latvec1);
+        // (SSG, SOC off) decouple the spin rotation from space: fit an independent R_spin per op,
+        // keeping more operations than the magnetic (Shubnikov) subgroup. Falls back to the magnetic
+        // path when SOC is on (spin locked to space) or symmetry_ssg is disabled -> bit-for-bit unchanged.
+        if (PARAM.inp.symmetry_ssg && !PARAM.inp.lspinorb)
+        {
+            this->analyze_spin_space_group_nspin4(atoms, st, latvec1);
+        }
+        else
+        {
+            this->analyze_magnetic_group_nspin4(atoms, st, latvec1);
+        }
     }
 
     // (nspin=2 collinear spin space group) split the full chemical group into the unitary magnetic

source/source_cell/module_symmetry/symm_magnetic.cpp

@@ -188,101 +188,6 @@ void Symmetry::analyze_magnetic_group_nspin4(const Atom* atoms, const Statistics
     ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "MAGNETIC SPACE GROUP OPERATIONS", this->nrotk);
 }
 
-void Symmetry::analyze_spin_space_group_nspin2(const Atom* atoms, const Statistics& st)
-{
-    // (nspin=2 collinear spin space group) The FULL chemical space group is already in
-    // gmatrix[0..nrotk). Split it, using the scalar collinear moments mag[iat], into:
-    //   - unitary magnetic subgroup : mag[iat] == mag[g(iat)]  for every atom  -> kept in gmatrix
-    //   - unitary spin-flip coset   : mag[iat] == -mag[g(iat)] for every atom  -> gmatrix_flip[]
-    // The spin-flip coset elements are [C2_perp||g]: they swap the up/down spin channels while
-    // rotating space by g. There is NO time reversal here (the collinear Hamiltonian is real, the two
-    // spin blocks are independent), so this is a plain channel swap - distinct from the antiunitary
-    // Theta*g coset built by analyze_magnetic_group_nspin4. Operations that neither preserve nor exactly
-    // flip the moment pattern are dropped.
-    std::vector<int> keep;   // unitary magnetic subgroup
-    keep.reserve(this->nrotk);
-    std::vector<int> flip;   // unitary spin-flip coset [C2_perp||g]
-    flip.reserve(this->nrotk);
-    int nrot_new = 0;
-    for (int isym = 0; isym < this->nrotk; ++isym)
-    {
-        bool preserve = true;
-        for (int iat = 0; iat < this->nat && preserve; ++iat)
-        {
-            const double mi = atoms[st.iat2it[iat]].mag[st.iat2ia[iat]];
-            const int jat = this->get_rotated_atom(isym, iat);
-            const double mj = atoms[st.iat2it[jat]].mag[st.iat2ia[jat]];
-            if (!this->equal(mi, mj)) { preserve = false; }
-        }
-        if (preserve)
-        {
-            keep.push_back(isym);
-            if (isym < this->nrot) { ++nrot_new; }   // pure point-group rotations are the first nrot ops
-            continue;
-        }
-        // g does not preserve the moment pattern; check whether it exactly FLIPS it, i.e.
-        // mag[iat] = -mag[g(iat)] for every atom. Then [C2_perp||g] (spatial g + up<->down swap) is a symmetry.
-        bool flip_ok = true;
-        for (int iat = 0; iat < this->nat && flip_ok; ++iat)
-        {
-            const double mi = atoms[st.iat2it[iat]].mag[st.iat2ia[iat]];
-            const int jat = this->get_rotated_atom(isym, iat);
-            const double mj = atoms[st.iat2it[jat]].mag[st.iat2ia[jat]];
-            if (!this->equal(mi, -mj)) { flip_ok = false; }
-        }
-        if (flip_ok) { flip.push_back(isym); }
-    }
-
-    // Capture the spin-flip coset BEFORE the unitary arrays are compacted in place below.
-    this->nrotk_flip = static_cast<int>(flip.size());
-    this->spin_flip_nspin2 = (this->nrotk_flip > 0);
-    if (this->nrotk_flip > 0)
-    {
-        this->isym_rotiat_flip_.resize(this->nrotk_flip);
-        for (int j = 0; j < this->nrotk_flip; ++j)
-        {
-            const int isym = flip[j];
-            this->gmatrix_flip[j] = this->gmatrix[isym];
-            this->kgmatrix_flip[j] = this->kgmatrix[isym];
-            this->gtrans_flip[j] = this->gtrans[isym];
-            this->isym_rotiat_flip_[j] = this->isym_rotiat_[isym];
-        }
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,
-            "SPIN-FLIP COSET OPERATIONS (nspin=2 SSG)", this->nrotk_flip);
-    }
-
-    const int nrotk_new = static_cast<int>(keep.size());
-    if (nrotk_new != this->nrotk)
-    {
-        // compact the unitary subgroup 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=2)", this->pgname);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "MAGNETIC SPACE GROUP OPERATIONS", this->nrotk);
-    }
-    if (this->spin_flip_nspin2)
-    {
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,
-            "SPIN SPACE GROUP OPERATIONS (unitary + spin-flip)", this->nrotk + this->nrotk_flip);
-    }
-}
-
 bool Symmetry::magmom_same_check(const Atom* atoms)const
 {
     ModuleBase::TITLE("Symmetry", "magmom_same_check");
@@ -334,32 +239,3 @@ int Symmetry::density_sym_ops(std::vector<ModuleBase::Matrix3>& kgmat,
     return nu + na;
 }
 
-int Symmetry::spin_flip_sym_ops(std::vector<ModuleBase::Matrix3>& kgmat,
-                                std::vector<ModuleBase::Vector3<double>>& gtr,
-                                std::vector<double>& flip_sign) const
-{
-    // (nspin=2 collinear SSG) Assemble the full spin space group used to symmetrize the collinear
-    // density and to fold the k-points: the nrotk unitary operations (flip_sign +1) followed by the
-    // nrotk_flip spatial parts of the spin-flip coset [C2_perp||g] (flip_sign -1). The combined set
-    // is the full chemical space group -- a group, closed under inverse -- so the invmap/grouping in
-    // rhog_symmetry* remains valid. In the (charge, mag) basis the charge is invariant (all ops act
-    // as ordinary space-group operations) and the magnetization flips sign under the coset.
-    const int nu = this->nrotk;
-    const int nf = (this->spin_flip_nspin2 ? this->nrotk_flip : 0);
-    kgmat.resize(nu + nf);
-    gtr.resize(nu + nf);
-    flip_sign.assign(nu + nf, 1.0);
-    for (int i = 0; i < nu; ++i)
-    {
-        kgmat[i] = this->kgmatrix[i];
-        gtr[i] = this->gtrans[i];
-    }
-    for (int j = 0; j < nf; ++j)
-    {
-        kgmat[nu + j] = this->kgmatrix_flip[j];
-        gtr[nu + j] = this->gtrans_flip[j];
-        flip_sign[nu + j] = -1.0;
-    }
-    return nu + nf;
-}
-

source/source_cell/module_symmetry/symm_rot_spin.cpp

@@ -125,6 +125,60 @@ ModuleBase::Matrix3 spin_so3(const ModuleBase::Matrix3& gmatc)
     return proper_part(gmatc).Transpose();
 }
 
+ModuleBase::Matrix3 fit_spin_rotation(const std::vector<ModuleBase::Vector3<double>>& from,
+                                      const std::vector<ModuleBase::Vector3<double>>& to,
+                                      bool& ok,
+                                      const double tol)
+{
+    ok = false;
+    const ModuleBase::Matrix3 identity(1, 0, 0, 0, 1, 0, 0, 0, 1);
+    const int n = static_cast<int>(from.size());
+    if (n == 0 || static_cast<int>(to.size()) != n) { return identity; }
+
+    // A norm below this is treated as a zero (non-magnetic) moment: it constrains nothing.
+    const double nrm_eps = 1e-6;
+
+    // 1st frame axis: the first non-zero moment.
+    int a = -1;
+    for (int i = 0; i < n; ++i) { if (from[i].norm() > nrm_eps) { a = i; break; } }
+    if (a < 0) { return identity; }   // all moments zero: no constraint (handled by has_moment gate upstream)
+    const ModuleBase::Vector3<double> e1 = from[a] * (1.0 / from[a].norm());
+    const ModuleBase::Vector3<double> f1 = to[a] * (1.0 / std::max(to[a].norm(), nrm_eps));
+
+    // 2nd frame axis: the first moment with a component perpendicular to e1 (rank >= 2).
+    int b = -1;
+    ModuleBase::Vector3<double> e2, f2;
+    for (int i = 0; i < n; ++i)
+    {
+        const ModuleBase::Vector3<double> perp = from[i] - e1 * (e1 * from[i]);
+        if (perp.norm() > nrm_eps) { b = i; e2 = perp * (1.0 / perp.norm()); break; }
+    }
+    if (b < 0) { return identity; }   // rank-1 (collinear moments): degenerate, out of nspin=4 SSG scope
+    const ModuleBase::Vector3<double> fperp = to[b] - f1 * (f1 * to[b]);
+    if (fperp.norm() <= nrm_eps) { return identity; }   // target lost its perpendicular part -> not a rotation
+    f2 = fperp * (1.0 / fperp.norm());
+
+    // 3rd axis by right-handed cross product -> both frames right-handed -> det(R)=+1.
+    const ModuleBase::Vector3<double> e3 = e1 ^ e2;
+    const ModuleBase::Vector3<double> f3 = f1 ^ f2;
+
+    // R maps e_k -> f_k :  R = F E^T,  R_ij = sum_k f_k[i] e_k[j].  Acts as column vectors: R * m.
+    const ModuleBase::Matrix3 R(
+        f1.x * e1.x + f2.x * e2.x + f3.x * e3.x, f1.x * e1.y + f2.x * e2.y + f3.x * e3.y, f1.x * e1.z + f2.x * e2.z + f3.x * e3.z,
+        f1.y * e1.x + f2.y * e2.x + f3.y * e3.x, f1.y * e1.y + f2.y * e2.y + f3.y * e3.y, f1.y * e1.z + f2.y * e2.z + f3.y * e3.z,
+        f1.z * e1.x + f2.z * e2.x + f3.z * e3.x, f1.z * e1.y + f2.z * e2.y + f3.z * e3.y, f1.z * e1.z + f2.z * e2.z + f3.z * e3.z);
+
+    // Verify R exactly maps EVERY moment (this is what makes the fit correct AND, for rank>=2
+    // moments, guarantees group closure: R is the unique proper map of the moment set).
+    for (int i = 0; i < n; ++i)
+    {
+        const ModuleBase::Vector3<double> d = R * from[i] - to[i];
+        if (std::fabs(d.x) > tol || std::fabs(d.y) > tol || std::fabs(d.z) > tol) { return identity; }
+    }
+    ok = true;
+    return R;
+}
+
 ModuleBase::Matrix3 pauli_rotation_matrix(const Su2& U)
 {
     // sigma matrices (row-major 2x2)

source/source_cell/module_symmetry/symm_rot_spin.h

@@ -52,6 +52,29 @@ Su2 so3_to_su2(const ModuleBase::Matrix3& gmatc, const double eps = 1e-6);
 ///        Computed directly from gmatc as W = R_proper^T (R_proper = proper part).
 ModuleBase::Matrix3 spin_so3(const ModuleBase::Matrix3& gmatc);
 
+/// @brief (nspin=4 SSG) Fit the INDEPENDENT proper spin rotation R (det=+1) that maps one
+///        set of moment vectors onto another: R * from[i] == to[i] for all i, using the
+///        same column-vector action `R * m` as spin_so3() (so it drops in for W).
+///
+/// Unlike spin_so3(gmatc), R is NOT derived from any spatial operation -- it is found from
+/// the moment configuration alone, which is exactly the spin/space decoupling of the spin
+/// space group (valid with SOC off). Built by an orthonormal-frame (Gram-Schmidt) construction
+/// from two linearly independent moment pairs, so no SVD is needed and det(R)=+1 is guaranteed:
+///   e1=norm(from[a]), e2=norm(from[b] perp e1), e3=e1 x e2;  f_k likewise from to[];  R = F E^T.
+/// Then R is VERIFIED to map every moment within `tol`; `ok` is set false if the verification
+/// fails or the moments are rank-deficient (collinear -> no independent second axis; that is the
+/// nspin=2 collinear regime, out of Phase-2 scope). A false `ok` means "reject this operation".
+///
+/// @param from  source moments m_i (all atoms)
+/// @param to    target moments (m_{g(i)} for the unitary test, -m_{g(i)} for the antiunitary test)
+/// @param[out] ok  true iff a proper rotation exactly mapping all moments (within tol) was found
+/// @param tol   absolute tolerance for the per-component moment match (default matches Symmetry::epsilon usage)
+/// @return the fitted rotation R (identity when ok is false; caller must check ok)
+ModuleBase::Matrix3 fit_spin_rotation(const std::vector<ModuleBase::Vector3<double>>& from,
+                                      const std::vector<ModuleBase::Vector3<double>>& to,
+                                      bool& ok,
+                                      const double tol = 1e-5);
+
 /// @brief The same W matrix computed independently from a given SU(2) matrix U via
 ///        W_ij = (1/2) Tr(sigma_i U sigma_j U^dagger). Used for verification.
 ModuleBase::Matrix3 pauli_rotation_matrix(const Su2& U);