Merge branch 'hr' into develop

Author: mohanchen

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

@@ -1773,8 +1773,27 @@ These variables are used to control the output of properties.
 
 - **Type**: Boolean \[Integer\](optional)
 - **Availability**: Numerical atomic orbital basis
-- **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).
+- **Description**: Whether to print the upper triangular part of the Hamiltonian matrices 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).
+  - For gamma only case:
+    - nspin = 1: `hks1_nao.txt` for the Hamiltonian matrix and `sks1_nao.txt` for the overlap matrix;
+    - nspin = 2: `hks1_nao.txt` and `hks2_nao.txt` for the Hamiltonian matrix and `sks1_nao.txt` for the overlap matrix. Note that the code will not output `sks2_nao.txt` because it is the same as `sks1_nao.txt`;
+    - nspin = 4: `hks12_nao.txt` for the Hamiltonian matrix and `sks12_nao.txt` for the overlap matrix.
+  - For multi-k points case:
+    - nspin = 1: `hks1k1_nao.txt` for the Hamiltonian matrix at the 1st k-point, and `sks1k1_nao.txt` for the overlap matrix for the 1st k-point, ...;
+    - nspin = 2: `hks1k1_nao.txt` and `hks2k1_nao.txt` for the two spin channels of the Hamiltonian matrix at the 1st k-point, and `sks1k1_nao.txt` for the overlap matrix for the 1st k-point. Note that the code will not output `sks2k1_nao.txt` because it is the same as `sks1k1_nao.txt`, ...;
+    - nspin = 4: `hks12k1_nao.txt` for the Hamiltonian matrix at the 1st k-point, and `sks12k1_nao.txt` for the overlap matrix for the 1st k-point, ...;
 - **Default**: False 8
+- **Unit**: Ry
+- **NOTE**: In the 3.10-LTS version, the file names are `data-0-H` and `data-0-S`, etc.
+
+### out_mat_hs2
+
+- **Type**: Boolean
+- **Availability**: Numerical atomic orbital basis (not gamma-only algorithm)
+- **Description**: Whether to print files containing the Hamiltonian matrix $H(R)$ (in Ry) and overlap matrix $S(R)$ into files in the directory `OUT.${suffix}`. For more information, please refer to [hs_matrix.md](../elec_properties/hs_matrix.md#out_mat_hs2).
+- **Default**: False
+- **Unit**: Ry
+- **NOTE**: In the 3.10-LTS version, the file names are `data-HR-sparse_SPIN0.csr` and `data-SR-sparse_SPIN0.csr`, etc.
 
 ### out_mat_tk
 
@@ -1790,13 +1809,6 @@ 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_hs2
-
-- **Type**: Boolean
-- **Availability**: Numerical atomic orbital basis (not gamma-only algorithm)
-- **Description**: Whether to print files containing the Hamiltonian matrix $H(R)$ (in Ry) and overlap matrix $S(R)$ into files in the directory `OUT.${suffix}`. For more information, please refer to [hs_matrix.md](../elec_properties/hs_matrix.md#out_mat_hs2).
-- **Default**: False
-
 ### out_mat_t
 
 - **Type**: Boolean
@@ -1822,7 +1834,7 @@ These variables are used to control the output of properties.
 
 - **Type**: Boolean
 - **Availability**: Numerical atomic orbital (NAO) and NAO-in-PW basis
-- **Description**: Whether to print the upper triangular part of the exchange-correlation matrices in **Kohn-Sham orbital representation** (unit: Ry): $\braket{\psi_i|V_\text{xc}^\text{(semi-)local}+V_\text{exx}+V_\text{DFTU}|\psi_j}$ for each k point into files in the directory `OUT.${suffix}`, which is useful for the subsequent GW calculation. (Note that currently DeePKS term is not included. ) The files are named `k-$k-Vxc`, the meaning of `$k`corresponding to k point and spin  is same as [hs_matrix.md](../elec_properties/hs_matrix.md#out_mat_hs).
+- **Description**: Whether to print the upper triangular part of the exchange-correlation matrices in **Kohn-Sham orbital representation** (unit: Ry): $\braket{\psi_i|V_\text{xc}^\text{(semi-)local}+V_\text{exx}+V_\text{DFTU}|\psi_j}$ for each k point into files in the directory `OUT.${suffix}`, which is useful for the subsequent GW calculation. (Note that currently DeePKS term is not included. ) The files are named `k-$k-Vxc`, the meaning of `$k`corresponding to k point and spin is the same as [hs_matrix.md](../elec_properties/hs_matrix.md#out_mat_hs).
 The band (KS orbital) energy for each (k-point, spin, band) will be printed in the file `OUT.${suffix}/vxc_out.dat`. If EXX is calculated, the local and EXX part of band energy will also be printed in `OUT.${suffix}/vxc_local_out.dat`and `OUT.${suffix}/vxc_exx_out.dat`, respectively. All the `vxc*_out.dat` files contains 3 integers (nk, nspin, nband) followed by nk\*nspin\*nband lines of energy Hartree and eV.
 - **Default**: False
 

source/module_esolver/lcao_after_scf.cpp

@@ -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_io/cal_r_overlap_R.cpp

@@ -459,7 +459,7 @@ void cal_r_overlap_R::out_rR(const UnitCell& ucell, const Grid_Driver& gd, const
             ssr << PARAM.globalv.global_out_dir << "data-rR-sparse.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)
@@ -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/filename.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/write_HS.h

@@ -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