Merge branch 'develop' into 2026-07-30-c

Author: mohanchen

Dockerfile.cuda

@@ -24,6 +24,13 @@ RUN cd /tmp && \
     ln -s /usr/local/include/elpa_openmp-$ELPA_VER/elpa /usr/local/include/ && \
     cd /tmp && rm -rf elpa-$ELPA_VER
 
+# RapidJSON
+RUN cd /tmp && wget --quiet https://codeload.github.com/Tencent/rapidjson/tar.gz/24b5e7a -O rapidjson-24b5e7a.tar.gz
+RUN tar -xzf rapidjson-24b5e7a.tar.gz && cd rapidjson-24b5e7a
+RUN cmake -B build -DRAPIDJSON_BUILD_DOC=OFF -DRAPIDJSON_BUILD_EXAMPLES=OFF -DRAPIDJSON_BUILD_TESTS=OFF
+RUN cmake --build build --target install
+RUN cd /tmp && rm -r rapidjson-24b5e7a
+
 ADD https://api.github.com/repos/deepmodeling/abacus-develop/git/refs/heads/develop /dev/null
 
 RUN git clone https://github.com/deepmodeling/abacus-develop.git --depth 1 && \

Dockerfile.gnu

@@ -26,6 +26,13 @@ RUN wget https://download.pytorch.org/libtorch/cpu/libtorch-cxx11-abi-shared-wit
         --no-check-certificate --quiet -O libtorch.zip && \
     unzip -q libtorch.zip -d /opt && rm libtorch.zip
 
+# RapidJSON
+RUN cd /tmp && wget --quiet https://codeload.github.com/Tencent/rapidjson/tar.gz/24b5e7a -O rapidjson-24b5e7a.tar.gz
+RUN tar -xzf rapidjson-24b5e7a.tar.gz && cd rapidjson-24b5e7a
+RUN cmake -B build -DRAPIDJSON_BUILD_DOC=OFF -DRAPIDJSON_BUILD_EXAMPLES=OFF -DRAPIDJSON_BUILD_TESTS=OFF
+RUN cmake --build build --target install
+RUN cd /tmp && rm -r rapidjson-24b5e7a
+
 ENV CMAKE_PREFIX_PATH=/opt/libtorch/share/cmake
 
 ADD https://api.github.com/repos/deepmodeling/abacus-develop/git/refs/heads/develop /dev/null

Dockerfile.intel

@@ -1,63 +1,12 @@
-FROM ubuntu:22.04
+FROM intel/oneapi-hpckit:2025.2.2-0-devel-ubuntu22.04
 
 RUN apt-get update && apt-get install -y \
     bc cmake git gnupg gcc g++ python3-numpy sudo wget vim unzip \
     libcereal-dev libxc-dev libgtest-dev libgmock-dev libbenchmark-dev \
     pkg-config build-essential autoconf automake libtool
 
