Merge branch 'lishui' into develop

Author: mohanchen

source/module_hamilt_general/operator.cpp

@@ -16,8 +16,10 @@ Operator<T, Device>::~Operator()
     {
         delete this->hpsi;
     }
+
     Operator* last = this->next_op;
     Operator* last_sub = this->next_sub_op;
+
     while (last != nullptr || last_sub != nullptr)
     {
         if (last_sub != nullptr)
@@ -58,7 +60,6 @@ typename Operator<T, Device>::hpsi_info Operator<T, Device>::hPsi(hpsi_info& inp
     T* hpsi_pointer = std::get<2>(input);
     if (this->in_place)
     {
-        // ModuleBase::GlobalFunc::COPYARRAY(this->hpsi->get_pointer(), hpsi_pointer, this->hpsi->size());
         syncmem_op()(hpsi_pointer, this->hpsi->get_pointer(), this->hpsi->size());
         delete this->hpsi;
         this->hpsi = new psi::Psi<T, Device>(hpsi_pointer, 
@@ -78,31 +79,35 @@ typename Operator<T, Device>::hpsi_info Operator<T, Device>::hPsi(hpsi_info& inp
                                         psi_input->get_nbasis(),
                                         true);
 
-        switch (op->get_act_type())
-        {
-        case 2:
-            op->act(psi_wrapper, *this->hpsi, nbands);
-            break;
-        default:
-            op->act(nbands,
-                    psi_input->get_nbasis(),
-                    psi_input->get_npol(),
-                    tmpsi_in,
-                    this->hpsi->get_pointer(),
-                    psi_input->get_current_nbas(),
-                    is_first_node);
-            break;
-        }
+		switch (op->get_act_type())
+		{
+			case 2:
+				op->act(psi_wrapper, *this->hpsi, nbands);
+				break;
+			default:
+				op->act(nbands,
+						psi_input->get_nbasis(),
+						psi_input->get_npol(),
+						tmpsi_in,
+						this->hpsi->get_pointer(),
+						psi_input->get_current_nbas(),
+						is_first_node);
+				break;
+		}
     };
 
     ModuleBase::timer::tick("Operator", "hPsi");
+
     call_act(this, true); // first node
+
     Operator* node((Operator*)this->next_op);
+
     while (node != nullptr)
     {
         call_act(node, false); // other nodes
         node = (Operator*)(node->next_op);
     }
+
     ModuleBase::timer::tick("Operator", "hPsi");
 
     return hpsi_info(this->hpsi, psi::Range(1, 0, 0, nbands / psi_input->get_npol()), hpsi_pointer);
@@ -157,6 +162,7 @@ void Operator<T, Device>::add(Operator* next)
     }
 }
 
+
 template <typename T, typename Device>
 T* Operator<T, Device>::get_hpsi(const hpsi_info& info) const
 {

source/module_hamilt_general/operator.h

@@ -73,8 +73,10 @@ class Operator
 
     /// developer-friendly interfaces for act() function
     /// interface type 2: input and change the Psi-type HPsi
-    // virtual void act(const psi::Psi<T, Device>& psi_in, psi::Psi<T, Device>& psi_out) const {};
-    virtual void act(const psi::Psi<T, Device>& psi_in, psi::Psi<T, Device>& psi_out, const int nbands) const {};
+	// virtual void act(const psi::Psi<T, Device>& psi_in, psi::Psi<T, Device>& psi_out) const {};
+	virtual void act(const psi::Psi<T, Device>& psi_in, 
+			psi::Psi<T, Device>& psi_out, 
+			const int nbands) const {};
 
     /// interface type 3: return a Psi-type HPsi
     // virtual psi::Psi<T> act(const psi::Psi<T,Device>& psi_in) const { return psi_in; };
@@ -125,4 +127,4 @@ class Operator
 
 } // end namespace hamilt
 
-#endif
+#endif

source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp

@@ -1,5 +1,5 @@
-#ifndef W_ABACUS_DEVELOP_ABACUS_DEVELOP_SOURCE_MODULE_HAMILT_LCAO_HAMILT_LCAODFT_HAMILT_LCAO_H
-#define W_ABACUS_DEVELOP_ABACUS_DEVELOP_SOURCE_MODULE_HAMILT_LCAO_HAMILT_LCAODFT_HAMILT_LCAO_H
+#ifndef HAMILT_LCAO_H 
+#define HAMILT_LCAO_H 
 
 #include "module_basis/module_nao/two_center_bundle.h"
 #include "module_cell/klist.h"
@@ -25,24 +25,30 @@ namespace hamilt
 {
 
 // template first for type of k space H matrix elements
-// template second for type of temporary matrix, gamma_only fix-gamma-matrix + S-gamma, multi-k fix-Real + S-Real
+// template second for type of temporary matrix, 
+// gamma_only fix-gamma-matrix + S-gamma, 
+// multi-k fix-Real + S-Real
 template <typename TK, typename TR>
 class HamiltLCAO : public Hamilt<TK>
 {
   public:
+
+
+    using TAC = std::pair<int, std::array<int, 3>>;
+
+
     /**
      * @brief Constructor of Hamiltonian for LCAO base
      * HR and SR will be allocated with Operators
      */
-    using TAC = std::pair<int, std::array<int, 3>>;
     HamiltLCAO(Gint_Gamma* GG_in,
                Gint_k* GK_in,
                const UnitCell& ucell,
                const Grid_Driver& grid_d,
-               const Parallel_Orbitals* paraV,
-               elecstate::Potential* pot_in,
-               const K_Vectors& kv_in,
-               const TwoCenterBundle& two_center_bundle,
+			   const Parallel_Orbitals* paraV,
+			   elecstate::Potential* pot_in,
+			   const K_Vectors& kv,
+			   const TwoCenterBundle& two_center_bundle,
                const LCAO_Orbitals& orb,
                elecstate::DensityMatrix<TK, double>* DM_in
 #ifdef __DEEPKS
@@ -57,13 +63,14 @@ class HamiltLCAO : public Hamilt<TK>
                std::vector<std::map<int, std::map<TAC, RI::Tensor<std::complex<double>>>>>* Hexxc = nullptr
 #endif
     );
+
     /**
      * @brief Constructor of vacuum Operators, only HR and SR will be initialed as empty HContainer
      */
     HamiltLCAO(const UnitCell& ucell,
                const Grid_Driver& grid_d,
                const Parallel_Orbitals* paraV,
-               const K_Vectors& kv_in,
+               const K_Vectors& kv,
                const TwoCenterIntegrator& intor_overlap_orb,
                const std::vector<double>& orb_cutoff);
 
@@ -80,82 +87,93 @@ class HamiltLCAO : public Hamilt<TK>
 
     /// get pointer of Operator<TK> ops
     Operator<TK>*& getOperator();
-    /// get hk-pointer
+
+    /// get H(k) pointer
     TK* getHk() const
     {
         return this->hsk->get_hk();
     }
-    /// get sk-pointer
+
+    /// get S(k) pointer
     TK* getSk() const
     {
         return this->hsk->get_sk();
     }
+
     int get_size_hsk() const
     {
         return this->hsk->get_size();
     }
+
     /// get HR pointer of *this->hR, which is a HContainer<TR> and contains H(R)
     HContainer<TR>*& getHR()
     {
         return this->hR;
     }
+
     /// get SR pointer of *this->sR, which is a HContainer<TR> and contains S(R)
     HContainer<TR>*& getSR()
     {
         return this->sR;
     }
+
 #ifdef __DEEPKS
     /// get V_delta_R pointer of *this->V_delta_R, which is a HContainer<TR> and contains V_delta(R)
     HContainer<TR>*& get_V_delta_R()
     {
         return this->V_delta_R;
     }
 #endif
+
     /// refresh the status of HR
     void refresh() override;
 
     // for target K point, update consequence of hPsi() and matrix()
-    virtual void updateHk(const int ik) override;
+    void updateHk(const int ik, const int isk);
 
     /**
      * @brief special for LCAO, update SK only
      *
      * @param ik target K point
+     * @param kvec_d: direct coordinates of k-points
      * @param hk_type 0: SK is row-major, 1: SK is collumn-major
      * @return void
      */
-    void updateSk(const int ik, const int hk_type = 0);
+	void updateSk(const int ik, 
+			const ModuleBase::Vector3<double>& kvec_d,
+			const int hk_type = 0);
 
     // core function: return H(k) and S(k) matrixs for direct solving eigenvalues.
     // not used in PW base
-    // void matrix(MatrixBlock<std::complex<double>> &hk_in, MatrixBlock<std::complex<double>> &sk_in) override;
     void matrix(MatrixBlock<TK>& hk_in, MatrixBlock<TK>& sk_in) override;
 
   private:
-    const K_Vectors* kv = nullptr;
 
-    // Real space Hamiltonian
+    //! Real space Hamiltonian H(R), where R is the Bravis lattice vector
     HContainer<TR>* hR = nullptr;
+
+    //! Real space overlap matrix S(R), where R is the Bravis lattice vector
     HContainer<TR>* sR = nullptr;
 
 #ifdef __DEEPKS
     HContainer<TR>* V_delta_R = nullptr;
 #endif
 
+    //! Hamiltonian and overlap matrices for a specific k point
     HS_Matrix_K<TK>* hsk = nullptr;
+
     // special case for NSPIN=2 , data of HR should be separated into two parts
     // save them in this->hRS2;
     std::vector<TR> hRS2;
+
     int refresh_times = 1;
 
-    /// current_spin for NSPIN=2, 0: hamiltonian for spin up, 1: hamiltonian for spin down
+    //! current_spin for NSPIN=2 case 
+    //! 0: Hamiltonian for spin up, 
+    //! 1: Hamiltonian for spin down
     int current_spin = 0;
 
     const int istep = 0;
-
-    // sk and hk will be refactored to HamiltLCAO later
-    // std::vector<TK> sk;
-    // std::vector<TK> hk;
 };
 
 } // namespace hamilt

source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h

@@ -1,5 +1,5 @@
-#ifndef ABACUS_HS_MATRIX_K_HPP
-#define ABACUS_HS_MATRIX_K_HPP
+#ifndef HS_MATRIX_K_HPP
+#define HS_MATRIX_K_HPP
 
 #include "module_basis/module_ao/parallel_orbitals.h"
 
@@ -9,6 +9,7 @@ namespace hamilt
     template <typename TK>
     class HS_Matrix_K
     {
+
         public:
             HS_Matrix_K(const Parallel_Orbitals* paraV, bool no_s=false){
                 hk.resize(paraV->nloc);
@@ -18,17 +19,27 @@ namespace hamilt
                 }
                 this->pv = paraV;
             }
+
             TK* get_hk() {return hk.data();}
+
             TK* get_sk() {return sk.data();}
+
             int get_size() {return hk.size();}
+
             void set_zero_hk() {hk.assign(hk.size(), 0);}
+
             void set_zero_sk() {sk.assign(sk.size(), 0);}
+
             const Parallel_Orbitals* get_pv() const {return this->pv;}
+
         private:
+
             std::vector<TK> hk;
+
             std::vector<TK> sk;
+
             const Parallel_Orbitals* pv = nullptr;
     };
 }
 
