Feature: output T(k) (#4982)

* output Tk * change header order
deepmodeling · Aug 24, 2024 · 2890890 · 2890890
1 parent 75b5cf7
commit 2890890
Show file tree

Hide file tree

Showing 4 changed files with 67 additions and 2 deletions.
diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
@@ -143,6 +143,7 @@
     - [out\_level](#out_level)
     - [out\_alllog](#out_alllog)
     - [out\_mat\_hs](#out_mat_hs)
+    - [out\_mat\_tk](#out_mat_tk)
     - [out\_mat\_r](#out_mat_r)
     - [out\_mat\_hs2](#out_mat_hs2)
     - [out\_mat\_t](#out_mat_t)
@@ -1667,6 +1668,13 @@ These variables are used to control the output of properties.
 - **Description**: Whether to print the upper triangular part of the Hamiltonian matrices (in Ry) and overlap matrices for each k point into files in the directory `OUT.${suffix}`. The second number controls precision. For more information, please refer to [hs_matrix.md](../elec_properties/hs_matrix.md#out_mat_hs). Also controled by [out_interval](#out_interval) and [out_app_flag](#out_app_flag).
 - **Default**: False 8
 
+### out_mat_tk
+
+- **Type**: Boolean \[Integer\](optional)
+- **Availability**: Numerical atomic orbital basis
+- **Description**: Whether to print the upper triangular part of the kinetic matrices (in Ry) for each k point into `OUT.${suffix}/data-i-T`, where i is the index of k points (see `OUT.${suffix}/kpoints`). One may optionally provide a second parameter to specify the precision. 
+- **Default**: False \[8\]
+
 ### out_mat_r
 
 - **Type**: Boolean

diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp
@@ -26,6 +26,7 @@
 #include "module_elecstate/module_charge/symmetry_rho.h"
 #include "module_elecstate/occupy.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" // need divide_HS_in_frag
+#include "module_hamilt_lcao/hamilt_lcaodft/hs_matrix_k.hpp"
 #include "module_hamilt_lcao/module_dftu/dftu.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "module_io/print_info.h"
@@ -1254,6 +1255,38 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
                        GlobalV::NPROC);
         tqo.calculate();
     }
+
+    if (PARAM.inp.out_mat_tk[0])
+    {
+        hamilt::HS_Matrix_K<TK> hsk(&pv, true);
+        hamilt::HContainer<TR> hR(&pv);
+        hamilt::Operator<TK>* ekinetic = 
+            new hamilt::EkineticNew<hamilt::OperatorLCAO<TK, TR>>(&hsk,
+                                                                  this->kv.kvec_d,
+                                                                  &hR,
+                                                                  &GlobalC::ucell,
+                                                                  &GlobalC::GridD,
+                                                                  two_center_bundle_.kinetic_orb.get());
+
+        const int nspin_k = (GlobalV::NSPIN == 2 ? 2 : 1);
+        for (int ik = 0; ik < this->kv.get_nks() / nspin_k; ++ik)
+        {
+            ekinetic->init(ik);
+            ModuleIO::save_mat(0,
+                               hsk.get_hk(),
+                               GlobalV::NLOCAL,
+                               false,
+                               PARAM.inp.out_mat_tk[1],
+                               1,
+                               GlobalV::out_app_flag,
+                               "T",
+                               "data-" + std::to_string(ik),
+                               this->pv,
+                               GlobalV::DRANK);
+        }
+
+        delete ekinetic;
+    }
 }
 
 //------------------------------------------------------------------------------

diff --git a/source/module_io/read_input_item_output.cpp b/source/module_io/read_input_item_output.cpp
@@ -254,6 +254,29 @@ void ReadInput::item_output()
         sync_intvec(input.out_mat_hs, 2, 0);
         this->add_item(item);
     }
+    {
+        Input_Item item("out_mat_tk");
+        item.annotation = "output T(k)";
+        item.read_value = [](const Input_Item& item, Parameter& para) {
+            size_t count = item.get_size();
+            if (count == 1)
+            {
+                para.input.out_mat_tk[0] = std::stoi(item.str_values[0]);
+                para.input.out_mat_tk[1] = 8;
+            }
+            else if (count == 2)
+            {
+                para.input.out_mat_tk[0] = std::stoi(item.str_values[0]);
+                para.input.out_mat_tk[1] = std::stoi(item.str_values[1]);
+            }
+            else
+            {
+                ModuleBase::WARNING_QUIT("ReadInput", "out_mat_tk should have 1 or 2 values");
+            }
+        };
+        sync_intvec(input.out_mat_tk, 2, 0);
+        this->add_item(item);
+    }
     {
         Input_Item item("out_mat_hs2");
         item.annotation = "output H(R) and S(R) matrix";
@@ -453,4 +476,4 @@ void ReadInput::item_output()
         this->add_item(item);
     }
 }
-} // namespace ModuleIO
+} // namespace ModuleIO
diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h
@@ -324,6 +324,7 @@ struct Input_para
     bool out_dm1 = false;                 ///< output density matrix. (multi-k points)
     bool out_bandgap = false;             ///< QO added for bandgap printing
     std::vector<int> out_mat_hs = {0, 8}; ///< output H matrix and S matrix in local basis.
+    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_dh = false;
@@ -567,4 +568,4 @@ struct Input_para
     bool test_stress = false;     ///< test the stress.
     bool test_skip_ewald = false; ///< variables for test only
 };
-#endif
+#endif