Feature: use new format of CSR file for out_mat_hs2

Author: dyzheng

docs/advanced/elec_properties/hs_matrix.md

@@ -26,24 +26,47 @@ The rest of the file contains the upper triangular part of the specified matrice
 
 The output of $H(R)$ and $S(R)$ matrices is controlled by the keyword [out_mat_hs2](../input_files/input-main.md#out_mat_hs2). This functionality is not available for gamma_only calculations. To generate such matrices for gamma only calculations, users should turn off [gamma_only](../input_files/input-main.md#gamma_only), and explicitly specify that gamma point is the only k point in the KPT file.
 
-For single-point SCF calculations, if nspin = 1 or nspin = 4, two files `hrs1_nao.csr` and `sr_nao.csr` are generated, which contain the Hamiltonian matrix $H(R)$ and overlap matrix $S(R)$ respectively. For nspin = 2, three files `hrs1_nao.csr` and `hrs2_nao.csr` and `sr_nao.csr` are created, where the first two files correspodn to $H(R)$ for spin up and spin down, respectively.
+### Output Format
 
-Each file or each section of the appended file starts with three lines, the first gives the current ion/md step, the second gives the dimension of the matrix, and the last indicates how many different `R` are in the file.
+The H(R) and S(R) matrices are output in standard Compressed Sparse Row (CSR) format, matching the format used by `out_dmr`.
 
-The rest of the files are arranged in blocks. Each block starts with a line giving the lattice vector `R` and the number of nonzero matrix elements, such as:
+For single-point SCF calculations:
+- **nspin = 1 or nspin = 4**: Two files `hrs1_nao.csr` and `srs1_nao.csr` are generated, containing the Hamiltonian matrix $H(R)$ and overlap matrix $S(R)$ respectively.
+- **nspin = 2**: Three files `hrs1_nao.csr`, `hrs2_nao.csr`, and `srs1_nao.csr` are created, where the first two files correspond to $H(R)$ for spin up and spin down, respectively.
 
+### File Structure
+
+Each file starts with a header:
 ```
--3 1 1 1020
+ --- Ionic Step 1 ---
+ # print H matrix in real space H(R)
+ 1 # number of spin directions
+ 1 # spin index
+ 100 # number of localized basis
+ 50 # number of Bravais lattice vector R
+
+[UnitCell information]
+
+#----------------------------------------------------------------------#
+#                               CSR Format                             #
+...
+ 0 0 0 5
+ # CSR values
+ 1.234e-01 2.345e-02 ...
+ # CSR column indices
+ 0 5 10 ...
+ # CSR row pointers
+ 0 3 7 ...
 ```
 
-which means there are 1020 nonzero elements in the (-3,1,1) cell.
+The CSR format stores a sparse m × n matrix M in row form using three arrays (values, column indices, row pointers). According to Wikipedia:
 
-If there is no nonzero matrix element, then the next block starts immediately on the next line. Otherwise, there will be 3 extra lines in the block, which gives the matrix in CSR format. According to Wikipedia:
+- The arrays **values** and **column indices** are of length NNZ (number of nonzero entries), and contain the non-zero values and the column indices of those values respectively.
+- The array **row pointers** is of length m + 1 and encodes the index where each row starts. The last element is NNZ.
 
-The CSR format stores a sparse m × n matrix M in row form using three (one-dimensional) arrays (V, COL_INDEX, ROW_INDEX). Let NNZ denote the number of nonzero entries in M. (Note that zero-based indices shall be used here.)
+### Precision Control
 
-  - The arrays V and COL_INDEX are of length NNZ, and contain the non-zero values and the column indices of those values respectively.
-  - The array ROW_INDEX is of length m + 1 and encodes the index in V and COL_INDEX where the given row starts. This is equivalent to ROW_INDEX[j] encoding the total number of nonzeros above row j. The last element is NNZ , i.e., the fictitious index in V immediately after the last valid index NNZ - 1.
+Use `out_mat_hs2 1 12` to output with 12-digit precision (default is 8).
 
 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

@@ -1916,13 +1916,15 @@
 
 ### out_mat_hs2
 
-- **Type**: Boolean
+- **Type**: Boolean \[Integer\](optional)
 - **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)*
-- **Description**: Whether to print files containing the Hamiltonian matrix and overlap matrix into files in the directory OUT.${suffix}. For more information, please refer to hs_matrix.md.
+- **Description**: Output Hamiltonian H(R) and overlap S(R) matrices in CSR format. Optional second parameter specifies output precision (default 8). For more information, please refer to hs_matrix.md.
 
   > Note: In the 3.10-LTS version, the file names are data-HR-sparse_SPIN0.csr and data-SR-sparse_SPIN0.csr, etc.