-#endif
+#endif

source/module_hamilt_lcao/hamilt_lcaodft/hs_matrix_k.hpp

@@ -100,6 +100,8 @@ class OperatorLCAO : public Operator<TK> {
     // protected:
     // Hamiltonian matrices which are calculated in OperatorLCAO
     HS_Matrix_K<TK>* hsk = nullptr;
+
+    // kvec_d
     const std::vector<ModuleBase::Vector3<double>>& kvec_d;
 
   protected:

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h

@@ -31,6 +31,7 @@ TDEkinetic<OperatorLCAO<TK, TR>>::TDEkinetic(HS_Matrix_K<TK>* hsk_in,
     // initialize HR to get adjs info.
     this->initialize_HR(Grid);
 }
+
 template <typename TK, typename TR>
 TDEkinetic<OperatorLCAO<TK, TR>>::~TDEkinetic()
 {
@@ -40,12 +41,14 @@ TDEkinetic<OperatorLCAO<TK, TR>>::~TDEkinetic()
     }
     TD_Velocity::td_vel_op = nullptr;
 }
+
 // term A^2*S
 template <typename TK, typename TR>
 void TDEkinetic<OperatorLCAO<TK, TR>>::td_ekinetic_scalar(std::complex<double>* Hloc,const TR& overlap, int nnr)
 {
     return;
 }
