Skip to content

Commit 7905e96

Browse files
maki49maki49
committed
fix: handle the colinear case at nspin=4
1 parent a350572 commit 7905e96

2 files changed

Lines changed: 98 additions & 183 deletions

File tree

‎ref/ssg/phase1_plan.md‎

Lines changed: 0 additions & 109 deletions
This file was deleted.

‎source/source_cell/module_symmetry/symm_spin_group.cpp‎

Lines changed: 98 additions & 74 deletions
Original file line numberDiff line numberDiff line change
@@ -12,104 +12,128 @@ void Symmetry::analyze_spin_space_group_nspin4(const Atom* atoms, const Statisti
1212
// else dropped. Arrays are populated EXACTLY like analyze_magnetic_group_nspin4 so all downstream
1313
// consumers (k-reduction, density, EXX) are unchanged; only the spin-rotation source differs.
1414

15-
// GUARD: the independent-spin fit is ill-posed for COLLINEAR moments (rank <= 1). A single moment
16-
// axis does not determine R_spin (any rotation about that axis maps the set), so fit_spin_rotation
17-
// rejects every operation -- INCLUDING identity -- which would collapse the group to nrotk=0 and
18-
// crash downstream. Collinear-in-nspin4 is also physically the nspin=2 regime and should be run as
19-
// such. Detect rank <= 1 (nonmagnetic, or all moments parallel/antiparallel, along ANY axis incl. z)
20-
// and fall back to the magnetic (Shubnikov) subgroup, which is always a correct (if conservative)
21-
// reduction and never involves the fit. spin_space_group_nspin4 stays false -> downstream uses the
22-
// magnetic path (build_wspin / cal_Ms take spin_so3(gmatc), not the fitted arrays).
23-
{
24-
ModuleBase::Vector3<double> axis(0.0, 0.0, 0.0);
25-
bool have_axis = false;
26-
bool rank_ge_2 = false; // found a moment not (anti)parallel to the first -> noncollinear
27-
for (int iat = 0; iat < this->nat; ++iat)
28-
{
29-
const ModuleBase::Vector3<double>& m = atoms[st.iat2it[iat]].m_loc_[st.iat2ia[iat]];
30-
if (m.norm() < this->epsilon) { continue; }
31-
if (!have_axis) { axis = m; have_axis = true; }
32-
else if ((m ^ axis).norm() > 100.0 * this->epsilon) { rank_ge_2 = true; break; } // cross!=0 => not (anti)parallel
33-
}
34-
if (!rank_ge_2)
35-
{
36-
GlobalV::ofs_running << "\n WARNING: symmetry_ssg=1 with nspin=4 requires NONCOLLINEAR moments"
37-
" (rank >= 2); the moment configuration is collinear or nonmagnetic, for"
38-
" which the decoupled-spin fit is ill-posed. Falling back to the magnetic"
39-
" (Shubnikov) subgroup. For full spin-space-group reduction of a COLLINEAR"
40-
" magnet, run nspin=2 (moment along z) -- the collinear SSG path handles it.\n";
41-
this->analyze_magnetic_group_nspin4(atoms, st, latvec);
42-
return;
43-
}
44-
}
45-
4615
const ModuleBase::Matrix3 ilatvec = latvec.Inverse();
47-
std::vector<int> keep;
16+
std::vector<int> keep; // unitary SSG operations
4817
keep.reserve(this->nrotk);
49-
std::vector<ModuleBase::Matrix3> keep_R; // fitted R_spin aligned with keep (-> spin_rotation_ssg)
18+
std::vector<ModuleBase::Matrix3> keep_R; // spin rotation R_spin aligned with keep (-> spin_rotation_ssg)
5019
keep_R.reserve(this->nrotk);
20+
std::vector<int> anti; // antiunitary SSG coset (rank>=2 only)
21+
anti.reserve(this->nrotk);
22+
std::vector<ModuleBase::Matrix3> anti_R; // R_spin aligned with anti (-> spin_rotation_anti_ssg)
23+
anti_R.reserve(this->nrotk);
5124
int nrot_new = 0;
5225

26+
// Moment rank: nonmagnetic (rank 0), collinear (rank 1, all moments (anti)parallel along one axis n_hat),
27+
// or noncollinear (rank >= 2). The independent-spin fit needs rank >= 2; the collinear and
28+
// nonmagnetic cases are handled separately below.
29+
ModuleBase::Vector3<double> axis(0.0, 0.0, 0.0); // n_hat (unnormalized) of the collinear moments
5330
bool has_moment = false;
54-
for (int iat = 0; iat < this->nat && !has_moment; ++iat)
31+
bool rank_ge_2 = false;
32+
for (int iat = 0; iat < this->nat; ++iat)
5533
{
5634
const ModuleBase::Vector3<double>& m = atoms[st.iat2it[iat]].m_loc_[st.iat2ia[iat]];
57-
if (!this->equal(m.x, 0.0) || !this->equal(m.y, 0.0) || !this->equal(m.z, 0.0)) { has_moment = true; }
35+
if (m.norm() < this->epsilon) { continue; }
36+
has_moment = true;
37+
if (axis.norm() < this->epsilon) { axis = m; }
38+
// "^" is the cross product operator of Vector3
39+
else if ((m ^ axis).norm() > 100.0 * this->epsilon) { rank_ge_2 = true; break; } // cross!=0 => not (anti)parallel
5840
}
59-
std::vector<int> anti;
60-
anti.reserve(this->nrotk);
61-
std::vector<ModuleBase::Matrix3> anti_R; // fitted R_spin aligned with anti (-> spin_rotation_anti_ssg)
62-
anti_R.reserve(this->nrotk);
6341

64-
// diagnostic: how many ops the magnetic (spin-locked-to-space) subgroup WOULD keep, for the
65-
// "SSG order > magnetic order" comparison printed below. Pure counting, no state change.
66-
int mag_unitary_count = 0;
42+
if (!has_moment)
43+
{
44+
// rank 0 (nonmagnetic): the fit is meaningless; the magnetic path keeps the full group. Fall back.
45+
this->analyze_magnetic_group_nspin4(atoms, st, latvec);
46+
return;
47+
}
6748

68-
const double fit_tol = 100.0 * this->epsilon; // moment-match tolerance for the fit verification
49+
// diagnostic: how many ops the magnetic (spin-locked-to-space) subgroup WOULD keep as unitary, for
50+
// the "SSG order > magnetic order" comparison printed below. Pure counting, no state change.
51+
int mag_unitary_count = 0;
6952
for (int isym = 0; isym < this->nrotk; ++isym)
7053
{
71-
// gather the moment pairs {m_i} and {m_{g(i)}} over all atoms.
72-
std::vector<ModuleBase::Vector3<double>> from(this->nat), to_plus(this->nat), to_minus(this->nat);
73-
for (int iat = 0; iat < this->nat; ++iat)
54+
const ModuleBase::Matrix3 gmatc = ilatvec * this->gmatrix[isym] * latvec;
55+
const ModuleBase::Matrix3 W = ModuleSymmetry::SpinRotation::spin_so3(gmatc);
56+
bool mag_ok = true;
57+
for (int iat = 0; iat < this->nat && mag_ok; ++iat)
7458
{
7559
const ModuleBase::Vector3<double>& m = atoms[st.iat2it[iat]].m_loc_[st.iat2ia[iat]];
60+
const ModuleBase::Vector3<double> mrot = W * m;
7661
const int jat = this->get_rotated_atom(isym, iat);
7762
const ModuleBase::Vector3<double>& mj = atoms[st.iat2it[jat]].m_loc_[st.iat2ia[jat]];
78-
from[iat] = m;
79-
to_plus[iat] = mj;
80-
to_minus[iat] = mj * (-1.0);
63+
if (!this->equal(mrot.x, mj.x) || !this->equal(mrot.y, mj.y) || !this->equal(mrot.z, mj.z)) { mag_ok = false; }
8164
}
65+
if (mag_ok) { ++mag_unitary_count; }
66+
}
8267

83-
// magnetic-subgroup diagnostic count (spin locked to space).
68+
if (rank_ge_2)
69+
{
70+
// NONCOLLINEAR: fit an independent proper R_spin per operation.
71+
const double fit_tol = 100.0 * this->epsilon; // moment-match tolerance for the fit verification
72+
for (int isym = 0; isym < this->nrotk; ++isym)
8473
{
85-
const ModuleBase::Matrix3 gmatc = ilatvec * this->gmatrix[isym] * latvec;
86-
const ModuleBase::Matrix3 W = ModuleSymmetry::SpinRotation::spin_so3(gmatc);
87-
bool mag_ok = true;
88-
for (int iat = 0; iat < this->nat && mag_ok; ++iat)
74+
std::vector<ModuleBase::Vector3<double>> from(this->nat), to_plus(this->nat), to_minus(this->nat);
75+
for (int iat = 0; iat < this->nat; ++iat)
8976
{
90-
const ModuleBase::Vector3<double> mrot = W * from[iat];
91-
if (!this->equal(mrot.x, to_plus[iat].x) || !this->equal(mrot.y, to_plus[iat].y)
92-
|| !this->equal(mrot.z, to_plus[iat].z)) { mag_ok = false; }
77+
const ModuleBase::Vector3<double>& m = atoms[st.iat2it[iat]].m_loc_[st.iat2ia[iat]];
78+
const int jat = this->get_rotated_atom(isym, iat);
79+
const ModuleBase::Vector3<double>& mj = atoms[st.iat2it[jat]].m_loc_[st.iat2ia[jat]];
80+
from[iat] = m;
81+
to_plus[iat] = mj;
82+
to_minus[iat] = mj * (-1.0);
83+
}
84+
// unitary SSG test: fit a proper R_spin mapping +m.
85+
bool ok = false;
86+
ModuleBase::Matrix3 R = ModuleSymmetry::SpinRotation::fit_spin_rotation(from, to_plus, ok, fit_tol);
87+
if (ok)
88+
{
89+
keep.push_back(isym);
90+
keep_R.push_back(R);
91+
if (isym < this->nrot) { ++nrot_new; }
92+
}
93+
else
94+
{
95+
// antiunitary SSG coset: R_spin maps -m, i.e. Theta*g (Theta = trs/sigma_y downstream).
96+
bool anti_ok = false;
97+
ModuleBase::Matrix3 Ra = ModuleSymmetry::SpinRotation::fit_spin_rotation(from, to_minus, anti_ok, fit_tol);
98+
if (anti_ok) { anti.push_back(isym); anti_R.push_back(Ra); }
9399
}
94-
if (mag_ok) { ++mag_unitary_count; }
95-
}
96-
97-
// unitary SSG test: fit an independent proper R_spin mapping +m.
98-
bool ok = false;
99-
ModuleBase::Matrix3 R = ModuleSymmetry::SpinRotation::fit_spin_rotation(from, to_plus, ok, fit_tol);
100-
if (ok)
101-
{
102-
keep.push_back(isym);
103-
keep_R.push_back(R);
104-
if (isym < this->nrot) { ++nrot_new; }
105100
}
106-
else if (has_moment)
101+
}
102+
else
103+
{
104+
// COLLINEAR (rank 1): the spin part is the discrete Z2 group {I, C2_perp}. For a moment
105+
// configuration +-M n_hat, an operation g that PRESERVES the sign pattern (m_{g(i)}=+m_i) gets
106+
// R_spin=I; one that REVERSES it (m_{g(i)}=-m_i) gets R_spin=C2 about a FIXED axis e_hat perp
107+
// n_hat (a proper 180-deg rotation sending n_hat -> -n_hat, purely UNITARY, no time reversal).
108+
// Using the SAME e_hat for every reversing op makes g->R_spin a valid Z2 homomorphism (C2_e^2=I),
109+
// so the group closes. Unlike the magnetic path (spin_so3 locks spin to space and DISCARDS ops
110+
// whose spatial rotation tilts n_hat off-axis), the decoupled R_spin recovers those ops.
111+
// All kept operations are UNITARY (anti stays empty); downstream is identical to the rank>=2 case.
112+
ModuleBase::Vector3<double> nhat = axis * (1.0 / axis.norm());
113+
const ModuleBase::Vector3<double> ref = (std::fabs(nhat.x) < 0.9)
114+
? ModuleBase::Vector3<double>(1.0, 0.0, 0.0) : ModuleBase::Vector3<double>(0.0, 1.0, 0.0);
115+
ModuleBase::Vector3<double> ehat = ref ^ nhat;
116+
ehat = ehat * (1.0 / ehat.norm()); // fixed unit vector perpendicular to n_hat
117+
// C2 about ehat: R = 2 ehat (x) ehat - I (symmetric; maps ehat->ehat, anything perp ehat -> -it,
118+
// in particular n_hat -> -n_hat since n_hat perp ehat).
119+
const ModuleBase::Matrix3 Rflip(
120+
2.0 * ehat.x * ehat.x - 1.0, 2.0 * ehat.x * ehat.y, 2.0 * ehat.x * ehat.z,
121+
2.0 * ehat.y * ehat.x, 2.0 * ehat.y * ehat.y - 1.0, 2.0 * ehat.y * ehat.z,
122+
2.0 * ehat.z * ehat.x, 2.0 * ehat.z * ehat.y, 2.0 * ehat.z * ehat.z - 1.0);
123+
const ModuleBase::Matrix3 Iden(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0);
124+
for (int isym = 0; isym < this->nrotk; ++isym)
107125
{
108-
// antiunitary SSG coset: R_spin maps -m, i.e. Theta*g (with Theta the trs/sigma_y factor
109-
// applied downstream) plus the fitted proper spin rotation.
110-
bool anti_ok = false;
111-
ModuleBase::Matrix3 Ra = ModuleSymmetry::SpinRotation::fit_spin_rotation(from, to_minus, anti_ok, fit_tol);
112-
if (anti_ok) { anti.push_back(isym); anti_R.push_back(Ra); }
126+
bool preserve = true, reverse = true;
127+
for (int iat = 0; iat < this->nat && (preserve || reverse); ++iat)
128+
{
129+
const ModuleBase::Vector3<double>& m = atoms[st.iat2it[iat]].m_loc_[st.iat2ia[iat]];
130+
const int jat = this->get_rotated_atom(isym, iat);
131+
const ModuleBase::Vector3<double>& mj = atoms[st.iat2it[jat]].m_loc_[st.iat2ia[jat]];
132+
if (!this->equal(mj.x, m.x) || !this->equal(mj.y, m.y) || !this->equal(mj.z, m.z)) { preserve = false; }
133+
if (!this->equal(mj.x, -m.x) || !this->equal(mj.y, -m.y) || !this->equal(mj.z, -m.z)) { reverse = false; }
134+
}
135+
if (preserve) { keep.push_back(isym); keep_R.push_back(Iden); if (isym < this->nrot) { ++nrot_new; } }
136+
else if (reverse) { keep.push_back(isym); keep_R.push_back(Rflip); if (isym < this->nrot) { ++nrot_new; } }
113137
}
114138
}
115139

0 commit comments

Comments
 (0)