Feature: add out_mat_hr0 to output H0 for first iteration

Author: dyzheng

docs/advanced/input_files/input-main.md

@@ -150,6 +150,7 @@
     - [out\_mat\_hs](#out_mat_hs)
     - [out\_mat\_tk](#out_mat_tk)
     - [out\_mat\_r](#out_mat_r)
+    - [out\_mat\_hr0](#out_mat_hr0)
     - [out\_mat\_hs2](#out_mat_hs2)
     - [out\_mat\_t](#out_mat_t)
     - [out\_mat\_dh](#out_mat_dh)
@@ -1756,6 +1757,13 @@ These variables are used to control the output of properties.
 - **Description**: Whether to print the matrix representation of the position matrix (in Bohr) into a file named `data-rR-tr` in the directory `OUT.${suffix}`. If [calculation](#calculation) is set to `get_S`, the position matrix can be obtained without scf iterations. For more information, please refer to [position_matrix.md](../elec_properties/position_matrix.md#extracting-position-matrices).
 - **Default**: False
 
+### out_mat_hr0
+
+- **Type**: Boolean
+- **Availability**: Numerical atomic orbital basis (not gamma-only algorithm)
+- **Description**: Whether to print files containing the initial Hamiltonian matrix $H(R)$ (in Ry) into files `data-HR0.dat` in the directory `OUT.${suffix}`. 
+- **Default**: False
+
 ### out_mat_hs2
 
 - **Type**: Boolean

source/module_esolver/esolver_ks_lcao.cpp

@@ -870,6 +870,26 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
         this->p_chgmix->mix_dmr(dm);
     }
 
+    // output the first H(R) when out_mat_hr0 is true
+    if(PARAM.inp.out_mat_hr0 == true && istep == 0 && iter == 1)
+    {
+        ModuleIO::output_mat_sparse(false,
+                                    PARAM.inp.out_mat_hr0,
+                                    false,
+                                    false,
+                                    false,
+                                    istep,
+                                    this->pelec->pot->get_effective_v(),
+                                    this->pv,
+                                    this->GK,
+                                    two_center_bundle_,
+                                    orb_,
+                                    ucell,
+                                    this->gd,
+                                    this->kv,
+                                    this->p_hamilt);
+    }
+
     // 2) save charge density
     // Peize Lin add 2020.04.04
     if (GlobalC::restart.info_save.save_charge)
@@ -1159,6 +1179,7 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep)
     {
         //! Print out sparse matrix
         ModuleIO::output_mat_sparse(PARAM.inp.out_mat_hs2,
+                                    false,
                                     PARAM.inp.out_mat_dh,
                                     PARAM.inp.out_mat_t,
                                     PARAM.inp.out_mat_r,

source/module_io/output_mat_sparse.cpp

@@ -8,6 +8,7 @@ namespace ModuleIO
 
 template <>
 void output_mat_sparse(const bool& out_mat_hsR,
+                       const bool& out_mat_hr0,
                        const bool& out_mat_dh,
                        const bool& out_mat_t,
                        const bool& out_mat_r,
@@ -26,6 +27,7 @@ void output_mat_sparse(const bool& out_mat_hsR,
 
 template <>
 void output_mat_sparse(const bool& out_mat_hsR,
+                       const bool& out_mat_hr0,
                        const bool& out_mat_dh,
                        const bool& out_mat_t,
                        const bool& out_mat_r,
@@ -47,6 +49,27 @@ void output_mat_sparse(const bool& out_mat_hsR,
     {
         output_HSR(ucell,istep, v_eff, pv, HS_Arrays, grid, kv, p_ham);
     }
+    else if (out_mat_hr0)
+    {
+        output_HSR(
+            ucell,
+            istep, 
+            v_eff, 
+            pv, 
+            HS_Arrays, 
+            grid, 
+            kv, 
+            p_ham,
+#ifdef _EXX
+            nullptr,
+            nullptr,
+#endif 
+            "data-SR-sparse_SPIN0.csr", 
+            "data-HR0_SPIN0.csr", 
+            "data-HR0_SPIN1.csr"
+        );
+    }
+    
 
     //! generate a file containing the kinetic energy matrix
     if (out_mat_t)

source/module_io/output_mat_sparse.h

@@ -12,6 +12,7 @@ namespace ModuleIO
 /// @brief the output interface to write the sparse matrix of H, S, T, and r
 template <typename T>
 void output_mat_sparse(const bool& out_mat_hsR,
+                       const bool& out_mat_hr0,
                        const bool& out_mat_dh,
                        const bool& out_mat_t,
                        const bool& out_mat_r,

source/module_io/read_input_item_output.cpp

@@ -267,14 +267,26 @@ void ReadInput::item_output()
         sync_intvec(input.out_mat_tk, 2, 0);
         this->add_item(item);
     }
+    {
+        Input_Item item("out_mat_hr0");
+        item.annotation = "output initial H(R) matrix";
+        read_sync_bool(input.out_mat_hr0);
+        item.check_value = [](const Input_Item& item, const Parameter& para) {
+            if (para.input.out_mat_hr0 && para.sys.gamma_only_local)
+            {
+                ModuleBase::WARNING_QUIT("ReadInput", "out_mat_r is not available for gamma only calculations");
+            }
+        };
+        this->add_item(item);
+    }
     {
         Input_Item item("out_mat_hs2");
         item.annotation = "output H(R) and S(R) matrix";
         read_sync_bool(input.out_mat_hs2);
         item.check_value = [](const Input_Item& item, const Parameter& para) {
-            if (para.input.out_mat_r && para.sys.gamma_only_local)
+            if (para.input.out_mat_hs2 && para.sys.gamma_only_local)
             {
-                ModuleBase::WARNING_QUIT("ReadInput", "out_mat_r is not available for gamma only calculations");
+                ModuleBase::WARNING_QUIT("ReadInput", "out_mat_hs2 is not available for gamma only calculations");
             }
         };
         this->add_item(item);

source/module_parameter/input_parameter.h

@@ -337,6 +337,7 @@ struct Input_para
     std::vector<int> out_mat_tk = {0, 8}; ///< output T(k) matrix in local basis.
     bool out_mat_hs2 = false;             ///< LiuXh add 2019-07-16, output H(R) matrix and
                                           ///< S(R) matrix in local basis.
+    bool out_mat_hr0 = false;             ///< output H(R) with atomic charge
     bool out_mat_dh = false;
     bool out_mat_xc = false;      ///< output exchange-correlation matrix in
                                   ///< KS-orbital representation.