+
 // term A^2*S
 template <>
 void TDEkinetic<OperatorLCAO<std::complex<double>, double>>::td_ekinetic_scalar(std::complex<double>* Hloc,
@@ -57,6 +60,7 @@ void TDEkinetic<OperatorLCAO<std::complex<double>, double>>::td_ekinetic_scalar(
     Hloc[nnr] += tmp;
     return;
 }
+
 // term A dot ∇
 template <typename TK, typename TR>
 void TDEkinetic<OperatorLCAO<TK, TR>>::td_ekinetic_grad(std::complex<double>* Hloc,
@@ -384,8 +388,11 @@ void TDEkinetic<OperatorLCAO<std::complex<double>, double>>::contributeHk(int ik
         const Parallel_Orbitals* paraV = this->hR_tmp->get_atom_pair(0).get_paraV();
         // save HR data for output
         int spin_tot = PARAM.inp.nspin;
+
         if (spin_tot == 4)
-            ;
+        {
+
+        }
         else if (!output_hR_done && TD_Velocity::out_mat_R)
         {
             for (int spin_now = 0; spin_now < spin_tot; spin_now++)
@@ -398,6 +405,7 @@ void TDEkinetic<OperatorLCAO<std::complex<double>, double>>::contributeHk(int ik
             }
             output_hR_done = true;
         }
+
         // folding inside HR to HK
         if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver))
         {