-- **Default**: False
+- **Default**: False [8]
 - **Unit**: Ry
+- **Output files**: `hrs1_nao.csr`, `hrs2_nao.csr` (if nspin=2), `srs1_nao.csr`
+- **Format**: Standard CSR format with UnitCell header, compatible with `out_dmr` format
 
 ### out_mat_tk
 

docs/advanced/interface/TB2J.md

@@ -122,15 +122,17 @@ After the key parameter `out_mat_hs2` is turned on, the Hamiltonian matrix $H(R)
 suffix                         Fe
 ```
 
-specifies the suffix of the output, in this calculation, we set the path to the directory of the DFT calculation, which is the current directory (".") and the suffix to Fe. 
+specifies the suffix of the output, in this calculation, we set the path to the directory of the DFT calculation, which is the current directory (".") and the suffix to Fe.
+
+> **Note (ABACUS v3.9.0.25+):** Starting from ABACUS v3.9.0.25, the output format has changed to standard CSR format with filenames `hrs1_nao.csr`, `hrs2_nao.csr` (for nspin=2), and `srs1_nao.csr`. The parameter `out_mat_hs2` now supports optional precision control: `out_mat_hs2 1 8` (default 8 digits). TB2J v0.9.0+ is required to read the new format. For older TB2J versions, please use ABACUS v3.8.x or earlier. 
 
 #### 2. Perform TB2J calculation:
 
 ```bash
 abacus2J.py --path . --suffix Fe --elements Fe  --kmesh 7 7 7
 ```
 
-This first read the atomic structures from th `STRU` file,  then read the Hamiltonian and the overlap matrices stored in the files named starting from `data-HR-*` and `data-SR-*` files.  It also read the fermi energy from the `OUT.Fe/running_scf.log` file.  
+This first reads the atomic structures from the `STRU` file, then reads the Hamiltonian and overlap matrices. For ABACUS v3.9.0.25+, the matrices are stored in `hrs1_nao.csr`, `hrs2_nao.csr` (nspin=2), and `srs1_nao.csr` files. For older versions, they are in `data-HR-*` and `data-SR-*` files. It also reads the fermi energy from the `OUT.Fe/running_scf.log` file.  
 
 With the command above, we can calculate the $J$ with a $7 \times 7 \times 7$ k-point grid. This allows for the calculation of exchange between spin pairs between $7 \times 7 \times 7$  supercell.  Note: the kmesh is not dense enough for a practical calculation. For a very dense k-mesh, the `--rcut` option can be used to set the maximum distance of the magnetic interactions and thus reduce the computation cost. But be sure that the cutoff is not too small. 
 

docs/advanced/interface/deeph.md

@@ -16,7 +16,13 @@ The first stage is during the data preparation phase, where we need to run a ser
 out_mat_hs2 1
 ```
 
-Files named data-HR-sparse_SPIN`${x}`.csr and data-SR-sparse_SPIN`${x}`.csr will be generated, which contain the Hamiltonian and overlap matrices respectively in csr format. `${x}` takes value of 0 or 1, based on the spin component. More details on this keyword can be found in the [list of input keywords](../input_files/input-main.md#out_mat_hs2).
+**For ABACUS v3.9.0.25+:** Files named `hrs1_nao.csr`, `hrs2_nao.csr` (for nspin=2), and `srs1_nao.csr` will be generated in `OUT.${suffix}/` directory, containing the Hamiltonian and overlap matrices in standard CSR format. You can optionally specify precision: `out_mat_hs2 1 8` (default 8 digits).
+
+**For ABACUS v3.8.x and earlier:** Files named `data-HR-sparse_SPIN${x}.csr` and `data-SR-sparse_SPIN${x}.csr` will be generated, where `${x}` takes value of 0 or 1 based on the spin component.
+
+> **Note:** DeepH v1.0.0+ is required to read the new CSR format from ABACUS v3.9.0.25+. For older DeepH versions, please use ABACUS v3.8.x or earlier.
+
+More details on this keyword can be found in the [list of input keywords](../input_files/input-main.md#out_mat_hs2).
 
 The second stage is during the inference phase. After DeepH training completes, we can apply the model to predict the Hamiltonian on other systems. For that purpose, we also need the overlap matrices from the new systems, but no SCF calculation is required.