Merge branch 'develop' into fft6

Author: A-006

docs/advanced/input_files/input-main.md

@@ -989,7 +989,7 @@ calculations.
 
 - **Type**: String
 - **Description**: In our package, the XC functional can either be set explicitly using the `dft_functional` keyword in `INPUT` file. If `dft_functional` is not specified, ABACUS will use the xc functional indicated in the pseudopotential file.
-  On the other hand, if dft_functional is specified, it will overwrite the functional from pseudopotentials and performs calculation with whichever functional the user prefers. We further offer two ways of supplying exchange-correlation functional. The first is using 'short-hand' names such as 'LDA', 'PBE', 'SCAN'. A complete list of 'short-hand' expressions can be found in [the source code](../../../source/module_hamilt_general/module_xc/xc_functional.cpp). The other way is only available when ***compiling with LIBXC***, and it allows for supplying exchange-correlation functionals as combinations of LIBXC keywords for functional components, joined by a plus sign, for example, 'dft_functional='LDA_X_1D_EXPONENTIAL+LDA_C_1D_CSC'. The list of LIBXC keywords can be found on its [website](https://www.tddft.org/programs/libxc/functionals/). In this way, **we support all the LDA,GGA and mGGA functionals provided by LIBXC**.
+  On the other hand, if dft_functional is specified, it will overwrite the functional from pseudopotentials and performs calculation with whichever functional the user prefers. We further offer two ways of supplying exchange-correlation functional. The first is using 'short-hand' names such as 'LDA', 'PBE', 'SCAN'. A complete list of 'short-hand' expressions can be found in [the source code](../../../source/module_hamilt_general/module_xc/xc_functional.cpp). The other way is only available when ***compiling with LIBXC***, and it allows for supplying exchange-correlation functionals as combinations of LIBXC keywords for functional components, joined by a plus sign, for example, dft_functional='LDA_X_1D_EXPONENTIAL+LDA_C_1D_CSC'. The list of LIBXC keywords can be found on its [website](https://libxc.gitlab.io/functionals/). In this way, **we support all the LDA,GGA and mGGA functionals provided by LIBXC**.
 
   Furthermore, the old INPUT parameter exx_hybrid_type for hybrid functionals has been absorbed into dft_functional. Options are `hf` (pure Hartree-Fock), `pbe0`(PBE0), `hse` (Note: in order to use HSE functional, LIBXC is required). Note also that HSE has been tested while PBE0 has NOT been fully tested yet, and the maximum CPU cores for running exx in parallel is $N(N+1)/2$, with N being the number of atoms. And forces for hybrid functionals are not supported yet.
 
@@ -1760,14 +1760,14 @@ The band (KS orbital) energy for each (k-point, spin, band) will be printed in t
 
 - **Type**: Boolean
 - **Availability**: Numerical atomic orbital basis
-- **Description**: Whether to print Hamiltonian matrices H(R)/density matrics DM(R) in npz format. This feature does not work for gamma-only calculations. Currently only intended for internal usage.
+- **Description**: Whether to print Hamiltonian matrices $H(R)$/density matrics $DM(R)$ in npz format. This feature does not work for gamma-only calculations. Currently only intended for internal usage.
 - **Default**: False
 
 ### dm_to_rho
 
 - **Type**: Boolean
 - **Availability**: Numerical atomic orbital basis
-- **Description**: Reads density matrix DM(R) in npz format and creates electron density on grids. This feature does not work for gamma-only calculations. Only supports serial calculations. Currently only intended for internal usage.
+- **Description**: Reads density matrix $DM(R)$ in npz format and creates electron density on grids. This feature does not work for gamma-only calculations. Only supports serial calculations. Currently only intended for internal usage.
 - **Default**: False
 
 ### out_app_flag

source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.hpp

@@ -1,8 +1,4 @@
-//=======================
-// AUTHOR : Peize Lin
 #include "module_parameter/parameter.h"
-// DATE :   2022-09-13
-//=======================
 
 #ifndef LCAO_HAMILT_HPP
 #define LCAO_HAMILT_HPP
@@ -23,6 +19,7 @@
 
 #ifdef __EXX
 // Peize Lin add 2022.09.13
+
 template <typename Tdata>
 void sparse_format::cal_HR_exx(
     const Parallel_Orbitals& pv,
@@ -31,8 +28,8 @@ void sparse_format::cal_HR_exx(
     const double& sparse_threshold,
     const int (&nmp)[3],
     const std::vector<std::map<int,
-                               std::map<std::pair<int, std::array<int, 3>>,
-                                        RI::Tensor<Tdata>>>>& Hexxs) {
+                      std::map<std::pair<int, std::array<int, 3>>,
+                      RI::Tensor<Tdata>>>>& Hexxs) {
     ModuleBase::TITLE("sparse_format", "cal_HR_exx");
     ModuleBase::timer::tick("sparse_format", "cal_HR_exx");
 
@@ -45,7 +42,7 @@ void sparse_format::cal_HR_exx(
                 .tau[GlobalC::ucell.iat2ia[iat]]);
     }
     const std::array<std::array<double, 3>, 3> latvec
-        = {RI_Util::Vector3_to_array3(GlobalC::ucell.a1),
+        = {RI_Util::Vector3_to_array3(GlobalC::ucell.a1), // too bad to use GlobalC here, 
            RI_Util::Vector3_to_array3(GlobalC::ucell.a2),
            RI_Util::Vector3_to_array3(GlobalC::ucell.a3)};
 
@@ -58,18 +55,22 @@ void sparse_format::cal_HR_exx(
                                          ? std::vector<int>{current_spin}
                                          : std::vector<int>{0, 1, 2, 3};
 
-    for (const int is: is_list) {
+    for (const int is: is_list) 
+    {
         int is0_b = 0;
         int is1_b = 0;
         std::tie(is0_b, is1_b) = RI_2D_Comm::split_is_block(is);
 
-        if (Hexxs.empty()) {
+        if (Hexxs.empty()) 
+        {
             break;
         }
 
-        for (const auto& HexxA: Hexxs[is]) {
+        for (const auto& HexxA: Hexxs[is]) 
+        {
             const int iat0 = HexxA.first;
-            for (const auto& HexxB: HexxA.second) {
+            for (const auto& HexxB: HexxA.second) 
+            {
                 const int iat1 = HexxB.first.first;
 
                 const Abfs::Vector3_Order<int> R = RI_Util::array3_to_Vector3(
@@ -81,47 +82,60 @@ void sparse_format::cal_HR_exx(
 
                 const RI::Tensor<Tdata>& Hexx = HexxB.second;
 
-                for (size_t iw0 = 0; iw0 < Hexx.shape[0]; ++iw0) {
+                for (size_t iw0 = 0; iw0 < Hexx.shape[0]; ++iw0) 
+                {
                     const int iwt0 = RI_2D_Comm::get_iwt(iat0, iw0, is0_b);
                     const int iwt0_local = pv.global2local_row(iwt0);
 
-                    if (iwt0_local < 0) {
+                    if (iwt0_local < 0) 
+                    {
                         continue;
                     }
 
-                    for (size_t iw1 = 0; iw1 < Hexx.shape[1]; ++iw1) {
+                    for (size_t iw1 = 0; iw1 < Hexx.shape[1]; ++iw1) 
+                    {
                         const int iwt1 = RI_2D_Comm::get_iwt(iat1, iw1, is1_b);
                         const int iwt1_local = pv.global2local_col(iwt1);
 
-                        if (iwt1_local < 0) {
+                        if (iwt1_local < 0) 
+                        {
                             continue;
                         }
 
-                        if (std::abs(Hexx(iw0, iw1)) > sparse_threshold) {
-                            if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) {
+                        if (std::abs(Hexx(iw0, iw1)) > sparse_threshold) 
+                        {
+                            if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) 
+                            {
                                 auto& HR_sparse_ptr
                                     = HS_Arrays
                                           .HR_sparse[current_spin][R][iwt0];
                                 double& HR_sparse = HR_sparse_ptr[iwt1];
                                 HR_sparse += RI::Global_Func::convert<double>(
                                     frac * Hexx(iw0, iw1));
-                                if (std::abs(HR_sparse) <= sparse_threshold) {
+                                if (std::abs(HR_sparse) <= sparse_threshold) 
+                                {
                                     HR_sparse_ptr.erase(iwt1);
                                 }
-                            } else if (PARAM.inp.nspin == 4) {
+                            } 
+                            else if (PARAM.inp.nspin == 4) 
+                            {
                                 auto& HR_sparse_ptr
                                     = HS_Arrays.HR_soc_sparse[R][iwt0];
+                                
                                 std::complex<double>& HR_sparse
                                     = HR_sparse_ptr[iwt1];
+                                
                                 HR_sparse += RI::Global_Func::convert<
-                                    std::complex<double>>(frac
-                                                          * Hexx(iw0, iw1));
-                                if (std::abs(HR_sparse) <= sparse_threshold) {
+                                    std::complex<double>>(frac * Hexx(iw0, iw1));
+                                
+                                if (std::abs(HR_sparse) <= sparse_threshold) 
+                                {
                                     HR_sparse_ptr.erase(iwt1);
                                 }
-                            } else {
-                                throw std::invalid_argument(
-                                    std::string(__FILE__) + " line "
+                            } 
+                            else 
+                            {
+                                throw std::invalid_argument(std::string(__FILE__) + " line "
                                     + std::to_string(__LINE__));
                             }
                         }

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h

@@ -11,16 +11,19 @@ namespace hamilt {
 template <typename TK, typename TR>
 class OperatorLCAO : public Operator<TK> {
   public:
+
     OperatorLCAO(
         HS_Matrix_K<TK>* hsk_in, 
-        const std::vector<ModuleBase::Vector3<double>>& kvec_d_in,
-        HContainer<TR>* hR_in)
+        const std::vector<ModuleBase::Vector3<double>>& kvec_d_in, //! k-point vectors
+        HContainer<TR>* hR_in)  //! H(R) matrix, R is the Bravis lattice vector
         : hsk(hsk_in), kvec_d(kvec_d_in), hR(hR_in){}
+
     virtual ~OperatorLCAO()
     {
-        if (this->allocated_smatrix) {
+        if (this->allocated_smatrix) 
+        {
             delete[] this->smatrix_k;
-}
+        }
     }
 
     /* Function init(k) is used for update HR and HK ,
@@ -47,7 +50,8 @@ class OperatorLCAO : public Operator<TK> {
     Gamma_only case (TK = double), SK would not changed during one SCF loop, a template triangle matrix SK_temp is used
     for accelerating. General case (TK = std::complex<double>), only pointers of HK and SK saved in OperatorLCAO
     */
-    void matrixHk(MatrixBlock<TK>& hk_in, MatrixBlock<TK>& sk_in) {
+    void matrixHk(MatrixBlock<TK>& hk_in, MatrixBlock<TK>& sk_in) 
+    {
         this->get_hs_pointers();
 #ifdef __MPI
         hk_in = MatrixBlock<TK>{hmatrix_k,
@@ -59,8 +63,15 @@ class OperatorLCAO : public Operator<TK> {
                                (size_t)this->hsk->get_pv()->ncol,
                                this->hsk->get_pv()->desc};
 #else
-        hk_in = MatrixBlock<TK>{hmatrix_k, (size_t)this->hsk->get_pv()->nrow, (size_t)this->hsk->get_pv()->ncol, nullptr};
-        sk_in = MatrixBlock<TK>{smatrix_k, (size_t)this->hsk->get_pv()->nrow, (size_t)this->hsk->get_pv()->ncol, nullptr};
+        hk_in = MatrixBlock<TK>{hmatrix_k, 
+                                (size_t)this->hsk->get_pv()->nrow, 
+                                (size_t)this->hsk->get_pv()->ncol, 
+                                nullptr};
+        
+        sk_in = MatrixBlock<TK>{smatrix_k, 
+                                (size_t)this->hsk->get_pv()->nrow, 
+                                (size_t)this->hsk->get_pv()->ncol, 
+                                nullptr};
 #endif
     }
 
@@ -77,7 +88,7 @@ class OperatorLCAO : public Operator<TK> {
     virtual void set_HR_fixed(void*) { return; }
 
     /**
-     * @brief reset hr_done status
+     * @brief reset the status of 'hr_done' (if H(R) is calculated)
      */
     void set_hr_done(bool hr_done_in);
 
@@ -87,32 +98,33 @@ class OperatorLCAO : public Operator<TK> {
     void set_current_spin(const int current_spin_in);
 
     // protected:
-    //  Hamiltonian matrix which are calculated in OperatorLCAO
+    // Hamiltonian matrices which are calculated in OperatorLCAO
     HS_Matrix_K<TK>* hsk = nullptr;
     const std::vector<ModuleBase::Vector3<double>>& kvec_d;
 
   protected:
     bool new_e_iteration = true;
 
-    // Real space Hamiltonian pointer
+    //! Real-space Hamiltonian pointer
     hamilt::HContainer<TR>* hR = nullptr;
 
-    // current spin index
+    //! current spin index
     int current_spin = 0;
 
-    // if HR is calculated
+    //! if H(R) is calculated
     bool hr_done = false;
 
   private:
+
     void get_hs_pointers();
 
-    // there are H and S matrix for each k point in reciprocal space
-    // type double for gamma_only case, type complex<double> for multi-k-points
-    // case
+    //! there are H and S matrix for each k point in reciprocal space
+    //! 'double' type for gamma_only case, 
+    //! 'complex<double>' type for multi k-points case
     TK* hmatrix_k = nullptr;
     TK* smatrix_k = nullptr;
 
-    // only used for Gamma_only case
+    //! only used for Gamma_only case
     bool allocated_smatrix = false;
 };
 

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.cpp

@@ -41,7 +41,8 @@ void hamilt::OverlapNew<hamilt::OperatorLCAO<TK, TR>>::initialize_SR(Grid_Driver
     for (int iat1 = 0; iat1 < ucell->nat; iat1++)
     {
         auto tau1 = ucell->get_tau(iat1);
-        int T1, I1;
+        int T1=0;
+        int I1=0;
         ucell->iat2iait(iat1, &I1, &T1);
         AdjacentAtomInfo adjs;
         GridD->Find_atom(*ucell, tau1, T1, I1, &adjs);
@@ -84,8 +85,8 @@ void hamilt::OverlapNew<hamilt::OperatorLCAO<TK, TR>>::calculate_SR()
     for (int iap = 0; iap < this->SR->size_atom_pairs(); ++iap)
     {
         hamilt::AtomPair<TR>& tmp = this->SR->get_atom_pair(iap);
-        int iat1 = tmp.get_atom_i();
-        int iat2 = tmp.get_atom_j();
+        const int iat1 = tmp.get_atom_i();
+        const int iat2 = tmp.get_atom_j();
         const Parallel_Orbitals* paraV = tmp.get_paraV();
 
         for (int iR = 0; iR < tmp.get_R_size(); ++iR)
@@ -116,9 +117,11 @@ void hamilt::OverlapNew<hamilt::OperatorLCAO<TK, TR>>::cal_SR_IJR(const int& iat
     // ---------------------------------------------
     // get info of orbitals of atom1 and atom2 from ucell
     // ---------------------------------------------
-    int T1, I1;
+    int T1=0;
+    int I1=0;
     this->ucell->iat2iait(iat1, &I1, &T1);
-    int T2, I2;
+    int T2=0;
+    int I2=0;
     this->ucell->iat2iait(iat2, &I2, &T2);
     Atom& atom1 = this->ucell->atoms[T1];
     Atom& atom2 = this->ucell->atoms[T2];
@@ -188,14 +191,15 @@ void hamilt::OverlapNew<hamilt::OperatorLCAO<TK, TR>>::contributeHR()
 template <typename TK, typename TR>
 void hamilt::OverlapNew<hamilt::OperatorLCAO<TK, TR>>::contributeHk(int ik)
 {
-    // if k vector is not changed, then do nothing and return, only for gamma_only case
+    //! if k vector is not changed, then do nothing and return, only for gamma_only case
     if (this->kvec_d[ik] == this->kvec_d_old && std::is_same<TK, double>::value)
     {
         return;
     }
     ModuleBase::TITLE("OverlapNew", "contributeHk");
     ModuleBase::timer::tick("OverlapNew", "contributeHk");
-    // set SK to zero and then calculate SK for each k vector
+    
+    //! set SK to zero and then calculate SK for each k vector
     this->hsk->set_zero_sk();
     if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver))
     {
@@ -207,6 +211,7 @@ void hamilt::OverlapNew<hamilt::OperatorLCAO<TK, TR>>::contributeHk(int ik)
         const int ncol = this->SR->get_atom_pair(0).get_paraV()->get_col_size();
         hamilt::folding_HR(*this->SR, this->hsk->get_sk(), this->kvec_d[ik], ncol, 0);
     }
+    
     // update kvec_d_old
     this->kvec_d_old = this->kvec_d[ik];