feat: support native binary H(R)/S(R) output (#7806)

Author: goodchong

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_io/module_ctrl/ctrl_scf_lcao.cpp

@@ -267,13 +267,13 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell,
     //------------------------------------------------------------------
     //! 7a) Output H(R) and S(R) matrices in CSR format
     //------------------------------------------------------------------
-    if (inp.out_hsr[0] == 1)
+    if (inp.out_hsr[0] == 1 || inp.out_hsr[0] == 2)
     {
         const int precision = inp.out_hsr[1];
         std::vector<hamilt::HContainer<TR>*> hr_vec = p_hamilt->getHR_vector();
         const hamilt::HContainer<TR>* sr = p_hamilt->getSR();
 
-        ModuleIO::write_hsr(hr_vec, sr, &ucell, precision, pv,
+        ModuleIO::write_hsr(hr_vec, sr, &ucell, inp.out_hsr[0], precision, pv,
                             out_app_flag, gamma_only, ucell.get_iat2iwt(), ucell.nat, istep);
     }
 

source/source_io/module_hs/write_HS_R.cpp

@@ -1,5 +1,6 @@
 #include "write_HS_R.h"
 
+#include "source_base/module_out/sparse_matrix.h"
 #include "source_base/timer.h"
 #include "source_base/tool_quit.h"
 #include "source_io/module_parameter/parameter.h"
@@ -9,6 +10,13 @@
 #include "source_lcao/spar_st.h"
 #include "write_HS_sparse.h"
 
+#include <algorithm>
+#include <complex>
+#include <fstream>
+#include <limits>
+#include <tuple>
+#include <vector>
+
 // if 'binary=true', output binary file.
 // The 'sparse_thr' is the accuracy of the sparse matrix.
 // If the absolute value of the matrix element is less than or equal to the
@@ -229,23 +237,39 @@ std::string ModuleIO::hsr_gen_fname(const std::string& prefix,
                                      const bool append,
                                      const int istep)
 {
+    return hsr_gen_fname(prefix, ispin, append, istep, 1);
+}
+
+std::string ModuleIO::hsr_gen_fname(const std::string& prefix,
+                                     const int ispin,
+                                     const bool append,
+                                     const int istep,
+                                     const int out_type)
+{
+    const std::string extension = out_type == 2 ? ".dat" : ".csr";
     if (!append && istep >= 0)
     {
-        return prefix + std::to_string(ispin + 1) + "g" + std::to_string(istep + 1) + "_nao.csr";
+        return prefix + std::to_string(ispin + 1) + "g" + std::to_string(istep + 1) + "_nao" + extension;
     }
     else
     {
-        return prefix + std::to_string(ispin + 1) + "_nao.csr";
+        return prefix + std::to_string(ispin + 1) + "_nao" + extension;
     }
 }
 
 std::string ModuleIO::sr_gen_fname(const bool append, const int istep)
 {
+    return sr_gen_fname(append, istep, 1);
+}
+
+std::string ModuleIO::sr_gen_fname(const bool append, const int istep, const int out_type)
+{
+    const std::string extension = out_type == 2 ? ".dat" : ".csr";
     if (!append && istep >= 0)
     {
-        return "srg" + std::to_string(istep + 1) + "_nao.csr";
+        return "srg" + std::to_string(istep + 1) + "_nao" + extension;
     }
-    return "sr_nao.csr";
+    return "sr_nao" + extension;
 }
 
 std::string ModuleIO::dhr_gen_fname(const std::string& prefix,
@@ -309,10 +333,159 @@ void ModuleIO::write_hcontainer_csr(const std::string& fname,
     ofs.close();
 }
 
+namespace
+{
+template <typename T>
+void write_native_value(std::ofstream& ofs, const T& value)
+{
+    ofs.write(reinterpret_cast<const char*>(&value), sizeof(T));
+}
+
+void write_native_value(std::ofstream& ofs, const std::complex<double>& value)
+{
+    const double real = value.real();
+    const double imag = value.imag();
+    write_native_value(ofs, real);
+    write_native_value(ofs, imag);
+}
+
+template <typename TR>
+ModuleIO::SparseMatrix<TR> make_sparse_R_block(hamilt::HContainer<TR>* mat_serial,
+                                               const int rx,
+                                               const int ry,
+                                               const int rz,
+                                               const double sparse_threshold)
+{
+    const int nbasis = mat_serial->get_nbasis();
+    ModuleIO::SparseMatrix<TR> sparse_matrix(nbasis, nbasis);
+    sparse_matrix.setSparseThreshold(sparse_threshold);
+    mat_serial->fix_R(rx, ry, rz);
+
+    for (int iap = 0; iap < mat_serial->size_atom_pairs(); ++iap)
+    {
+        const auto atom_pair = mat_serial->get_atom_pair(iap);
+        const int r_index = atom_pair.find_R(rx, ry, rz);
+        if (r_index < 0)
+        {
+            continue;
+        }
+        const auto matrix_info = atom_pair.get_matrix_values(r_index);
+        const int* index = std::get<0>(matrix_info).data();
+        TR* data = std::get<1>(matrix_info);
+        for (int irow = index[0]; irow < index[0] + index[1]; ++irow)
+        {
+            for (int icol = index[2]; icol < index[2] + index[3]; ++icol)
+            {
+                sparse_matrix.insert(irow, icol, *data);
+                ++data;
+            }
+        }
+    }
+
+    mat_serial->unfix_R();
+    return sparse_matrix;
+}
+} // namespace
+
+template <typename TR>
+void ModuleIO::write_hcontainer_csr_binary(const std::string& fname,
+                                            hamilt::HContainer<TR>* mat_serial,
+                                            const int istep,
+                                            const bool append)
+{
+    std::ios_base::openmode mode = std::ios::out | std::ios::binary;
+    if (append && istep > 0)
+    {
+        mode |= std::ios::app;
+    }
+    std::ofstream ofs(fname.c_str(), mode);
+    if (!ofs.is_open())
+    {
+        ModuleBase::WARNING_QUIT("ModuleIO::write_hcontainer_csr_binary",
+                                 "Cannot open HContainer binary CSR file: " + fname);
+    }
+
+    const size_t nR_size = mat_serial->size_R_loop();
+    if (nR_size > static_cast<size_t>(std::numeric_limits<int>::max()))
+    {
+        ModuleBase::WARNING_QUIT("ModuleIO::write_hcontainer_csr_binary",
+                                 "Too many R blocks for native binary CSR output: " + fname);
+    }
+
+    std::vector<std::tuple<int, int, int>> R_coordinates;
+    R_coordinates.reserve(nR_size);
+    for (size_t iR = 0; iR < nR_size; ++iR)
+    {
+        int rx = 0;
+        int ry = 0;
+        int rz = 0;
+        mat_serial->loop_R(iR, rx, ry, rz);
+        R_coordinates.push_back(std::make_tuple(rx, ry, rz));
+    }
+    std::sort(R_coordinates.begin(), R_coordinates.end());
+
+    const int step = std::max(istep, 0);
+    const int nbasis = mat_serial->get_nbasis();
+    const int nR = static_cast<int>(R_coordinates.size());
+    write_native_value(ofs, step);
+    write_native_value(ofs, nbasis);
+    write_native_value(ofs, nR);
+
+    const double sparse_threshold = 1e-10;
+    for (const auto& R_coordinate: R_coordinates)
+    {
+        const int rx = std::get<0>(R_coordinate);
+        const int ry = std::get<1>(R_coordinate);
+        const int rz = std::get<2>(R_coordinate);
+        const SparseMatrix<TR> sparse_matrix
+            = make_sparse_R_block(mat_serial, rx, ry, rz, sparse_threshold);
+        const auto& elements = sparse_matrix.getElements();
+        if (elements.size() > static_cast<size_t>(std::numeric_limits<int>::max()))
+        {
+            ModuleBase::WARNING_QUIT("ModuleIO::write_hcontainer_csr_binary",
+                                     "Too many nonzero values for native binary CSR output: " + fname);
+        }
+        const int nnz = static_cast<int>(elements.size());
+
+        write_native_value(ofs, rx);
+        write_native_value(ofs, ry);
+        write_native_value(ofs, rz);
+        write_native_value(ofs, nnz);
+
+        for (const auto& element: elements)
+        {
+            write_native_value(ofs, element.second);
+        }
+
+        std::vector<long long> row_ptr(nbasis + 1, 0);
+        for (const auto& element: elements)
+        {
+            write_native_value(ofs, element.first.second);
+            ++row_ptr[element.first.first + 1];
+        }
+        for (int irow = 1; irow <= nbasis; ++irow)
+        {
+            row_ptr[irow] += row_ptr[irow - 1];
+        }
+        for (const long long pointer: row_ptr)
+        {
+            write_native_value(ofs, pointer);
+        }
+    }
+
+    ofs.close();
+    if (!ofs)
+    {
+        ModuleBase::WARNING_QUIT("ModuleIO::write_hcontainer_csr_binary",
+                                 "Failed to write HContainer binary CSR file: " + fname);
+    }
+}
+
 template <typename TR>
 void ModuleIO::write_hsr(const std::vector<hamilt::HContainer<TR>*>& hr_vec,
                           const hamilt::HContainer<TR>* sr,
                           const UnitCell* ucell,
+                          const int out_type,
                           const int precision,
                           const Parallel_2D& paraV,
                           const bool append,
@@ -321,6 +494,10 @@ void ModuleIO::write_hsr(const std::vector<hamilt::HContainer<TR>*>& hr_vec,
                           const int nat,
                           const int istep)
 {
+    if (out_type != 1 && out_type != 2)
+    {
+        ModuleBase::WARNING_QUIT("ModuleIO::write_hsr", "out_type must be 1 or 2");
+    }
     const int nspin = hr_vec.size();
     assert(nspin > 0);
     const std::string representation_note
@@ -346,9 +523,16 @@ void ModuleIO::write_hsr(const std::vector<hamilt::HContainer<TR>*>& hr_vec,
         if (GlobalV::MY_RANK == 0)
         {
             std::string fname = PARAM.globalv.global_out_dir
-                                + hsr_gen_fname("hrs", ispin, append, istep);
-            write_hcontainer_csr(
-                fname, ucell, precision, &hr_serial, istep, ispin, nspin, "H", representation_note);
+                                + hsr_gen_fname("hrs", ispin, append, istep, out_type);
+            if (out_type == 2)
+            {
+                write_hcontainer_csr_binary(fname, &hr_serial, istep, append);
+            }
+            else
+            {
+                write_hcontainer_csr(
+                    fname, ucell, precision, &hr_serial, istep, ispin, nspin, "H", representation_note);
+            }
         }
     }
 
@@ -369,9 +553,16 @@ void ModuleIO::write_hsr(const std::vector<hamilt::HContainer<TR>*>& hr_vec,
         if (GlobalV::MY_RANK == 0)
         {
             std::string fname = PARAM.globalv.global_out_dir
-                                + sr_gen_fname(append, istep);
-            write_hcontainer_csr(
-                fname, ucell, precision, &sr_serial, istep, 0, 1, "S", representation_note);
+                                + sr_gen_fname(append, istep, out_type);
+            if (out_type == 2)
+            {
+                write_hcontainer_csr_binary(fname, &sr_serial, istep, append);
+            }
+            else
+            {
+                write_hcontainer_csr(
+                    fname, ucell, precision, &sr_serial, istep, 0, 1, "S", representation_note);
+            }
         }
     }
 }
@@ -384,15 +575,20 @@ template void ModuleIO::write_hcontainer_csr<std::complex<double>>(
     const std::string&, const UnitCell*, const int,
     hamilt::HContainer<std::complex<double>>*, const int, const int, const int, const std::string&, const std::string&);
 
+template void ModuleIO::write_hcontainer_csr_binary<double>(
+    const std::string&, hamilt::HContainer<double>*, const int, const bool);
+template void ModuleIO::write_hcontainer_csr_binary<std::complex<double>>(
+    const std::string&, hamilt::HContainer<std::complex<double>>*, const int, const bool);
+
 template void ModuleIO::write_hsr<double>(
     const std::vector<hamilt::HContainer<double>*>&,
     const hamilt::HContainer<double>*,
-    const UnitCell*, const int, const Parallel_2D&,
+    const UnitCell*, const int, const int, const Parallel_2D&,
     const bool, const bool, const int*, const int, const int);
 template void ModuleIO::write_hsr<std::complex<double>>(
     const std::vector<hamilt::HContainer<std::complex<double>>*>&,
     const hamilt::HContainer<std::complex<double>>*,
-    const UnitCell*, const int, const Parallel_2D&,
+    const UnitCell*, const int, const int, const Parallel_2D&,
     const bool, const bool, const int*, const int, const int);