Merge branch 'develop' into ml4dft

Author: haozhihan

docs/advanced/elec_properties/hs_matrix.md

@@ -12,7 +12,7 @@ and
 
 ## out_mat_hs
 
-Users may set the keyword [out_mat_hs](../input_files/input-main.md#out_mat_hs) to true for outputting the upper triangular part of the Hamiltonian matrices and overlap matrices for each k point into files in the directory `OUT.${suffix}`. It is available for both gamma_only and multi-k calculations. 
+Users can set the keyword [out_mat_hs](../input_files/input-main.md#out_mat_hs) to true for outputting the upper triangular part of the Hamiltonian matrices and overlap matrices for each k point into files in the directory `OUT.${suffix}`. It is available for both gamma_only and multi-k calculations. 
 
 The files are named `data-$k-H` and `data-$k-S`, where `$k` is a composite index consisting of the k point index as well as the spin index. The corresponding sequence of the orbitals can be seen in [Basis Set](../pp_orb.md#basis-set).
 
@@ -35,7 +35,7 @@ The rest of the file contains the upper triangular part of the specified matrice
 
 The output of R-space 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 `data-HR-sparse_SPIN0.csr` and `data-SR-sparse_SPIN0.csr` are generated, which contain the Hamiltonian matrix $H(R)$ and overlap matrix $S(R)$ respectively. For nspin = 2, three files `data-HR-sparse_SPIN0.csr` and `data-HR-sparse_SPIN1.csr` and `data-SR-sparse_SPIN0.csr` are created, where the first two contain $H(R)$ for spin up and spin down, respectively.
+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.
 
 As for molecular dynamics calculations, the format is controlled by [out_interval](../input_files/input-main.md#out_interval) and [out_app_flag](../input_files/input-main.md#out_app_flag) in the same manner as the position matrix as detailed in [out_mat_r](../input_files/input-main.md#out_mat_r).
 
@@ -74,4 +74,4 @@ We provide [examples](https://github.com/deepmodeling/abacus-develop/tree/develo
 - out_hs_multik : writing H(k) and S(k) for multi-k calculation
 - out_s_multik : running get_S for multi-k calculation
 
-Reference output files are provided in each directory.
+Reference output files are provided in each directory.

docs/advanced/input_files/input-main.md

@@ -126,51 +126,26 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep, const
 #endif
 
     //------------------------------------------------------------------
-    // 5) write Hamiltonian and Overlap matrix in LCAO basis
-    //------------------------------------------------------------------
-    for (int ik = 0; ik < this->kv.get_nks(); ++ik)
-    {
-        if (PARAM.inp.out_mat_hs[0])
-        {
-            this->p_hamilt->updateHk(ik);
-        }
-        bool bit = false; // LiuXh, 2017-03-21
-        // if set bit = true, there would be error in soc-multi-core
-        // calculation, noted by zhengdy-soc
-        if (this->psi != nullptr && (istep % PARAM.inp.out_interval == 0))
-        {
-            hamilt::MatrixBlock<TK> h_mat;
-            hamilt::MatrixBlock<TK> s_mat;
-
-            this->p_hamilt->matrix(h_mat, s_mat);
-
-            if (PARAM.inp.out_mat_hs[0])
-            {
-                ModuleIO::save_mat(istep,
-                                   h_mat.p,
-                                   PARAM.globalv.nlocal,
-                                   bit,
-                                   PARAM.inp.out_mat_hs[1],
-                                   1,
-                                   PARAM.inp.out_app_flag,
-                                   "H",
-                                   "data-" + std::to_string(ik),
-                                   this->pv,
-                                   GlobalV::DRANK);
-                ModuleIO::save_mat(istep,
-                                   s_mat.p,
-                                   PARAM.globalv.nlocal,
-                                   bit,
-                                   PARAM.inp.out_mat_hs[1],
-                                   1,
-                                   PARAM.inp.out_app_flag,
-                                   "S",
-                                   "data-" + std::to_string(ik),
-                                   this->pv,
-                                   GlobalV::DRANK);
-            }
-        }
-    }
+    // 5) write H(k) and S(k) in NAO basis set 
+    //------------------------------------------------------------------
+	if (PARAM.inp.out_mat_hs[0])
+	{
+		if (this->psi != nullptr && (istep % PARAM.inp.out_interval == 0))
+		{
+			ModuleIO::write_hsk(PARAM.globalv.global_out_dir,
+					PARAM.inp.nspin,
+					this->kv.get_nks(), 
+					this->kv.get_nkstot(), 
+					this->kv.ik2iktot, 
+                    this->kv.isk,
+					this->p_hamilt, 
+					this->pv, 
+					PARAM.globalv.gamma_only_local,
+					PARAM.inp.out_app_flag,
+					istep,
+					GlobalV::ofs_running);
+		}
+	}
 
     //------------------------------------------------------------------
     // 6) write electronic wavefunctions in LCAO basis
@@ -182,9 +157,9 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep, const
 				this->psi[0],
 				this->pelec->ekb,
 				this->pelec->wg,
-				this->pelec->klist->kvec_c,
-				this->pelec->klist->ik2iktot,
-				this->pelec->klist->get_nkstot(),
+				this->kv.kvec_c,
+				this->kv.ik2iktot,
+				this->kv.get_nkstot(),
 				this->pv,
 				PARAM.inp.nspin,
 				istep);
@@ -201,11 +176,11 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep, const
         LCAO_Deepks_Interface<TK, TR> deepks_interface(ld_shared_ptr);
 
         deepks_interface.out_deepks_labels(this->pelec->f_en.etot,
-                                           this->pelec->klist->get_nks(),
+                                           this->kv.get_nks(),
                                            ucell.nat,
                                            PARAM.globalv.nlocal,
                                            this->pelec->ekb,
-                                           this->pelec->klist->kvec_d,
+                                           this->kv.kvec_d,
                                            ucell,
                                            orb_,
                                            this->gd,
@@ -370,15 +345,23 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep, const
         for (int ik = 0; ik < this->kv.get_nks() / nspin_k; ++ik)
         {
             ekinetic->init(ik);
-            ModuleIO::save_mat(0,
+
+            const int out_label = 1; // 1: .txt, 2: .dat
+
+			std::string t_fn = ModuleIO::filename_output(PARAM.globalv.global_out_dir,
+					"tk","nao",ik,this->kv.ik2iktot,
+					PARAM.inp.nspin,this->kv.get_nkstot(),
+					out_label,PARAM.inp.out_app_flag,
+                    PARAM.globalv.gamma_only_local,istep);
+
+            ModuleIO::save_mat(istep,
                                hsk.get_hk(),
                                PARAM.globalv.nlocal,
-                               false,
+                               false, // bit
                                PARAM.inp.out_mat_tk[1],
-                               1,
+                               1, // true for upper triangle matrix
                                PARAM.inp.out_app_flag,
-                               "T",
-                               "data-" + std::to_string(ik),
+                               t_fn, 
                                this->pv,
                                GlobalV::DRANK);
         }

source/module_esolver/lcao_after_scf.cpp

@@ -274,7 +274,7 @@ void cal_r_overlap_R::out_rR(const UnitCell& ucell, const Grid_Driver& gd, const
     int output_R_number = 0;
 
     std::stringstream tem1;
-    tem1 << PARAM.globalv.global_out_dir << "temp-data-rR-sparse.dat";
+    tem1 << PARAM.globalv.global_out_dir << "tmp-rr.csr";
     std::ofstream ofs_tem1;
     std::ifstream ifs_tem1;
 
@@ -451,15 +451,15 @@ void cal_r_overlap_R::out_rR(const UnitCell& ucell, const Grid_Driver& gd, const
         std::stringstream ssr;
         if (PARAM.inp.calculation == "md" && !PARAM.inp.out_app_flag)
         {
-            ssr << PARAM.globalv.global_matrix_dir << step << "_"
-                << "data-rR-sparse.csr";
+            ssr << PARAM.globalv.global_matrix_dir
+                << "rrg" << step << ".csr";
         }
         else
         {
-            ssr << PARAM.globalv.global_out_dir << "data-rR-sparse.csr";
+            ssr << PARAM.globalv.global_out_dir << "rr.csr";
         }
 
-        if (binary)
+        if (binary) // .dat
         {
             ofs_tem1.close();
             int nlocal = PARAM.globalv.nlocal;
@@ -480,7 +480,7 @@ void cal_r_overlap_R::out_rR(const UnitCell& ucell, const Grid_Driver& gd, const
             ifs_tem1.close();
             out_r.close();
         }
-        else
+        else // .txt
         {
             ofs_tem1.close();
             if (PARAM.inp.calculation == "md" && PARAM.inp.out_app_flag && step)
@@ -524,12 +524,12 @@ void cal_r_overlap_R::out_rR_other(const UnitCell& ucell, const int& istep, cons
     std::stringstream ssr;
     if (PARAM.inp.calculation == "md" && !PARAM.inp.out_app_flag)
     {
-        ssr << PARAM.globalv.global_matrix_dir << step << "_"
-            << "data-rR-sparse.csr";
+        ssr << PARAM.globalv.global_matrix_dir
+            << "rrg" << step << ".csr";
     }
     else
     {
-        ssr << PARAM.globalv.global_out_dir << "data-rR-sparse.csr";
+        ssr << PARAM.globalv.global_out_dir << "rr.csr";
     }
 
     if (GlobalV::DRANK == 0)
@@ -676,13 +676,13 @@ void cal_r_overlap_R::out_rR_other(const UnitCell& ucell, const int& istep, cons
 
         Parallel_Reduce::reduce_all(rR_nonzero_num, 3);
 
-        if (binary)
+        if (binary) // .dat
         {
             out_r.write(reinterpret_cast<char*>(&dRx), sizeof(int));
             out_r.write(reinterpret_cast<char*>(&dRy), sizeof(int));
             out_r.write(reinterpret_cast<char*>(&dRz), sizeof(int));
         }
-        else
+        else // .txt
         {
             out_r << dRx << " " << dRy << " " << dRz << std::endl;
         }

source/module_io/cal_r_overlap_R.cpp

@@ -22,16 +22,16 @@ std::string filename_output(
     // {k(optinal)}{k-point index}{g(optional)}{geometry index1}{_basis(nao|pw)} 
     // + {".txt"/".dat"}"
 
-	std::set<std::string> valid_properties = {"wf", "chg", "h", "s"};
+	std::set<std::string> valid_properties = {"wf", "chg", "hk", "sk", "tk", "vxc"};
 	if (valid_properties.find(property) == valid_properties.end()) 
 	{
-		ModuleBase::WARNING_QUIT("ModuleIO::filename_output", "unknown property");
+		ModuleBase::WARNING_QUIT("ModuleIO::filename_output", "unknown property in filename function");
 	}
 
 	std::set<std::string> valid_basis = {"pw", "nao"};
 	if (valid_basis.find(basis) == valid_basis.end()) 
 	{
-		ModuleBase::WARNING_QUIT("ModuleIO::filename_output", "unknown basis");
+		ModuleBase::WARNING_QUIT("ModuleIO::filename_output", "unknown basis in filename function");
 	}
 
     assert(ik_local>=0);

source/module_io/filename.cpp

@@ -1,21 +1,38 @@
 #ifndef WRITE_HS_H
 #define WRITE_HS_H
 
-#include "module_base/global_function.h"
-#include "module_base/global_variable.h"
-
 #include <string>
+#include <vector>
+
+//#include "module_base/global_function.h"
+//#include "module_base/global_variable.h"
+#include "module_basis/module_ao/parallel_orbitals.h" // use Parallel_Orbitals
+
 
 // mohan add this file 2010-09-10
 namespace ModuleIO
 {
-    /// @brief save a square matrix
+	template<typename T>
+		void write_hsk(
+				const std::string &global_out_dir,
+				const int nspin,
+				const int nks, 
+				const int nkstot, 
+				const std::vector<int> &ik2iktot,
+				const std::vector<int> &isk,
+				hamilt::Hamilt<T>* p_hamilt,
+				const Parallel_Orbitals &pv,
+				const bool gamma_only,
+				const bool out_app_flag,
+				const int istep,
+				std::ofstream &ofs_running);	
+
+    /// @brief save a square matrix, such as H(k) and S(k)
     /// @param[in] istep : the step of the calculation
     /// @param[in] mat : the local matrix
     /// @param[in] bit : true for binary, false for decimal
     /// @param[in] tri : true for upper triangle, false for full matrix
     /// @param[in] app : true for append, false for overwrite
-    /// @param[in] label : the symbol of the matrix, like "H", "S"
     /// @param[in] file_name : the name of the output file
     /// @param[in] pv : the 2d-block parallelization information
     /// @param[in] drank : the rank of the current process
@@ -27,16 +44,11 @@ namespace ModuleIO
         const int precision,
         const bool tri,
         const bool app,
-        const std::string label,
         const std::string& file_name,
         const Parallel_2D& pv,
         const int drank,
         const bool reduce = true);
 
-
-// mohan comment out 2021-02-10
-// void save_HS_ccf(const int &iter, const int &Hnnz, const int *colptr_H, const int *rowind_H, 
-// const double *nzval_H, const double *nzval_S, bool bit);
 }
 #include "write_HS.hpp"
 #endif