Feature: backport SOC magnetic symmetry and extend LibRPA exports

docs/advanced/input_files/input-main.md

@@ -182,6 +182,8 @@
     - [out\_element\_info](#out_element_info)
     - [restart\_save](#restart_save)
     - [rpa](#rpa)
+    - [out\_librpa\_reader\_version](#out_librpa_reader_version)
+    - [out\_librpa\_abf\_overlap](#out_librpa_abf_overlap)
     - [out\_pchg](#out_pchg)
     - [out\_wfc\_norm](#out_wfc_norm)
     - [out\_wfc\_re\_im](#out_wfc_re_im)
@@ -595,7 +597,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
@@ -2076,6 +2078,20 @@
   > Note: If symmetry is set to 1, additional files containing the necessary information for exploiting symmetry in the subsequent rpa calculation will be output: irreducible_sector.txt, symrot_k.txt and symrot_R.txt.
 - **Default**: False
 
+### out_librpa_reader_version
+
+- **Type**: Integer
+- **Availability**: *Numerical atomic orbital basis with rpa=True.*
+- **Description**: Select the ABACUS output format for files consumed by LibRPA. 0 writes the legacy text files, and 1 writes LibRPA reader-v1 binary files directly.
+- **Default**: 0
+
+### out_librpa_abf_overlap
+
+- **Type**: Boolean
+- **Availability**: *Numerical atomic orbital basis with rpa=True, reader version 1, and shrink ABFs.*
+- **Description**: Writes dense raw active-ABF q-space overlap matrices as `v1_abf_overlap_active_iq_<iq>.dat` for offline PSD diagnostics; this can be high-cost in memory and output size. The writer fails closed if an actually present post-communication `(I,J,R)` key has duplicate MPI contributors, and reports the offending integer key. This check proves uniqueness only for keys present in the post-communication map; it does not assert complete R coverage or pre-communication ownership. The analyzer uses `basis_aux_shrink_out` only to validate per-type shell layout and its declared total; that file does not provide atom-to-type mapping. Full-unshrunk overlap output is not provided because that basis lifecycle is not available safely here.
+- **Default**: False
+
 ### out_pchg
 
 - **Type**: String

docs/advanced/interface/LibRPA.md

@@ -0,0 +1,97 @@
+# LibRPA exports with SOC and magnetic symmetry
+
+The LCAO RI exporter writes the ABACUS states and auxiliary-basis data needed
+by a compatible LibRPA reader. For SOC, the consumer must support spinor
+wavefunctions and the optional `spin_symmetry` block described below.
+Enabling an export does not certify symmetry equivalence of a downstream
+RPA/GW calculation.
+
+## Producer settings
+
+For an otherwise complete, converged LCAO SOC input, the relevant settings are:
+
+```text
+basis_type                  lcao
+nspin                       4
+lspinorb                    1
+symmetry                    1
+rpa                         1
+out_ri_cv                   1
+out_wfc_lcao                1
+out_mat_xc                  1
+out_librpa_reader_version   1
+```
+
+Set the intended magnetic moments in `STRU`. With `symmetry=1`, an all-zero
+initial magnetic configuration stays zero; it is not automatically replaced
+by nonzero starting moments. A nonmagnetic SOC system uses a grey group,
+whereas an ordered magnetic system uses the operations preserving that
+configuration, possibly combined with time reversal. For an on/off comparison,
+use a common converged density and keep the structure, basis, k mesh and
+numerical thresholds fixed.
+
+## Magnetic operation block in `stru_out`
+
+The existing physical lattice/position units and spatial `row` convention
+are retained. After the atom records, symmetry-enabled output contains:
+
+```text
+N row
+<9 integer rotation entries and 3 fractional translations, repeated N times>
+spin_symmetry GREY 2
+<one antiunitary flag per operation>
+```
+
+For a nonmagnetic SOC system, `GREY=1`; the spatial operations are exported
+once, all with flag `0`, and the consumer generates their time-reversed
+partners. For a magnetic SOC system, `GREY=0`; the unitary operations come
+first with flag `0`, followed by the spatial parts of antiunitary operations
+with flag `1`. The total `N` includes both blocks.
+
+The spin-action source `2` tells the compatible consumer to derive spin
+rotations from the spatial axial-vector rotation. Explicit SU(2) matrix
+entries are therefore not written. The trailer is only added for `nspin=4`
+when the symmetry operation table is exported; ordinary scalar exports retain
+their existing layout.
+
+## Optional auxiliary-overlap diagnostic
+
+To inspect the active, shrunk auxiliary basis, also set:
+
+```text
+shrink_abfs_pca_thr          1e-6
+out_librpa_abf_overlap       1
+```
+
+The threshold is an example, not a convergence recommendation. The option
+requires `rpa=1`, reader version `1`, and a nonnegative shrink threshold.
+It writes `v1_abf_overlap_active_iq_<iq>.dat` in the same active auxiliary
+basis as the associated Coulomb matrices. It does not export the full
+unshrunk overlap. Dense matrices can require substantial memory and disk.
+
+From the ABACUS source tree, inspect a completed export with:
+
+```bash
+python3 tools/analyze_librpa_abf_overlap.py /path/to/export
+python3 tools/test_analyze_librpa_abf_overlap.py
+```
+
+The analyzer checks metadata and matrix properties. The writer rejects
+duplicate MPI contributors for present `(I,J,R)` blocks. Neither check
+establishes complete real-space coverage or physical convergence.
+
+## Occupation precision and validation
+
+`band_out` writes the Fermi energy in Hartree and each occupation as the
+native k-weighted occupation multiplied by the number of exported k points.
+The stream precision is initialized before these values, including the first
+occupation. Previously the first occupation and Fermi energy retained the
+stream's default six significant digits; removing an irreducible k weight
+could turn an occupied state into `0.999999` or `1.000002`.
+
+An export check can divide the printed occupation by the exported k-point
+count and compare it with `OUT.<suffix>/eig_occ.txt`. Validate the first row
+as well as later rows. A successful LibRPA reader/EXX run establishes that the
+files can be consumed; full/reduced-grid RPA or GW equivalence requires a
+separate numerical comparison. This producer check does not change LibRPA's
+archived-input regression workflow.

docs/advanced/interface/index.rst

@@ -19,3 +19,4 @@ Interfaces to Other Softwares
    ShengBTE
    candela
    TB2J
+   LibRPA

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
@@ -3146,6 +3146,22 @@ parameters:
     default_value: "False"
     unit: ""
     availability: ""
+  - name: out_librpa_reader_version
+    category: Output information
+    type: Integer
+    description: |
+      Select the ABACUS output format for files consumed by LibRPA. 0 writes the legacy text files, and 1 writes LibRPA reader-v1 binary files directly.
+    default_value: "0"
+    unit: ""
+    availability: Numerical atomic orbital basis with rpa=True.
+  - name: out_librpa_abf_overlap
+    category: Output information
+    type: Boolean
+    description: |
+      Write dense raw active-ABF q-space overlap matrices as v1_abf_overlap_active_iq_<iq>.dat for PSD diagnostics; this can be high-cost in memory and output size. The writer fail-closes if an actually present post-communication (I,J,R) key has duplicate MPI contributors, and reports the offending integer key. This check proves uniqueness only for keys present in the post-communication map; it does not assert complete R coverage or pre-communication ownership. The analyzer uses basis_aux_shrink_out only for per-type shell-layout and declared-total checks because it has no atom-to-type mapping. Requires rpa=True, out_librpa_reader_version=1, and shrink ABFs; full-unshrunk overlap output is not provided.
+    default_value: "False"
+    unit: ""
+    availability: Numerical atomic orbital basis with rpa=True, reader version 1, and shrink ABFs.
   - name: out_pchg
     category: Output information
     type: String

python/pyabacus/src/ModuleDriver/py_driver.cpp

@@ -425,7 +425,7 @@ CalculationResult PyDriver::run(
         );
 
         // Read structure
-        impl_->ucell_->setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running);
+        impl_->ucell_->setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running, 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

@@ -548,7 +548,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/run_symmetry.cpp

@@ -38,7 +38,7 @@ void calculate()
 	output out;
 	ucell.setup_cell(
 	"STRU", 
-	ofs_running);
+	ofs_running, 0);
 	std::cout << "set up cell classic done." << std::endl;
     symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running);
 	ofs_running.close();

source/source_cell/module_symmetry/symm_analysis.cpp

@@ -288,6 +288,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 (PARAM.inp.nspin == 4)
+    {
+        this->analyze_magnetic_group_nspin4(atoms, st, latvec1);
+    }
+
     // Do this here for debug
     if (PARAM.inp.calculation == "relax")
     {