solve conflicts

Author: abacus_fixer

docs/advanced/elec_properties/hs_matrix.md

@@ -65,12 +65,12 @@ The output of $H(R)$ and $S(R)$ matrices is controlled by [out_hsr](../input_fil
 | --- | --- |
 | `0` | Disabled |
 | `1` | Text CSR; an optional second value controls precision, for example `out_hsr 1 12` |
-| `2` | Reserved for future binary output; not implemented |
+| `2` | Native binary CSR using `.dat` files |
 | `3` | NPZ: `hrs1_nao.npz`, `hrs2_nao.npz` when needed, and `sr_nao.npz` |
 
 The legacy keywords `out_mat_hs2 1 [precision]` and `out_hsr_npz 1` remain supported as aliases for text and NPZ output respectively. If `out_hsr` is present together with either legacy keyword, `out_hsr` takes precedence.
 
-For a multi-k calculation, the files contain the individual real-space blocks stored for the Bravais lattice vectors $R$. For a gamma-only calculation, ABACUS stores the real-space contributions in a folded representation. Both text CSR and NPZ output write this internal representation directly: all stored $R$-space contributions are summed into a single block labelled `R = (0, 0, 0)`.
+For a multi-k calculation, the files contain the individual real-space blocks stored for the Bravais lattice vectors $R$. For a gamma-only calculation, ABACUS stores the real-space contributions in a folded representation. Text CSR, native binary, and NPZ output write this internal representation directly: all stored $R$-space contributions are summed into a single block labelled `R = (0, 0, 0)`.
 
 The folded gamma-only output is sufficient to inspect the matrix used by the gamma-only real-space container, but it does not retain the original lattice-vector resolution and cannot be used to interpolate matrices at arbitrary k points. Terms that are added only while constructing $H(k)$, rather than stored in the internal $H(R)$ container, are not guaranteed to be present. Use [out_hsk](../input_files/input-main.md#out_hsk) when the final $H(\Gamma)$ and $S(\Gamma)$ matrices are required.
 
@@ -89,6 +89,22 @@ In gamma-only mode, every generated file reports one Bravais lattice vector and
 # representation: gamma-only folded matrix; stored R-space contributions are summed into R = (0, 0, 0)
 ```
 
+### Native Binary CSR Format
+
+Set `out_hsr 2` to write the same H(R) and S(R) matrix sets with a `.dat` suffix. For example, an `nspin = 2` calculation writes `hrs1_nao.dat`, `hrs2_nao.dat`, and `sr_nao.dat`. When `out_app_flag` is false, the one-based ionic step is included before `_nao`, for example `hrs1g1_nao.dat` and `srg1_nao.dat`.
+
+Each ionic step is a complete record with no padding or self-describing header:
+
+1. Native `int`: zero-based ionic step, matrix dimension, number of R blocks.
+2. For every R block in lexicographic `(Rx, Ry, Rz)` order, four native `int` values: `Rx`, `Ry`, `Rz`, and `nnz`.
+3. `nnz` matrix values in CSR order. Real matrices use one native `double`; complex matrices use consecutive real and imaginary `double` values.
+4. `nnz` native `int` column indices.
+5. `dimension + 1` native `long long` row pointers.
+
+All R blocks stored by the internal HContainer are present, including blocks with zero nonzero values. The sparse threshold is `1e-10`, matching text CSR output. The format uses the host integer representation and byte order and therefore requires a compatible ABI. It does not contain unit-cell, spin, or matrix-label metadata; those are determined by the calculation input and filename.
+
+When `out_app_flag` is true, the first ionic step truncates the shared file and later ionic steps append complete records. Otherwise every ionic step is written to its own file.
+
 ### NPZ Format
 
 Set `out_hsr 3` to write `hrs1_nao.npz`, `hrs2_nao.npz` when a second spin channel is present, and `sr_nao.npz`. Matrix entry names include the atom-pair indices and the three components of $R$. Multi-k calculations retain the stored $R$ blocks, while gamma-only calculations contain only matrix entry names ending in `_0_0_0`.
@@ -125,7 +141,7 @@ The CSR format stores a sparse m × n matrix M in row form using three arrays (v
 
 ### Precision Control
 
-Use `out_hsr 1 12` to output text CSR files with 12-digit precision (default is 8). Precision is ignored for NPZ output.
+Use `out_hsr 1 12` to output text CSR files with 12-digit precision (default is 8). Precision is ignored for native binary and NPZ output.
 
 For calculations involving ionic movements, the output frequency of the matrix is controlled by [out_freq_ion](../input_files/input-main.md#out_freq_ion) and [out_app_flag](../input_files/input-main.md#out_app_flag). 
 

docs/advanced/input_files/input-main.md

@@ -2037,11 +2037,13 @@
 - **Description**: Output Hamiltonian and overlap matrices in real space, indexed by the Bravais lattice vector R, in the directory OUT.${suffix}. The first integer selects the format:
   - 0: disabled;
   - 1: text CSR output; the optional second integer controls precision and defaults to 8;
-  - 2: reserved for binary output, which is not implemented yet;
+  - 2: native binary CSR output using .dat files;
   - 3: NPZ output using hrs1_nao.npz, hrs2_nao.npz when needed, and sr_nao.npz.
 
   For multi-k calculations, the output contains the individual real-space blocks stored for the Bravais lattice vectors R. For gamma-only calculations, the internal real-space contributions are folded into a single R = (0, 0, 0) block. This folded result cannot recover the original R-resolved contributions or interpolate arbitrary k points. Terms added only while constructing H(k) are not guaranteed to be present.
 
+  For binary output, each file uses the same basename as text output with a .dat suffix. Every native record contains the zero-based ionic step, matrix dimension, and number of R blocks as ints. Each R block contains three int coordinates, an int nonzero count, native double values (real/imaginary double pairs for complex matrices), int column indices, and long long row pointers. Native integer representation and byte order are used. When out_app_flag is true, the first ionic step truncates the file and later steps append complete records.
+
   > Note: In the 3.10-LTS version, the file names are data-HR-sparse_SPIN0.csr and data-SR-sparse_SPIN0.csr, etc.
 - **Default**: 0 8
 - **Unit**: Ry

docs/parameters.yaml

@@ -3145,11 +3145,13 @@ parameters:
       Output Hamiltonian and overlap matrices in real space, indexed by the Bravais lattice vector R, in the directory OUT.${suffix}. The first integer selects the format:
       * 0: disabled;
       * 1: text CSR output; the optional second integer controls precision and defaults to 8;
-      * 2: reserved for binary output, which is not implemented yet;
+      * 2: native binary CSR output using .dat files;
       * 3: NPZ output using hrs1_nao.npz, hrs2_nao.npz when needed, and sr_nao.npz.
 
       For multi-k calculations, the output contains the individual real-space blocks stored for the Bravais lattice vectors R. For gamma-only calculations, the internal real-space contributions are folded into a single R = (0, 0, 0) block. This folded result cannot recover the original R-resolved contributions or interpolate arbitrary k points. Terms added only while constructing H(k) are not guaranteed to be present.
 
+      For binary output, each file uses the same basename as text output with a .dat suffix. Every native record contains the zero-based ionic step, matrix dimension, and number of R blocks as ints. Each R block contains three int coordinates, an int nonzero count, native double values (real/imaginary double pairs for complex matrices), int column indices, and long long row pointers. Native integer representation and byte order are used. When out_app_flag is true, the first ionic step truncates the file and later steps append complete records.
+
       [NOTE] In the 3.10-LTS version, the file names are data-HR-sparse_SPIN0.csr and data-SR-sparse_SPIN0.csr, etc.
     default_value: 0 8
     unit: Ry

source/source_cell/module_neighlist/bin_manager.cpp

@@ -5,6 +5,15 @@
 #include <stdexcept>
 #include "bin_manager.h"
 
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+namespace
+{
+constexpr int neighbor_build_openmp_threshold = 256;
+}
+
 // ========== Bin class implementation ==========
 
 const std::vector<ModuleNeighList::LocalAtomIndex>& Bin::get_atom_indices() const {
@@ -176,6 +185,63 @@ int BinManager::bin_index(int ix, int iy, int iz) const {
     return ix * nbiny_ * nbinz_ + iy * nbinz_ + iz;
 }
 
+template <typename Emit>
+void BinManager::visit_neighbors(const NeighborAtom& atom,
+                                 const std::vector<NeighborAtom>& binned_atoms,
+                                 double sradius2,
+                                 const Emit& emit) const
+{
+    const int ix = std::min(
+        std::max(int((atom.position_x - x_min_) / bin_sizex_), 0),
+        nbinx_ - 1
+    );
+
+    const int iy = std::min(
+        std::max(int((atom.position_y - y_min_) / bin_sizey_), 0),
+        nbiny_ - 1
+    );
+
+    const int iz = std::min(
+        std::max(int((atom.position_z - z_min_) / bin_sizez_), 0),
+        nbinz_ - 1
+    );
+
+    for (int dx = -1; dx <= 1; dx++)
+    {
+        for (int dy = -1; dy <= 1; dy++)
+        {
+            for (int dz = -1; dz <= 1; dz++)
+            {
+                const int jx = ix + dx;
+                const int jy = iy + dy;
+                const int jz = iz + dz;
+
+                if (jx < 0 || jx >= nbinx_ ||
+                    jy < 0 || jy >= nbiny_ ||
+                    jz < 0 || jz >= nbinz_)
+                {
+                    continue;
+                }
+
+                const int nidx = bin_index(jx, jy, jz);
+                for (const ModuleNeighList::LocalAtomIndex binned_atom_index : bins_[nidx].get_atom_indices())
+                {
+                    const NeighborAtom& natom = binned_atoms[static_cast<std::size_t>(binned_atom_index)];
+                    const double delta_x = atom.position_x - natom.position_x;
+                    const double delta_y = atom.position_y - natom.position_y;
+                    const double delta_z = atom.position_z - natom.position_z;
+                    const double dist2 = delta_x * delta_x + delta_y * delta_y + delta_z * delta_z;
+
+                    if (natom.atom_id != atom.atom_id && dist2 <= sradius2)
+                    {
+                        emit(natom.atom_id);
+                    }
+                }
+            }
+        }
+    }
+}
+
 void BinManager::build_atom_neighbors(
     NeighborList& neighbor_list,
     const std::vector<NeighborAtom>& atoms,
@@ -184,71 +250,63 @@ void BinManager::build_atom_neighbors(
 {
     assert(atoms.size() == static_cast<size_t>(neighbor_list.get_nlocal()));
 
-    double sradius2 = sradius_ * sradius_;
+    const double sradius2 = sradius_ * sradius_;
 
     neighbor_list.reset();
 
-    std::vector<int> neigh_tmp;
-
     const int nlocal = neighbor_list.get_nlocal();
-    for (int i = 0; i < nlocal; i++)
-    {
-        neigh_tmp.clear();
-        const NeighborAtom& atom = atoms[i];
-
-        int ix = std::min(
-            std::max(int((atom.position_x - x_min_) / bin_sizex_), 0),
-            nbinx_ - 1
-        );
 
-        int iy = std::min(
-            std::max(int((atom.position_y - y_min_) / bin_sizey_), 0),
-            nbiny_ - 1
-        );
-
-        int iz = std::min(
-            std::max(int((atom.position_z - z_min_) / bin_sizez_), 0),
-            nbinz_ - 1
-        );
+#ifdef _OPENMP
+    const bool use_parallel = nlocal >= neighbor_build_openmp_threshold && omp_get_max_threads() > 1;
+    if (use_parallel)
+    {
+        std::vector<std::size_t> neighbor_counts(static_cast<std::size_t>(nlocal), 0);
 
-        for (int dx = -1; dx <= 1; dx++)
+#pragma omp parallel for schedule(static)
+        for (int i = 0; i < nlocal; i++)
         {
-            for (int dy = -1; dy <= 1; dy++)
-            {
-                for (int dz = -1; dz <= 1; dz++)
-                {
-                    int jx = ix + dx;
-                    int jy = iy + dy;
-                    int jz = iz + dz;
-
-                    if (jx < 0 || jx >= nbinx_ ||
-                        jy < 0 || jy >= nbiny_ ||
-                        jz < 0 || jz >= nbinz_)
-                        continue;
-
-                    int nidx = bin_index(jx, jy, jz);
+            std::size_t count = 0;
+            visit_neighbors(atoms[i], binned_atoms, sradius2,
+                            [&count](ModuleNeighList::LocalAtomIndex) { ++count; });
+            neighbor_counts[static_cast<std::size_t>(i)] = count;
+        }
 
-                    for (const ModuleNeighList::LocalAtomIndex binned_atom_index : bins_[nidx].get_atom_indices())
-                    {
-                        const NeighborAtom& natom = binned_atoms[static_cast<std::size_t>(binned_atom_index)];
-                        double dx = atom.position_x - natom.position_x;
-                        double dy = atom.position_y - natom.position_y;
-                        double dz = atom.position_z - natom.position_z;
+        for (int i = 0; i < nlocal; i++)
+        {
+            const int n = ModuleNeighList::checked_int_size(
+                neighbor_counts[static_cast<std::size_t>(i)],
+                "BinManager neighbor count"
+            );
+            neighbor_list.firstneigh_[i] = neighbor_list.allocator_.allocate(n);
+            neighbor_list.numneigh_[i] = n;
+        }
 
-                        double dist2 = dx * dx + dy * dy + dz * dz;
+#pragma omp parallel for schedule(static)
+        for (int i = 0; i < nlocal; i++)
+        {
+            int* ptr = neighbor_list.firstneigh_[i];
+            int k = 0;
+            visit_neighbors(atoms[i], binned_atoms, sradius2,
+                            [&](ModuleNeighList::LocalAtomIndex atom_id)
+                            {
+                                assert(ptr != nullptr);
+                                ptr[k++] = atom_id;
+                            });
+            assert(k == neighbor_list.numneigh_[i]);
+        }
+        return;
+    }
+#endif
 
-                        if (natom.atom_id == atom.atom_id)
-                        {
-                            continue;
-                        }
-                        if (dist2 <= sradius2)
+    std::vector<int> neigh_tmp;
+    for (int i = 0; i < nlocal; i++)
+    {
+        neigh_tmp.clear();
+        visit_neighbors(atoms[i], binned_atoms, sradius2,
+                        [&neigh_tmp](ModuleNeighList::LocalAtomIndex atom_id)
                         {
-                            neigh_tmp.push_back(natom.atom_id);
-                        }
-                    }
-                }
-            }
-        }
+                            neigh_tmp.push_back(atom_id);
+                        });
 
         const int n = ModuleNeighList::checked_int_size(neigh_tmp.size(), "BinManager neighbor count");
 

source/source_cell/module_neighlist/bin_manager.h

@@ -198,6 +198,21 @@ class BinManager
      * @return Flat index in the bins_ array.
      */
     int bin_index(int ix, int iy, int iz) const;
+
+    /**
+     * @brief Visit neighbors of one atom in the existing deterministic bin order.
+     *
+     * @tparam Emit Callable accepting a rank-local neighbor atom ID.
+     * @param atom Atom used as the neighbor-list center.
+     * @param binned_atoms All atoms assigned to bins by do_binning().
+     * @param sradius2 Squared search radius.
+     * @param emit Callback invoked once for every accepted neighbor.
+     */
+    template <typename Emit>
+    void visit_neighbors(const NeighborAtom& atom,
+                         const std::vector<NeighborAtom>& binned_atoms,
+                         double sradius2,
+                         const Emit& emit) const;
 };
 
 #endif // BIN_MANAGER_H