Background
The SCPH/QHA code now identifies acoustic modes at the Gamma point from the eigenvectors, by projection onto the mass-weighted rigid-translation subspace (Dynamical::detect_acoustic_modes_at_gamma, ScphQhaCommon::classify_acoustic_modes_from_cmat). This assumes exactly three zero modes at Gamma, which is correct for 3D and 2D periodic systems.
For lower-dimensional systems the count is different:
- 1D (chain/wire along the periodic axis): 3 translations + 1 rotation about the periodic axis = 4 zero modes.
- 0D (molecules/clusters): 3 translations + 3 rotations = 6 zero modes (5 for linear molecules, where the rotation about the molecular axis is trivial).
Current behavior
Neither the previous nor the current implementation handles these cases:
- The previous magnitude-threshold assignment in
zerofill_elements_acoustic_at_gamma exited with "Could not assign acoustic modes at Gamma" as soon as more than three near-zero modes were found.
- The current eigenvector-based detection correctly picks the 3 translations, but a rotational zero mode is then classified as an optical mode: its (numerically noisy) V3/V4 couplings are kept, and in
Scph::compute_qmat_and_dmat the occupation-factor floor for collapsed non-acoustic modes (omega_floor_gamma = smallest harmonic optical frequency) degenerates to eps8 because the rotational mode itself is the smallest "optical" frequency. The resulting near-divergent occupation factor destabilizes the SCPH iteration.
Proposed design
-
Explicit declaration, not auto-detection. Add an input variable (e.g. SYSTEM_DIM = 3/1/0, default 3) declaring the periodicity. Detecting vacuum padding or inferring from KMESH is fragile; an explicit setting also documents intent in the input file. 2D needs no special treatment (in-plane periodicity breaks all rotations; the ZA mode at Gamma is the rigid out-of-plane translation).
-
Generalize the detector to the rigid-body subspace. In addition to the three mass-weighted translations, build rotation generators sqrt(m_j) * (e_alpha x (r_j - r_cm)) — only the periodic-axis rotation for 1D, all three for 0D — orthonormalized against the translations with a rank tolerance. The rank reduction automatically yields 5 instead of 6 zero modes for linear molecules. Downstream code is already flag-driven, so the changes are localized:
zerofill_elements_acoustic_at_gamma: count check == 3 -> == N;setup_structural_opt_buffers and the relaxation code: ns - 3 -> ns - N (needs an audit of the structural-optimization loops);classify_acoustic_modes_from_cmat: works unchanged (overlap against the flagged harmonic subspace);omega_floor_gamma: automatically excludes rotations once they are flagged.
-
Physical prerequisite. Zeroing V3/V4 elements of rotational modes is only legitimate when the force constants satisfy the Born-Huang rotational sum rules (ALM's ICONST options). The projection provides a free diagnostic: a mode with high rotational projection but a harmonic frequency well away from zero indicates broken rotational invariance, and the code should warn or stop rather than proceed. Excluding rigid rotations from the thermal average (Q = 0, like translations) mirrors the standard treatment of the six rigid degrees of freedom in molecular VSCF; rovibrational coupling remains beyond SCPH's scope.
-
Phased implementation.
- Phase (a), cheap safety: after detecting the translations, compute the rotational projections of any remaining near-zero Gamma modes and stop with a clear message ("N zero modes detected at Gamma; 1D/0D systems with rotational zero modes are not yet supported") instead of the current obscure misbehavior.
- Phase (b), full support: implement
SYSTEM_DIM together with a regression fixture (a diatomic molecule or a carbyne-like chain).
Related
tools/GenDisplacement.py (_identify_gamma_acoustic) has the same blind spot: the displacement-uniformity score identifies only translations, so for a cluster the random-displacement mode would thermally excite rigid rotations via the omega_min floor. The projection method generalizes to rotations; the uniformity test does not — porting the projection approach to the Python side would let both codes share one definition of the rigid-body subspace.
Background
The SCPH/QHA code now identifies acoustic modes at the Gamma point from the eigenvectors, by projection onto the mass-weighted rigid-translation subspace (
Dynamical::detect_acoustic_modes_at_gamma,ScphQhaCommon::classify_acoustic_modes_from_cmat). This assumes exactly three zero modes at Gamma, which is correct for 3D and 2D periodic systems.For lower-dimensional systems the count is different:
Current behavior
Neither the previous nor the current implementation handles these cases:
zerofill_elements_acoustic_at_gammaexited with "Could not assign acoustic modes at Gamma" as soon as more than three near-zero modes were found.Scph::compute_qmat_and_dmatthe occupation-factor floor for collapsed non-acoustic modes (omega_floor_gamma= smallest harmonic optical frequency) degenerates to eps8 because the rotational mode itself is the smallest "optical" frequency. The resulting near-divergent occupation factor destabilizes the SCPH iteration.Proposed design
Explicit declaration, not auto-detection. Add an input variable (e.g.
SYSTEM_DIM = 3/1/0, default 3) declaring the periodicity. Detecting vacuum padding or inferring from KMESH is fragile; an explicit setting also documents intent in the input file. 2D needs no special treatment (in-plane periodicity breaks all rotations; the ZA mode at Gamma is the rigid out-of-plane translation).Generalize the detector to the rigid-body subspace. In addition to the three mass-weighted translations, build rotation generators
sqrt(m_j) * (e_alpha x (r_j - r_cm))— only the periodic-axis rotation for 1D, all three for 0D — orthonormalized against the translations with a rank tolerance. The rank reduction automatically yields 5 instead of 6 zero modes for linear molecules. Downstream code is already flag-driven, so the changes are localized:zerofill_elements_acoustic_at_gamma: count check== 3->== N;setup_structural_opt_buffersand the relaxation code:ns - 3->ns - N(needs an audit of the structural-optimization loops);classify_acoustic_modes_from_cmat: works unchanged (overlap against the flagged harmonic subspace);omega_floor_gamma: automatically excludes rotations once they are flagged.Physical prerequisite. Zeroing V3/V4 elements of rotational modes is only legitimate when the force constants satisfy the Born-Huang rotational sum rules (ALM's ICONST options). The projection provides a free diagnostic: a mode with high rotational projection but a harmonic frequency well away from zero indicates broken rotational invariance, and the code should warn or stop rather than proceed. Excluding rigid rotations from the thermal average (Q = 0, like translations) mirrors the standard treatment of the six rigid degrees of freedom in molecular VSCF; rovibrational coupling remains beyond SCPH's scope.
Phased implementation.
SYSTEM_DIMtogether with a regression fixture (a diatomic molecule or a carbyne-like chain).Related
tools/GenDisplacement.py(_identify_gamma_acoustic) has the same blind spot: the displacement-uniformity score identifies only translations, so for a cluster the random-displacement mode would thermally excite rigid rotations via the omega_min floor. The projection method generalizes to rotations; the uniformity test does not — porting the projection approach to the Python side would let both codes share one definition of the rigid-body subspace.