deepmodeling · mohanchen · Jul 2, 2026 · Jun 16, 2026 · Jun 29, 2026 · Jun 30, 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/cal_mag.h"                          // use cal_mag()
 #include "../module_wannier/to_wannier90_lcao.h"                   // use toWannier90_LCAO
@@ -227,6 +228,37 @@ 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 = PARAM.globalv.global_out_dir + "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
+                = PARAM.globalv.global_out_dir + "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
+                = PARAM.globalv.global_out_dir + "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
@@ -401,6 +401,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)
 {

diff --git a/tests/03_NAO_multik/CASES_CNPY.txt b/tests/03_NAO_multik/CASES_CNPY.txt
@@ -0,0 +1,3 @@
+scf_out_hr_npz
+scf_out_hsr_npz
+scf_out_dm_npz
diff --git a/tests/03_NAO_multik/CMakeLists.txt b/tests/03_NAO_multik/CMakeLists.txt
@@ -14,3 +14,11 @@ else()
         WORKING_DIRECTORY ${ABACUS_TEST_DIR}/03_NAO_multik
     )
 endif()
+
+if(ENABLE_CNPY)
+    add_test(
+        NAME 03_NAO_multik_npz
+        COMMAND ${BASH} ../integrate/Autotest.sh -a ${ABACUS_BIN_PATH} -n 4 -f CASES_CNPY.txt
+        WORKING_DIRECTORY ${ABACUS_TEST_DIR}/03_NAO_multik
+    )
+endif()
diff --git a/tests/03_NAO_multik/scf_out_dm_npz/INPUT b/tests/03_NAO_multik/scf_out_dm_npz/INPUT
@@ -0,0 +1,30 @@
+INPUT_PARAMETERS
+#Parameters (1.General)
+suffix            autotest
+calculation       scf
+
+nbands            6
+symmetry          0
+pseudo_dir        ../../PP_ORB
+orbital_dir       ../../PP_ORB
+gamma_only        0
+
+#Parameters (2.Iteration)
+ecutwfc           20
+scf_thr           1e-8
+scf_nmax          1
+
+#Parameters (3.Basis)
+basis_type        lcao
+
+#Parameters (4.Smearing)
+smearing_method   gauss
+smearing_sigma    0.002
+
+#Parameters (5.Mixing)
+mixing_type       broyden
+mixing_beta       0.7
+mixing_gg0        0.0
+
+out_dm_npz        1
+ks_solver         scalapack_gvx
diff --git a/tests/03_NAO_multik/scf_out_dm_npz/KPT b/tests/03_NAO_multik/scf_out_dm_npz/KPT
@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+2 1 1 0 0 0
diff --git a/tests/03_NAO_multik/scf_out_dm_npz/README b/tests/03_NAO_multik/scf_out_dm_npz/README
@@ -0,0 +1 @@
+test the output of DM(R) matrix in NPZ format under OUT.autotest
diff --git a/tests/03_NAO_multik/scf_out_dm_npz/STRU b/tests/03_NAO_multik/scf_out_dm_npz/STRU
@@ -0,0 +1,22 @@
+ATOMIC_SPECIES
+Si 14 Si_dojo_nsoc.upf upf201
+
+NUMERICAL_ORBITAL
+Si_dojo_6au_sz.orb
+
+LATTICE_CONSTANT
+15.3  // add lattice constant
+
+LATTICE_VECTORS
+0.0 0.5 0.5
+0.5 0.0 0.5
+0.5 0.5 0.0
+
+ATOMIC_POSITIONS
+Direct
+
+Si	// Element type
+0.0	// magnetism
+2
+0.00 0.00 0.00 1 1 1
+0.25 0.25 0.25 1 1 1
diff --git a/tests/03_NAO_multik/scf_out_dm_npz/result.ref b/tests/03_NAO_multik/scf_out_dm_npz/result.ref
@@ -0,0 +1,4 @@
+etotref -174.1078103590385
+etotperatomref -87.0539051795
+OutputDMNPZ_pass 0
+totaltimeref 0.25
Original file line number	Diff line number	Diff line change
		@@ -0,0 +1 @@
		test the output of DM(R) matrix in NPZ format under OUT.autotest