deepmodeling · pincher-chen · Jun 16, 2026
diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
@@ -177,6 +177,9 @@
     - [out\_mat\_l](#out_mat_l)
     - [out\_xc\_r](#out_xc_r)
     - [out\_eband\_terms](#out_eband_terms)
+    - [out\_hr\_npz](#out_hr_npz)
+    - [out\_hsr\_npz](#out_hsr_npz)
+    - [out\_dm\_npz](#out_dm_npz)
     - [out\_mul](#out_mul)
     - [out\_app\_flag](#out_app_flag)
     - [out\_ndigits](#out_ndigits)
@@ -2049,6 +2052,28 @@
 - **Description**: Whether to print the band energy terms separately in the file OUT.{term}_out.dat. The terms include the kinetic, pseudopotential (local + nonlocal), Hartree and exchange-correlation (including exact exchange if calculated).
 - **Default**: False
 
+### out_hr_npz
+
+- **Type**: Boolean
+- **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)*
+- **Description**: Whether to print Hamiltonian matrices H(R) in npz format. The output files are named output_HR0.npz, output_HR1.npz, and so on according to spin channel. This feature requires ABACUS to be built with CNPY.
+- **Default**: False
+- **Unit**: Ry
+
+### out_hsr_npz
+
+- **Type**: Boolean
+- **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)*
+- **Description**: Whether to print Hamiltonian matrices H(R) and overlap matrix S(R) in npz format. The output files are named output_SR.npz, output_HR0.npz, output_HR1.npz, and so on according to spin channel. This feature requires ABACUS to be built with CNPY.
+- **Default**: False
+
+### out_dm_npz
+
+- **Type**: Boolean
+- **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)*
+- **Description**: Whether to print density matrices DM(R) in npz format. The output files are named output_DM0.npz, output_DM1.npz, and so on according to spin channel. This feature requires ABACUS to be built with CNPY.
+- **Default**: False
+
 ### out_mul
 
 - **Type**: Boolean

diff --git a/docs/parameters.yaml b/docs/parameters.yaml
@@ -3147,6 +3147,30 @@ parameters:
     default_value: "False"
     unit: ""
     availability: Numerical atomic orbital basis
+  - name: out_hr_npz
+    category: Output information
+    type: Boolean
+    description: |
+      Whether to print Hamiltonian matrices H(R) in npz format. The output files are named output_HR0.npz, output_HR1.npz, and so on according to spin channel. This feature requires ABACUS to be built with CNPY.
+    default_value: "False"
+    unit: Ry
+    availability: Numerical atomic orbital basis (not gamma-only algorithm)
+  - name: out_hsr_npz
+    category: Output information
+    type: Boolean
+    description: |
+      Whether to print Hamiltonian matrices H(R) and overlap matrix S(R) in npz format. The output files are named output_SR.npz, output_HR0.npz, output_HR1.npz, and so on according to spin channel. This feature requires ABACUS to be built with CNPY.
+    default_value: "False"
+    unit: ""
+    availability: Numerical atomic orbital basis (not gamma-only algorithm)
+  - name: out_dm_npz
+    category: Output information
+    type: Boolean
+    description: |
+      Whether to print density matrices DM(R) in npz format. The output files are named output_DM0.npz, output_DM1.npz, and so on according to spin channel. This feature requires ABACUS to be built with CNPY.
+    default_value: "False"
+    unit: ""
+    availability: Numerical atomic orbital basis (not gamma-only algorithm)
   - name: out_mul
     category: Output information
     type: Boolean

diff --git a/source/source_io/module_ctrl/ctrl_scf_lcao.cpp b/source/source_io/module_ctrl/ctrl_scf_lcao.cpp
@@ -11,6 +11,7 @@
 #include "../module_unk/berryphase.h"                          // use berryphase
 #include "../module_hs/cal_pLpR.h"                            // use AngularMomentumCalculator()
 #include "source_io/module_hs/output_mat_sparse.h"                   // use ModuleIO::output_mat_sparse()
+#include "source_io/module_ml/io_npz.h"                       // use ModuleIO::output_mat_npz()
 #include "../module_hs/write_HS_R.h"                          // use ModuleIO::write_hsr()
 #include "../module_mulliken/output_mulliken.h"                     // use cal_mag()
 #include "../module_wannier/to_wannier90_lcao.h"                   // use toWannier90_LCAO
@@ -228,6 +229,35 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell,
                             out_app_flag, ucell.get_iat2iwt(), ucell.nat, istep);
     }
 
+    //------------------------------------------------------------------
+    //! 7a.1) Output H(R), S(R), and DM(R) matrices in NPZ format
+    //------------------------------------------------------------------
+    if (inp.out_hsr_npz)
+    {
+        std::string zipname = "output_SR.npz";
+        ModuleIO::output_mat_npz(ucell, zipname, *(p_hamilt->getSR()));
+    }
+
+    if (inp.out_hr_npz || inp.out_hsr_npz)
+    {
+        std::vector<hamilt::HContainer<TR>*> hr_vec = p_hamilt->getHR_vector();
+        for (int ispin = 0; ispin < hr_vec.size(); ++ispin)
+        {
+            std::string zipname = "output_HR" + std::to_string(ispin) + ".npz";
+            ModuleIO::output_mat_npz(ucell, zipname, *(hr_vec[ispin]));
+        }
+    }
+
+    if (inp.out_dm_npz)
+    {
+        const std::vector<hamilt::HContainer<double>*>& dmr_vec = dm->get_DMR_vector();
+        for (int ispin = 0; ispin < dmr_vec.size(); ++ispin)
+        {
+            std::string zipname = "output_DM" + std::to_string(ispin) + ".npz";
+            ModuleIO::output_mat_npz(ucell, zipname, *(dmr_vec[ispin]));
+        }
+    }
+
     //------------------------------------------------------------------
     //! 7b) Output dH, dS, T, r matrices (old sparse path, without H/S)
     //------------------------------------------------------------------

diff --git a/source/source_io/module_ml/io_npz.cpp b/source/source_io/module_ml/io_npz.cpp
@@ -321,7 +321,8 @@ void read_mat_npz(const Parallel_Orbitals* paraV,
 #endif
 }
 
-void output_mat_npz(const UnitCell& ucell, std::string& zipname, const hamilt::HContainer<double>& hR)
+template <typename T>
+void output_mat_npz_impl(const UnitCell& ucell, std::string& zipname, const hamilt::HContainer<T>& hR)
 {
     ModuleBase::TITLE("ModuleIO", "output_mat_npz");
 
@@ -412,13 +413,13 @@ void output_mat_npz(const UnitCell& ucell, std::string& zipname, const hamilt::H
 
 //fourth block: hr(i0,jR)
 #ifdef __MPI
-    hamilt::HContainer<double>* HR_serial;
+    hamilt::HContainer<T>* HR_serial;
     Parallel_Orbitals serialV;
     serialV.set_serial(PARAM.globalv.nlocal, PARAM.globalv.nlocal);
     serialV.set_atomic_trace(ucell.get_iat2iwt(), ucell.nat, PARAM.globalv.nlocal);
     if(GlobalV::MY_RANK == 0)
     {
-        HR_serial = new hamilt::HContainer<double>(&serialV);
+        HR_serial = new hamilt::HContainer<T>(&serialV);
     }
     hamilt::gatherParallels(hR, HR_serial, 0);
 
@@ -471,4 +472,16 @@ void output_mat_npz(const UnitCell& ucell, std::string& zipname, const hamilt::H
 #endif
 }
 
+void output_mat_npz(const UnitCell& ucell, std::string& zipname, const hamilt::HContainer<double>& hR)
+{
+    output_mat_npz_impl(ucell, zipname, hR);
+}
+
+void output_mat_npz(const UnitCell& ucell,
+                    std::string& zipname,
+                    const hamilt::HContainer<std::complex<double>>& hR)
+{
+    output_mat_npz_impl(ucell, zipname, hR);
+}
+
 } // namespace ModuleIO
diff --git a/source/source_io/module_ml/io_npz.h b/source/source_io/module_ml/io_npz.h
@@ -5,6 +5,7 @@
 #include "source_cell/unitcell.h"
 #include "source_lcao/module_hcontainer/hcontainer.h"
 
+#include <complex>
 #include <string>
 #include <vector>
 
@@ -17,6 +18,9 @@ void read_mat_npz(const Parallel_Orbitals* paraV,
                   hamilt::HContainer<double>& hR);
 
 void output_mat_npz(const UnitCell& ucell, std::string& zipname, const hamilt::HContainer<double>& hR);
+void output_mat_npz(const UnitCell& ucell,
+                    std::string& zipname,
+                    const hamilt::HContainer<std::complex<double>>& hR);
 
 } // namespace ModuleIO
 

diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h
@@ -400,6 +400,9 @@ struct Input_para
                                   ///< KS-orbital representation.
     std::vector<int> out_mat_xc2 = {0, 8};  ///< output Vxc(R) matrix with precision
     bool out_eband_terms = false; ///< output the band energy terms separately
+    bool out_hr_npz = false;      ///< output H(R) matrix in npz format
+    bool out_hsr_npz = false;     ///< output H(R) and S(R) matrices in npz format
+    bool out_dm_npz = false;      ///< output DM(R) matrix in npz format
     int out_interval = 1;
     bool out_app_flag = true; ///< whether output r(R), H(R), S(R), T(R), and dH(R) matrices
                               ///< in an append manner during MD liuyu 2023-03-20

diff --git a/source/source_io/module_parameter/read_input_item_output.cpp b/source/source_io/module_parameter/read_input_item_output.cpp
@@ -594,7 +594,7 @@ Also controled by out_freq_ion and out_app_flag.
         };
         item.check_value = [](const Input_Item& item, const Parameter& para) {
             if ((para.inp.out_mat_r[0] || para.inp.out_mat_hs2[0] || para.inp.out_mat_t[0] || para.inp.out_mat_dh[0]
-                 || para.inp.dm_to_rho)
+                 || para.inp.out_hr_npz || para.inp.out_hsr_npz || para.inp.out_dm_npz || para.inp.dm_to_rho)
                 && para.sys.gamma_only_local)
             {
                 ModuleBase::WARNING_QUIT("ReadInput",
@@ -816,6 +816,72 @@ The circle order of the charge density on real space grids is: x is the outer lo
         read_sync_bool(input.out_eband_terms);
         this->add_item(item);
     }
+    {
+        Input_Item item("out_hr_npz");
+        item.annotation = "output H(R) matrix in npz format";
+        item.category = "Output information";
+        item.type = "Boolean";
+        item.description = "Whether to print Hamiltonian matrices H(R) in npz format. This feature does not work for gamma-only calculations.";
+        item.default_value = "False";
+        item.unit = "Ry";
+        item.availability = "Numerical atomic orbital basis (not gamma-only algorithm)";
+        read_sync_bool(input.out_hr_npz);
+        item.check_value = [](const Input_Item& item, const Parameter& para) {
+            if (para.input.out_hr_npz)
+            {
+#ifndef __USECNPY
+                ModuleBase::WARNING_QUIT("ReadInput",
+                                         "to write in npz format, please "
+                                         "recompile with -DENABLE_CNPY=1");
+#endif
+            }
+        };
+        this->add_item(item);
+    }
+    {
+        Input_Item item("out_hsr_npz");
+        item.annotation = "output H(R) and S(R) matrices in npz format";
+        item.category = "Output information";
+        item.type = "Boolean";
+        item.description = "Whether to print Hamiltonian matrices H(R) and overlap matrix S(R) in npz format. This feature does not work for gamma-only calculations.";
+        item.default_value = "False";
+        item.unit = "Ry";
+        item.availability = "Numerical atomic orbital basis (not gamma-only algorithm)";
+        read_sync_bool(input.out_hsr_npz);
+        item.check_value = [](const Input_Item& item, const Parameter& para) {
+            if (para.input.out_hsr_npz)
+            {
+#ifndef __USECNPY
+                ModuleBase::WARNING_QUIT("ReadInput",
+                                         "to write in npz format, please "
+                                         "recompile with -DENABLE_CNPY=1");
+#endif
+            }
+        };
+        this->add_item(item);
+    }
+    {
+        Input_Item item("out_dm_npz");
+        item.annotation = "output DM(R) matrix in npz format";
+        item.category = "Output information";
+        item.type = "Boolean";
+        item.description = "Whether to print density matrices DM(R) in npz format. This feature does not work for gamma-only calculations.";
+        item.default_value = "False";
+        item.unit = "";
+        item.availability = "Numerical atomic orbital basis (not gamma-only algorithm)";
+        read_sync_bool(input.out_dm_npz);
+        item.check_value = [](const Input_Item& item, const Parameter& para) {
+            if (para.input.out_dm_npz)
+            {
+#ifndef __USECNPY
+                ModuleBase::WARNING_QUIT("ReadInput",
+                                         "to write in npz format, please "
+                                         "recompile with -DENABLE_CNPY=1");
+#endif
+            }
+        };
+        this->add_item(item);
+    }
     {
         Input_Item item("out_mul");
         item.annotation = "mulliken charge or not";

diff --git a/source/source_io/test_serial/read_input_item_test.cpp b/source/source_io/test_serial/read_input_item_test.cpp
@@ -973,6 +973,24 @@ TEST_F(InputTest, Item_test)
         it->second.reset_value(it->second, param);
         EXPECT_EQ(param.input.out_mat_hs[0], 1);
     }
+    { // out_hr_npz
+        auto it = find_label("out_hr_npz", readinput.input_lists);
+        it->second.str_values = {"1"};
+        it->second.read_value(it->second, param);
+        EXPECT_EQ(param.input.out_hr_npz, true);
+    }
+    { // out_hsr_npz
+        auto it = find_label("out_hsr_npz", readinput.input_lists);
+        it->second.str_values = {"1"};
+        it->second.read_value(it->second, param);
+        EXPECT_EQ(param.input.out_hsr_npz, true);
+    }
+    { // out_dm_npz
+        auto it = find_label("out_dm_npz", readinput.input_lists);
+        it->second.str_values = {"1"};
+        it->second.read_value(it->second, param);
+        EXPECT_EQ(param.input.out_dm_npz, true);
+    }
 }
 TEST_F(InputTest, Item_test2)
 {