Keep updating esolver and related source files

source/driver_run.cpp

@@ -19,16 +19,16 @@
  * Esolver::Run takes in a configuration and provides force and stress, 
  * the configuration-changing subroutine takes force and stress and updates the configuration
  */
-void Driver::driver_run()
+void Driver::driver_run(void)
 {
     ModuleBase::TITLE("Driver", "driver_line");
     ModuleBase::timer::tick("Driver", "driver_line");
 
-    // 1. Determine type of Esolver
+    //! 1: initialize the ESolver 
     ModuleESolver::ESolver *p_esolver = nullptr;
     ModuleESolver::init_esolver(p_esolver);
 
-    // 2. Setup cell and atom information
+    //! 2: setup cell and atom information
 #ifndef __LCAO
     if(GlobalV::BASIS_TYPE == "lcao_in_pw" || GlobalV::BASIS_TYPE == "lcao")
     {
@@ -37,30 +37,29 @@ void Driver::driver_run()
 #endif
     GlobalC::ucell.setup_cell(GlobalV::stru_file, GlobalV::ofs_running);
 
-    // 3. For these two types of calculations
+    //! 3: for these two types of calculations
     // nothing else need to be initialized
-    if(GlobalV::CALCULATION == "test_neighbour" || GlobalV::CALCULATION == "test_memory")
+    if(GlobalV::CALCULATION == "test_neighbour" 
+    || GlobalV::CALCULATION == "test_memory")
     {
-        p_esolver->Run(0, GlobalC::ucell);
+        p_esolver->run(0, GlobalC::ucell);
         ModuleBase::QUIT();
     }
 
-    // 4. Initialize Esolver,and fill json-structure 
-    p_esolver->Init(INPUT, GlobalC::ucell);
+    //! 4: initialize Esolver and fill json-structure 
+    p_esolver->init(INPUT, GlobalC::ucell);
 
 
 #ifdef __RAPIDJSON
     Json::gen_stru_wrapper(&GlobalC::ucell);
 #endif
 
-    //------------------------------------------------------------
-    // This part onward needs to be refactored.
-    //---------------------------MD/Relax-------------------------
+    //! 5: md or relax calculations 
     if(GlobalV::CALCULATION == "md")
     {
         Run_MD::md_line(GlobalC::ucell, p_esolver, INPUT.mdp);
     }
-    else // scf; cell relaxation; nscf; etc
+    else //! scf; cell relaxation; nscf; etc
     {
         if (GlobalV::precision_flag == "single")
         {
@@ -73,10 +72,9 @@ void Driver::driver_run()
             rl_driver.relax_driver(p_esolver);
         }
     }
-    //---------------------------MD/Relax------------------
 
-    // 6. clean up esolver
-    p_esolver->postprocess();
+    //! 6: clean up esolver
+    p_esolver->post_process();
     ModuleESolver::clean_esolver(p_esolver);
 
     ModuleBase::timer::tick("Driver", "driver_line");

source/module_elecstate/elecstate_lcao.cpp

@@ -14,8 +14,10 @@ namespace elecstate
 template <>
 void ElecStateLCAO<double>::print_psi(const psi::Psi<double>& psi_in, const int istep)
 {
-    if (!ElecStateLCAO<double>::out_wfc_lcao)
-        return;
+	if (!ElecStateLCAO<double>::out_wfc_lcao)
+	{
+		return;
+	}
 
     // output but not do  "2d-to-grid" conversion
     double** wfc_grid = nullptr;
@@ -28,8 +30,11 @@ void ElecStateLCAO<double>::print_psi(const psi::Psi<double>& psi_in, const int
 template <>
 void ElecStateLCAO<std::complex<double>>::print_psi(const psi::Psi<std::complex<double>>& psi_in, const int istep)
 {
-    if (!ElecStateLCAO<std::complex<double>>::out_wfc_lcao && !ElecStateLCAO<std::complex<double>>::need_psi_grid)
-        return;
+	if (!ElecStateLCAO<std::complex<double>>::out_wfc_lcao 
+			&& !ElecStateLCAO<std::complex<double>>::need_psi_grid)
+	{
+		return;
+	}
 
     // output but not do "2d-to-grid" conversion
     std::complex<double>** wfc_grid = nullptr;
@@ -38,6 +43,7 @@ void ElecStateLCAO<std::complex<double>>::print_psi(const psi::Psi<std::complex<
     {
         wfc_grid = this->lowf->wfc_k_grid[ik];
     }
+
 #ifdef __MPI
     this->lowf->wfc_2d_to_grid(istep,
                                ElecStateLCAO<std::complex<double>>::out_wfc_flag,
@@ -252,4 +258,4 @@ double ElecStateLCAO<std::complex<double>>::get_spin_constrain_energy()
 template class ElecStateLCAO<double>; // Gamma_only case
 template class ElecStateLCAO<std::complex<double>>; // multi-k case
 
-} // namespace elecstate
+} // namespace elecstate

source/module_elecstate/module_dm/cal_dm_psi.cpp

@@ -44,8 +44,10 @@ void cal_dm_psi(const Parallel_Orbitals* ParaV,
                     ModuleBase::WARNING_QUIT("ElecStateLCAO::cal_dm", "please check global2local_col!");
                 }
             }
-            if (ib_global >= wg.nc)
-                continue;
+			if (ib_global >= wg.nc)
+			{
+				continue;
+			}
             const double wg_local = wg(ik, ib_global);
             double* wg_wfc_pointer = &(wg_wfc(0, ib_local, 0));
             BlasConnector::scal(nbasis_local, wg_local, wg_wfc_pointer, 1);
@@ -107,9 +109,11 @@ void cal_dm_psi(const Parallel_Orbitals* ParaV,
                     ModuleBase::WARNING_QUIT("ElecStateLCAO::cal_dm", "please check global2local_col!");
                 }
             }
-            if (ib_global >= wg.nc)
-                continue;
-            const double wg_local = wg(ik, ib_global);
+			if (ib_global >= wg.nc)
+			{
+				continue;
+			}
+			const double wg_local = wg(ik, ib_global);
             std::complex<double>* wg_wfc_pointer = &(wg_wfc(0, ib_local, 0));
             BlasConnector::scal(nbasis_local, wg_local, wg_wfc_pointer, 1);
         }
@@ -233,7 +237,8 @@ void psiMulPsi(const psi::Psi<std::complex<double>>& psi1,
     const char N_char = 'N', T_char = 'T';
     const int nlocal = psi1.get_nbasis();
     const int nbands = psi1.get_nbands();
-    const std::complex<double> one_complex = {1.0, 0.0}, zero_complex = {0.0, 0.0};
+    const std::complex<double> one_complex = {1.0, 0.0};
+    const std::complex<double> zero_complex = {0.0, 0.0};
     zgemm_(&N_char,
            &T_char,
            &nlocal,

source/module_elecstate/module_dm/density_matrix.cpp

@@ -440,11 +440,11 @@ void DensityMatrix<TK,TR>::cal_DMR_test()
                 const int* r_index = tmp_ap.get_R_index(ir);
                 hamilt::BaseMatrix<TR>* tmp_matrix = tmp_ap.find_matrix(r_index[0], r_index[1], r_index[2]);
 #ifdef __DEBUG
-            if (tmp_matrix == nullptr)
-            {
-                std::cout << "tmp_matrix is nullptr" << std::endl;
-                continue;
-            }
+				if (tmp_matrix == nullptr)
+				{
+					std::cout << "tmp_matrix is nullptr" << std::endl;
+					continue;
+				}
 #endif
                 std::complex<TR> tmp_res;
                 // loop over k-points
@@ -515,44 +515,47 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR()
 #endif
                 // loop over k-points
                 if(GlobalV::NSPIN !=4 )
-                for (int ik = 0; ik < this->_nks; ++ik)
-                {
-                    // cal k_phase
-                    // if TK==std::complex<double>, kphase is e^{ikR}
-                    const ModuleBase::Vector3<double> dR(r_index[0], r_index[1], r_index[2]);
-                    const double arg = (this->_kv->kvec_d[ik] * dR) * ModuleBase::TWO_PI;
-                    double sinp, cosp;
-                    ModuleBase::libm::sincos(arg, &sinp, &cosp);
-                    std::complex<double> kphase = std::complex<double>(cosp, sinp);
-                    // set DMR element
-                    double* tmp_DMR_pointer = tmp_matrix->get_pointer();
-                    std::complex<double>* tmp_DMK_pointer = this->_DMK[ik + ik_begin].data();
-                    double* DMK_real_pointer = nullptr;
-                    double* DMK_imag_pointer = nullptr;
-                    // jump DMK to fill DMR
-                    // DMR is row-major, DMK is column-major
-                    tmp_DMK_pointer += col_ap * this->_paraV->nrow + row_ap;
-                    for (int mu = 0; mu < this->_paraV->get_row_size(iat1); ++mu)
-                    {
-                        DMK_real_pointer = (double*)tmp_DMK_pointer;
-                        DMK_imag_pointer = DMK_real_pointer + 1;
-                        BlasConnector::axpy(this->_paraV->get_col_size(iat2),
-                                            kphase.real(),
-                                            DMK_real_pointer,
-                                            ld_hk2,
-                                            tmp_DMR_pointer,
-                                            1);
-                        // "-" since i^2 = -1
-                        BlasConnector::axpy(this->_paraV->get_col_size(iat2),
-                                            -kphase.imag(),
-                                            DMK_imag_pointer,
-                                            ld_hk2,
-                                            tmp_DMR_pointer,
-                                            1);
-                        tmp_DMK_pointer += 1;
-                        tmp_DMR_pointer += this->_paraV->get_col_size(iat2);
-                    }
-                }
+				{
+					for (int ik = 0; ik < this->_nks; ++ik)
+					{
+						// cal k_phase
+						// if TK==std::complex<double>, kphase is e^{ikR}
+						const ModuleBase::Vector3<double> dR(r_index[0], r_index[1], r_index[2]);
+						const double arg = (this->_kv->kvec_d[ik] * dR) * ModuleBase::TWO_PI;
+						double sinp, cosp;
+						ModuleBase::libm::sincos(arg, &sinp, &cosp);
+						std::complex<double> kphase = std::complex<double>(cosp, sinp);
+						// set DMR element
+						double* tmp_DMR_pointer = tmp_matrix->get_pointer();
+						std::complex<double>* tmp_DMK_pointer = this->_DMK[ik + ik_begin].data();
+						double* DMK_real_pointer = nullptr;
+						double* DMK_imag_pointer = nullptr;
+						// jump DMK to fill DMR
+						// DMR is row-major, DMK is column-major
+						tmp_DMK_pointer += col_ap * this->_paraV->nrow + row_ap;
+						for (int mu = 0; mu < this->_paraV->get_row_size(iat1); ++mu)
+						{
+							DMK_real_pointer = (double*)tmp_DMK_pointer;
+							DMK_imag_pointer = DMK_real_pointer + 1;
+							BlasConnector::axpy(this->_paraV->get_col_size(iat2),
+									kphase.real(),
+									DMK_real_pointer,
+									ld_hk2,
+									tmp_DMR_pointer,
+									1);
+							// "-" since i^2 = -1
+							BlasConnector::axpy(this->_paraV->get_col_size(iat2),
+									-kphase.imag(),
+									DMK_imag_pointer,
+									ld_hk2,
+									tmp_DMR_pointer,
+									1);
+							tmp_DMK_pointer += 1;
+							tmp_DMR_pointer += this->_paraV->get_col_size(iat2);
+						}
+					}
+				}
+
                 // treat DMR as pauli matrix when NSPIN=4
                 if(GlobalV::NSPIN==4)
                 {

source/module_elecstate/module_dm/density_matrix.h

@@ -15,28 +15,33 @@ namespace elecstate
  * <TK,TR> = <double,double> for Gamma-only calculation
  * <TK,TR> = <std::complex<double>,double> for multi-k calculation
  */
-    template<typename T> struct ShiftRealComplex
-    {
-        using type = void;
-    };
-    template<>
-    struct ShiftRealComplex<double> {
-        using type = std::complex<double>;
-    };
-    template<>
-    struct ShiftRealComplex<std::complex<double>> {
-        using type = double;
-    };
-
-    template <typename TK, typename TR>
-    class DensityMatrix
-    {
-        using TRShift = typename ShiftRealComplex<TR>::type;
-  public:
-    /**
-     * @brief Destructor of class DensityMatrix
-     */
-    ~DensityMatrix();
+template<typename T> struct ShiftRealComplex
+{
+    using type = void;
+};
+
+template<>
+struct ShiftRealComplex<double> 
+{
+	using type = std::complex<double>;
+};
+
+template<>
+struct ShiftRealComplex<std::complex<double>> 
+{
+	using type = double;
+};
+
+template <typename TK, typename TR>
+class DensityMatrix
+{
+	using TRShift = typename ShiftRealComplex<TR>::type;
+
+	public:
+	/**
+	 * @brief Destructor of class DensityMatrix
+	 */
+	~DensityMatrix();
 
     /**
      * @brief Constructor of class DensityMatrix for multi-k calculation
@@ -160,6 +165,7 @@ namespace elecstate
      * @brief get pointer of paraV
      */
     const Parallel_Orbitals* get_paraV_pointer() const;
+
     const K_Vectors* get_kv_pointer() const;
 
     /**