-# Following steps by https://software.intel.com/content/www/us/en/develop/documentation/installation-guide-for-intel-oneapi-toolkits-linux/top/installation/install-using-package-managers/apt.html .
-RUN wget -O- https://apt.repos.intel.com/intel-gpg-keys/GPG-PUB-KEY-INTEL-SW-PRODUCTS.PUB \
-    | gpg --dearmor | sudo tee /usr/share/keyrings/oneapi-archive-keyring.gpg > /dev/null && \
-    echo "deb [signed-by=/usr/share/keyrings/oneapi-archive-keyring.gpg] https://apt.repos.intel.com/oneapi all main" \
-    | sudo tee /etc/apt/sources.list.d/oneAPI.list
-
-# To save disk space, only install the essential components, but not the whole toolkit.
-RUN apt-get update && \
-    apt-get install -y \
-        intel-oneapi-compiler-dpcpp-cpp \
-        intel-oneapi-compiler-fortran \
-        intel-oneapi-mkl-devel \
-        intel-oneapi-mkl-sycl-devel \
-        intel-oneapi-mkl-sycl-distributed-dft-devel \
-        intel-oneapi-mpi-devel \
-        intel-oneapi-vtune && \
-    apt-get clean && \
-    rm -rf /var/lib/apt/lists/*
-
-RUN ls -la /opt/intel/oneapi/mkl/latest/lib/intel64/*sycl*dft*
-RUN ls -la /opt/intel/oneapi/mkl/latest/lib/intel64/*dis*
-
-# Set oneAPI environment variables
-ENV ONEAPI_ROOT=/opt/intel/oneapi
-ENV I_MPI_ROOT=${ONEAPI_ROOT}/mpi/latest
-ENV MKLROOT=${ONEAPI_ROOT}/mkl/latest
-ENV CMPLR_ROOT=${ONEAPI_ROOT}/compiler/latest
-
-# Set library paths and include paths
-ENV LIBRARY_PATH=${ONEAPI_ROOT}/tbb/latest/lib/intel64/gcc4.8:${ONEAPI_ROOT}/mpi/latest/lib:${MKLROOT}/lib/intel64:${ONEAPI_ROOT}/compiler/latest/lib/:${LIBRARY_PATH}
-ENV LD_LIBRARY_PATH=${ONEAPI_ROOT}/tbb/latest/lib/intel64/gcc4.8:${ONEAPI_ROOT}/mpi/latest/lib:${MKLROOT}/lib/intel64:${ONEAPI_ROOT}/compiler/latest/lib/:${LD_LIBRARY_PATH}
-ENV PATH=${ONEAPI_ROOT}/vtune/latest/bin64:${ONEAPI_ROOT}/mpi/latest/bin:${MKLROOT}/bin/intel64:${ONEAPI_ROOT}/compiler/latest/bin:${PATH}
-ENV CPATH=${MKLROOT}/include:${ONEAPI_ROOT}/mpi/latest/include:${CPATH}
-ENV PKG_CONFIG_PATH=${MKLROOT}/lib/pkgconfig:${ONEAPI_ROOT}/mpi/latest/lib/pkgconfig:${PKG_CONFIG_PATH}
-
-# Set CMAKE related paths
-ENV CMAKE_PREFIX_PATH=${ONEAPI_ROOT}/tbb/latest:${MKLROOT}/lib/cmake:${ONEAPI_ROOT}/dpl/latest/lib/cmake:${ONEAPI_ROOT}/dnnl/latest/lib/cmake:${ONEAPI_ROOT}/dal/latest:${ONEAPI_ROOT}/compiler/latest:${CMAKE_PREFIX_PATH}
-
-SHELL ["/bin/bash", "-c"]
-ENV CC=mpiicx CXX=mpiicpx FC=mpiifx
-
-# Verify oneAPI installation
-RUN source ${ONEAPI_ROOT}/setvars.sh && \
-    echo "=== Verify compiler ===" && \
-    which mpiicx && mpiicx --version && \
-    echo "=== Verify MKL ===" && \
-    ls ${MKLROOT}/lib/intel64/ && \
-    echo "=== Verify MPI ===" && \
-    which mpirun && mpirun --version
-
 # https://elpa.mpcdf.mpg.de/software/tarball-archive/ELPA_TARBALL_ARCHIVE.html
-RUN source /opt/intel/oneapi/setvars.sh && \
-    cd /tmp && \
+RUN cd /tmp && \
     ELPA_VER=2022.11.001 && \
     wget -q https://elpa.mpcdf.mpg.de/software/tarball-archive/Releases/$ELPA_VER/elpa-$ELPA_VER.tar.gz && \
     tar xzf elpa-$ELPA_VER.tar.gz  && rm elpa-$ELPA_VER.tar.gz && \
@@ -68,8 +17,14 @@ RUN source /opt/intel/oneapi/setvars.sh && \
     ln -s /usr/local/include/elpa_openmp-$ELPA_VER/elpa /usr/local/include/ && \
     cd /tmp && rm -rf elpa-$ELPA_VER
 
-# rapidjson and libtorch
-RUN cd /tmp && git clone --depth 1 https://github.com/Tencent/rapidjson.git && cp -r rapidjson/include/rapidjson /usr/include/ && rm -rf rapidjson
+# RapidJSON
+RUN cd /tmp && wget --quiet https://codeload.github.com/Tencent/rapidjson/tar.gz/24b5e7a -O rapidjson-24b5e7a.tar.gz
+RUN tar -xzf rapidjson-24b5e7a.tar.gz && cd rapidjson-24b5e7a
+RUN cmake -B build -DRAPIDJSON_BUILD_DOC=OFF -DRAPIDJSON_BUILD_EXAMPLES=OFF -DRAPIDJSON_BUILD_TESTS=OFF
+RUN cmake --build build --target install
+RUN cd /tmp && rm -r rapidjson-24b5e7a
+
+# LibTorch (Note: Using pre-built Torch library with MKL might cause issues)
 RUN wget -q https://download.pytorch.org/libtorch/cpu/libtorch-cxx11-abi-shared-with-deps-2.0.0%2Bcpu.zip -O /tmp/libtorch.zip && \
     unzip -q /tmp/libtorch.zip -d /opt && rm -f /tmp/libtorch.zip
 ENV CMAKE_PREFIX_PATH=/opt/libtorch/share/cmake:${CMAKE_PREFIX_PATH}

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;
+}
+