diff --git a/source/driver_run.cpp b/source/driver_run.cpp index 5222be6aa5..de6cb18193 100644 --- a/source/driver_run.cpp +++ b/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"); diff --git a/source/module_elecstate/elecstate_lcao.cpp b/source/module_elecstate/elecstate_lcao.cpp index df9b377376..6dd00a7f4f 100644 --- a/source/module_elecstate/elecstate_lcao.cpp +++ b/source/module_elecstate/elecstate_lcao.cpp @@ -14,8 +14,10 @@ namespace elecstate template <> void ElecStateLCAO::print_psi(const psi::Psi& psi_in, const int istep) { - if (!ElecStateLCAO::out_wfc_lcao) - return; + if (!ElecStateLCAO::out_wfc_lcao) + { + return; + } // output but not do "2d-to-grid" conversion double** wfc_grid = nullptr; @@ -28,8 +30,11 @@ void ElecStateLCAO::print_psi(const psi::Psi& psi_in, const int template <> void ElecStateLCAO>::print_psi(const psi::Psi>& psi_in, const int istep) { - if (!ElecStateLCAO>::out_wfc_lcao && !ElecStateLCAO>::need_psi_grid) - return; + if (!ElecStateLCAO>::out_wfc_lcao + && !ElecStateLCAO>::need_psi_grid) + { + return; + } // output but not do "2d-to-grid" conversion std::complex** wfc_grid = nullptr; @@ -38,6 +43,7 @@ void ElecStateLCAO>::print_psi(const psi::Psilowf->wfc_k_grid[ik]; } + #ifdef __MPI this->lowf->wfc_2d_to_grid(istep, ElecStateLCAO>::out_wfc_flag, @@ -252,4 +258,4 @@ double ElecStateLCAO>::get_spin_constrain_energy() template class ElecStateLCAO; // Gamma_only case template class ElecStateLCAO>; // multi-k case -} // namespace elecstate \ No newline at end of file +} // namespace elecstate diff --git a/source/module_elecstate/module_dm/cal_dm_psi.cpp b/source/module_elecstate/module_dm/cal_dm_psi.cpp index 4164294b75..58d58ea959 100644 --- a/source/module_elecstate/module_dm/cal_dm_psi.cpp +++ b/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* 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>& psi1, const char N_char = 'N', T_char = 'T'; const int nlocal = psi1.get_nbasis(); const int nbands = psi1.get_nbands(); - const std::complex one_complex = {1.0, 0.0}, zero_complex = {0.0, 0.0}; + const std::complex one_complex = {1.0, 0.0}; + const std::complex zero_complex = {0.0, 0.0}; zgemm_(&N_char, &T_char, &nlocal, diff --git a/source/module_elecstate/module_dm/density_matrix.cpp b/source/module_elecstate/module_dm/density_matrix.cpp index 865245f905..bdd1947d62 100644 --- a/source/module_elecstate/module_dm/density_matrix.cpp +++ b/source/module_elecstate/module_dm/density_matrix.cpp @@ -440,11 +440,11 @@ void DensityMatrix::cal_DMR_test() const int* r_index = tmp_ap.get_R_index(ir); hamilt::BaseMatrix* 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 tmp_res; // loop over k-points @@ -515,44 +515,47 @@ void DensityMatrix, 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, kphase is e^{ikR} - const ModuleBase::Vector3 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 kphase = std::complex(cosp, sinp); - // set DMR element - double* tmp_DMR_pointer = tmp_matrix->get_pointer(); - std::complex* 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, kphase is e^{ikR} + const ModuleBase::Vector3 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 kphase = std::complex(cosp, sinp); + // set DMR element + double* tmp_DMR_pointer = tmp_matrix->get_pointer(); + std::complex* 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) { diff --git a/source/module_elecstate/module_dm/density_matrix.h b/source/module_elecstate/module_dm/density_matrix.h index 519f798f5d..1c758ae139 100644 --- a/source/module_elecstate/module_dm/density_matrix.h +++ b/source/module_elecstate/module_dm/density_matrix.h @@ -15,28 +15,33 @@ namespace elecstate * = for Gamma-only calculation * = ,double> for multi-k calculation */ - template struct ShiftRealComplex - { - using type = void; - }; - template<> - struct ShiftRealComplex { - using type = std::complex; - }; - template<> - struct ShiftRealComplex> { - using type = double; - }; - - template - class DensityMatrix - { - using TRShift = typename ShiftRealComplex::type; - public: - /** - * @brief Destructor of class DensityMatrix - */ - ~DensityMatrix(); +template struct ShiftRealComplex +{ + using type = void; +}; + +template<> +struct ShiftRealComplex +{ + using type = std::complex; +}; + +template<> +struct ShiftRealComplex> +{ + using type = double; +}; + +template +class DensityMatrix +{ + using TRShift = typename ShiftRealComplex::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; /** diff --git a/source/module_esolver/esolver.cpp b/source/module_esolver/esolver.cpp index 336ed32b1b..fb366dfaf3 100644 --- a/source/module_esolver/esolver.cpp +++ b/source/module_esolver/esolver.cpp @@ -13,152 +13,173 @@ namespace ModuleESolver { - void ESolver::printname() - { - std::cout << classname << std::endl; - } - std::string determine_type() - { - std::string esolver_type = "none"; - if (GlobalV::BASIS_TYPE == "pw") - { - if(GlobalV::ESOLVER_TYPE == "sdft") - { - esolver_type = "sdft_pw"; - } - else if(GlobalV::ESOLVER_TYPE == "ofdft") - { - esolver_type = "ofdft"; - } - else if(GlobalV::ESOLVER_TYPE == "ksdft") - { - esolver_type = "ksdft_pw"; - } - } - else if (GlobalV::BASIS_TYPE == "lcao_in_pw") - { +void ESolver::printname(void) +{ + std::cout << classname << std::endl; +} + +std::string determine_type(void) +{ + std::string esolver_type = "none"; + if (GlobalV::BASIS_TYPE == "pw") + { + if(GlobalV::ESOLVER_TYPE == "sdft") + { + esolver_type = "sdft_pw"; + } + else if(GlobalV::ESOLVER_TYPE == "ofdft") + { + esolver_type = "ofdft"; + } + else if(GlobalV::ESOLVER_TYPE == "ksdft") + { + esolver_type = "ksdft_pw"; + } + } + else if (GlobalV::BASIS_TYPE == "lcao_in_pw") + { #ifdef __LCAO - if(GlobalV::ESOLVER_TYPE == "sdft") - { - esolver_type = "sdft_pw"; - } - else if(GlobalV::ESOLVER_TYPE == "ksdft") - { - esolver_type = "ksdft_pw"; - } + if(GlobalV::ESOLVER_TYPE == "sdft") + { + esolver_type = "sdft_pw"; + } + else if(GlobalV::ESOLVER_TYPE == "ksdft") + { + esolver_type = "ksdft_pw"; + } #else - ModuleBase::WARNING_QUIT("ESolver", "LCAO basis type must be compiled with __LCAO"); + ModuleBase::WARNING_QUIT("ESolver", "LCAO basis type must be compiled with __LCAO"); #endif - } - else if (GlobalV::BASIS_TYPE == "lcao") - { + } + else if (GlobalV::BASIS_TYPE == "lcao") + { #ifdef __LCAO - if(GlobalV::ESOLVER_TYPE == "tddft") - { - esolver_type = "ksdft_lcao_tddft"; - } - else if(GlobalV::ESOLVER_TYPE == "ksdft") - { - esolver_type = "ksdft_lcao"; - } + if(GlobalV::ESOLVER_TYPE == "tddft") + { + esolver_type = "ksdft_lcao_tddft"; + } + else if(GlobalV::ESOLVER_TYPE == "ksdft") + { + esolver_type = "ksdft_lcao"; + } #else - ModuleBase::WARNING_QUIT("ESolver", "LCAO basis type must be compiled with __LCAO"); + ModuleBase::WARNING_QUIT("ESolver", "LCAO basis type must be compiled with __LCAO"); #endif - } + } + + if(GlobalV::ESOLVER_TYPE == "lj") + { + esolver_type = "lj_pot"; + } + else if(GlobalV::ESOLVER_TYPE == "dp") + { + esolver_type = "dp_pot"; + } + else if(esolver_type == "none") + { + ModuleBase::WARNING_QUIT("ESolver", "No such esolver_type combined with basis_type"); + } + + GlobalV::ofs_running << " The esolver type has been set to : " << esolver_type << std::endl; + + auto device_info = GlobalV::device_flag; + + for (char &c : device_info) + { + if (std::islower(c)) + { + c = std::toupper(c); + } + } + if (GlobalV::MY_RANK == 0) + { + std::cout << " RUNNING WITH DEVICE : " << device_info << " / " + << psi::device::get_device_info(GlobalV::device_flag) << std::endl; + } - if(GlobalV::ESOLVER_TYPE == "lj") - { - esolver_type = "lj_pot"; - } - else if(GlobalV::ESOLVER_TYPE == "dp") - { - esolver_type = "dp_pot"; - } - else if(esolver_type == "none") - { - ModuleBase::WARNING_QUIT("ESolver", "No such esolver_type combined with basis_type"); - } + GlobalV::ofs_running << "\n RUNNING WITH DEVICE : " << device_info << " / " + << psi::device::get_device_info(GlobalV::device_flag) << std::endl; + + return esolver_type; +} - GlobalV::ofs_running << " The esolver type has been set to : " << esolver_type << std::endl; - auto device_info = GlobalV::device_flag; - for (char &c : device_info) { - if (std::islower(c)) { - c = std::toupper(c); - } - } - if (GlobalV::MY_RANK == 0) { - std::cout << " RUNNING WITH DEVICE : " << device_info << " / " - << psi::device::get_device_info(GlobalV::device_flag) << std::endl; - } - GlobalV::ofs_running << "\n RUNNING WITH DEVICE : " << device_info << " / " - << psi::device::get_device_info(GlobalV::device_flag) << std::endl; - return esolver_type; - } - //Some API to operate E_Solver - void init_esolver(ESolver*& p_esolver) - { - //determine type of esolver based on INPUT information - std::string esolver_type = determine_type(); +//Some API to operate E_Solver +void init_esolver(ESolver*& p_esolver) +{ + //determine type of esolver based on INPUT information + std::string esolver_type = determine_type(); - //initialize the corresponding Esolver child class - if (esolver_type == "ksdft_pw") - { - #if ((defined __CUDA) || (defined __ROCM)) - if (GlobalV::device_flag == "gpu") { - if (GlobalV::precision_flag == "single") { - p_esolver = new ESolver_KS_PW, psi::DEVICE_GPU>(); - } - else { - p_esolver = new ESolver_KS_PW, psi::DEVICE_GPU>(); - } - return; - } - #endif - if (GlobalV::precision_flag == "single") { - p_esolver = new ESolver_KS_PW, psi::DEVICE_CPU>(); - } - else { - p_esolver = new ESolver_KS_PW, psi::DEVICE_CPU>(); - } - } + //initialize the corresponding Esolver child class + if (esolver_type == "ksdft_pw") + { +#if ((defined __CUDA) || (defined __ROCM)) + if (GlobalV::device_flag == "gpu") + { + if (GlobalV::precision_flag == "single") + { + p_esolver = new ESolver_KS_PW, psi::DEVICE_GPU>(); + } + else + { + p_esolver = new ESolver_KS_PW, psi::DEVICE_GPU>(); + } + return; + } +#endif + if (GlobalV::precision_flag == "single") + { + p_esolver = new ESolver_KS_PW, psi::DEVICE_CPU>(); + } + else + { + p_esolver = new ESolver_KS_PW, psi::DEVICE_CPU>(); + } + } #ifdef __LCAO - else if (esolver_type == "ksdft_lcao") - { - if (GlobalV::GAMMA_ONLY_LOCAL) - p_esolver = new ESolver_KS_LCAO(); - else if (GlobalV::NSPIN < 4) - p_esolver = new ESolver_KS_LCAO, double>(); - else - p_esolver = new ESolver_KS_LCAO, std::complex>(); - } - else if (esolver_type == "ksdft_lcao_tddft") - { - p_esolver = new ESolver_KS_LCAO_TDDFT(); - } + else if (esolver_type == "ksdft_lcao") + { + if (GlobalV::GAMMA_ONLY_LOCAL) + { + p_esolver = new ESolver_KS_LCAO(); + } + else if (GlobalV::NSPIN < 4) + { + p_esolver = new ESolver_KS_LCAO, double>(); + } + else + { + p_esolver = new ESolver_KS_LCAO, std::complex>(); + } + } + else if (esolver_type == "ksdft_lcao_tddft") + { + p_esolver = new ESolver_KS_LCAO_TDDFT(); + } #endif - else if (esolver_type == "sdft_pw") - { - p_esolver = new ESolver_SDFT_PW(); - } - else if(esolver_type == "ofdft") - { - p_esolver = new ESolver_OF(); - } - else if (esolver_type == "lj_pot") - { - p_esolver = new ESolver_LJ(); - } - else if (esolver_type == "dp_pot") - { - p_esolver = new ESolver_DP(INPUT.mdp.pot_file); - } - } + else if (esolver_type == "sdft_pw") + { + p_esolver = new ESolver_SDFT_PW(); + } + else if(esolver_type == "ofdft") + { + p_esolver = new ESolver_OF(); + } + else if (esolver_type == "lj_pot") + { + p_esolver = new ESolver_LJ(); + } + else if (esolver_type == "dp_pot") + { + p_esolver = new ESolver_DP(INPUT.mdp.pot_file); + } +} + - void clean_esolver(ESolver*& pesolver) - { - delete pesolver; - } +void clean_esolver(ESolver*& pesolver) +{ + delete pesolver; +} } diff --git a/source/module_esolver/esolver.h b/source/module_esolver/esolver.h index 88345617ce..142efa2a4a 100644 --- a/source/module_esolver/esolver.h +++ b/source/module_esolver/esolver.h @@ -9,8 +9,6 @@ namespace ModuleESolver { class ESolver { - // protected: - // ModuleBase::matrix lattice_v; public: ESolver() { @@ -21,24 +19,27 @@ class ESolver { } - // virtual void Init(Input_EnSolver &inp, matrix &lattice_v)=0 - virtual void Init(Input& inp, UnitCell& cell) = 0; + //! initialize the energy solver by using input parameters and cell modules + virtual void init(Input& inp, UnitCell& cell) = 0; - // They shoud be add after atom class is refactored - // virtual void UpdateLatAtom(ModuleBase::matrix &lat_in, Atom &atom_in); - // virtual void UpdateLat(ModuleBase::matrix &lat_in); - // virtual void UpdateAtom(Atom &atom_in); + //! run energy solver + virtual void run(int istep, UnitCell& cell) = 0; - virtual void Run(int istep, UnitCell& cell) = 0; + //! deal with exx and other calculation than scf/md/relax: + //! such as nscf, get_wf and get_pchg + virtual void others(const int istep){}; - // Deal with exx and other calculation than scf/md/relax: - // such as nscf, get_wf and get_pchg - virtual void othercalculation(const int istep){}; + //! calculate total energy of a given system + virtual double cal_energy() = 0; - virtual double cal_Energy() = 0; - virtual void cal_Force(ModuleBase::matrix& force) = 0; - virtual void cal_Stress(ModuleBase::matrix& stress) = 0; - virtual void postprocess(){}; + //! calcualte forces for the atoms in the given cell + virtual void cal_force(ModuleBase::matrix& force) = 0; + + //! calcualte stress of given cell + virtual void cal_stress(ModuleBase::matrix& stress) = 0; + + //! perform post processing calculations + virtual void post_process(){}; // Print current classname. void printname(); @@ -59,7 +60,7 @@ class ESolver * * @return [out] std::string The type of ESolver */ -std::string determine_type(); +std::string determine_type(void); /** * @brief Determine and initialize an ESolver based on input information. diff --git a/source/module_esolver/esolver_dp.cpp b/source/module_esolver/esolver_dp.cpp index 8551ead5ff..4406a38ef8 100644 --- a/source/module_esolver/esolver_dp.cpp +++ b/source/module_esolver/esolver_dp.cpp @@ -25,7 +25,7 @@ namespace ModuleESolver { - void ESolver_DP::Init(Input& inp, UnitCell& ucell) + void ESolver_DP::init(Input& inp, UnitCell& ucell) { ucell_ = &ucell; dp_potential = 0; @@ -59,7 +59,7 @@ namespace ModuleESolver assert(ucell.nat == iat); } - void ESolver_DP::Run(const int istep, UnitCell& ucell) + void ESolver_DP::run(const int istep, UnitCell& ucell) { ModuleBase::TITLE("ESolver_DP", "Run"); ModuleBase::timer::tick("ESolver_DP", "Run"); @@ -122,18 +122,18 @@ namespace ModuleESolver ModuleBase::timer::tick("ESolver_DP", "Run"); } - double ESolver_DP::cal_Energy() + double ESolver_DP::cal_energy() { return dp_potential; } - void ESolver_DP::cal_Force(ModuleBase::matrix& force) + void ESolver_DP::cal_force(ModuleBase::matrix& force) { force = dp_force; ModuleIO::print_force(GlobalV::ofs_running, *ucell_, "TOTAL-FORCE (eV/Angstrom)", force, false); } - void ESolver_DP::cal_Stress(ModuleBase::matrix& stress) + void ESolver_DP::cal_stress(ModuleBase::matrix& stress) { stress = dp_virial; @@ -148,7 +148,7 @@ namespace ModuleESolver ModuleIO::print_stress("TOTAL-STRESS", stress, true, false); } - void ESolver_DP::postprocess() + void ESolver_DP::post_process(void) { GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; GlobalV::ofs_running << std::setprecision(16); diff --git a/source/module_esolver/esolver_dp.h b/source/module_esolver/esolver_dp.h index 8e61661c38..54008eb038 100644 --- a/source/module_esolver/esolver_dp.h +++ b/source/module_esolver/esolver_dp.h @@ -36,7 +36,7 @@ class ESolver_DP : public ESolver * @param inp input parameters * @param cell unitcell information */ - void Init(Input& inp, UnitCell& cell) override; + void init(Input& inp, UnitCell& cell) override; /** * @brief Run the DP solver for a given ion/md step and unit cell @@ -44,7 +44,7 @@ class ESolver_DP : public ESolver * @param istep the current ion/md step * @param cell unitcell information */ - void Run(const int istep, UnitCell& cell) override; + void run(const int istep, UnitCell& cell) override; /** * @brief get the total energy without ion kinetic energy @@ -52,28 +52,28 @@ class ESolver_DP : public ESolver * @param etot the computed energy * @return total energy without ion kinetic energy */ - double cal_Energy() override; + double cal_energy() override; /** * @brief get the computed atomic forces * * @param force the computed atomic forces */ - void cal_Force(ModuleBase::matrix& force) override; + void cal_force(ModuleBase::matrix& force) override; /** * @brief get the computed lattice virials * * @param stress the computed lattice virials */ - void cal_Stress(ModuleBase::matrix& stress) override; + void cal_stress(ModuleBase::matrix& stress) override; /** * @brief Prints the final total energy of the DP model to the output file * * This function prints the final total energy of the DP model in eV to the output file along with some formatting. */ - void postprocess() override; + void post_process() override; private: /** diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index 566c97b877..076196941e 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -40,7 +40,7 @@ ESolver_FP::~ESolver_FP() } -void ESolver_FP::Init(Input& inp, UnitCell& cell) +void ESolver_FP::init(Input& inp, UnitCell& cell) { if(!GlobalV::use_paw) { diff --git a/source/module_esolver/esolver_fp.h b/source/module_esolver/esolver_fp.h index b05c400245..38e492d581 100644 --- a/source/module_esolver/esolver_fp.h +++ b/source/module_esolver/esolver_fp.h @@ -40,7 +40,7 @@ namespace ModuleESolver virtual ~ESolver_FP(); //! Initialize of the first-principels energy solver - virtual void Init(Input& inp, UnitCell& cell) override; + virtual void init(Input& inp, UnitCell& cell) override; virtual void init_after_vc(Input& inp, UnitCell& cell); // liuyu add 2023-03-09 diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index f0ec3d49c2..419dc48ca3 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -73,9 +73,9 @@ ESolver_KS::~ESolver_KS() } template -void ESolver_KS::Init(Input& inp, UnitCell& ucell) +void ESolver_KS::init(Input& inp, UnitCell& ucell) { - ESolver_FP::Init(inp,ucell); + ESolver_FP::init(inp,ucell); //------------------Charge Mixing------------------ p_chgmix->set_mixing(GlobalV::MIXING_MODE, @@ -94,12 +94,12 @@ void ESolver_KS::Init(Input& inp, UnitCell& ucell) #ifdef USE_PAW if(GlobalV::use_paw) { - int * atom_type; - double ** atom_coord; + int * atom_type = nullptr; + double ** atom_coord = nullptr; std::vector filename_list; atom_type = new int [ucell.nat]; - atom_coord = new double * [ucell.nat]; + atom_coord = new double *[ucell.nat]; filename_list.resize(ucell.ntype); for(int ia = 0; ia < ucell.nat; ia ++) @@ -199,12 +199,14 @@ void ESolver_KS::Init(Input& inp, UnitCell& ucell) #ifdef __MPI this->pw_wfc->initmpi(GlobalV::NPROC_IN_POOL, GlobalV::RANK_IN_POOL, POOL_WORLD); #endif + this->pw_wfc->initgrids(inp.ref_cell_factor * ucell.lat0, ucell.latvec, this->pw_rho->nx, this->pw_rho->ny, this->pw_rho->nz); this->pw_wfc->initparameters(false, inp.ecutwfc, this->kv.nks, this->kv.kvec_d.data()); + #ifdef __MPI if (INPUT.pw_seed > 0) { @@ -212,8 +214,11 @@ void ESolver_KS::Init(Input& inp, UnitCell& ucell) } // qianrui add 2021-8-13 to make different kpar parameters can get the same results #endif + this->pw_wfc->ft.fft_mode = inp.fft_mode; + this->pw_wfc->setuptransform(); + for (int ik = 0; ik < this->kv.nks; ++ik) { this->kv.ngk[ik] = this->pw_wfc->npwk[ik]; @@ -285,6 +290,7 @@ void ESolver_KS::Init(Input& inp, UnitCell& ucell) #endif } + template void ESolver_KS::init_after_vc(Input& inp, UnitCell& ucell) { @@ -322,6 +328,7 @@ void ESolver_KS::hamilt2density(const int istep, const int iter, cons ModuleBase::timer::tick(this->classname, "hamilt2density"); } + template void ESolver_KS::print_wfcfft(Input& inp, std::ofstream &ofs) { @@ -352,6 +359,7 @@ void ESolver_KS::print_wfcfft(Input& inp, std::ofstream &ofs) ofs << "\n PARALLEL PW FOR WAVE FUNCTIONS" << std::endl; ofs <<" "<< std::setw(8) << "PROC"<< std::setw(15) << "COLUMNS(POT)"<< std::setw(15) << "PW" << std::endl; + for (int i = 0; i < GlobalV::NPROC_IN_POOL ; ++i) { ofs <<" "<::print_wfcfft(Input& inp, std::ofstream &ofs) << std::setw(15) << this->pw_wfc->npw_per[i] << std::endl; } + ofs << " --------------- sum -------------------" << std::endl; ofs << " " << std::setw(8) << GlobalV::NPROC_IN_POOL @@ -367,25 +376,26 @@ void ESolver_KS::print_wfcfft(Input& inp, std::ofstream &ofs) ModuleBase::GlobalFunc::DONE(ofs, "INIT PLANEWAVE"); } + template -void ESolver_KS::Run(const int istep, UnitCell& ucell) +void ESolver_KS::run(const int istep, UnitCell& ucell) { if (!(GlobalV::CALCULATION == "scf" || GlobalV::CALCULATION == "md" || GlobalV::CALCULATION == "relax" || GlobalV::CALCULATION == "cell-relax")) { - this->othercalculation(istep); + this->others(istep); } else { - ModuleBase::timer::tick(this->classname, "Run"); + ModuleBase::timer::tick(this->classname, "run"); - this->beforescf(istep); //Something else to do before the iter loop + this->before_scf(istep); //Something else to do before the iter loop ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT SCF"); if(this->maxniter > 0) { - this->printhead(); //print the headline on the screen. + this->print_head(); //print the headline on the screen. } bool firstscf = true; @@ -393,7 +403,7 @@ void ESolver_KS::Run(const int istep, UnitCell& ucell) this->niter = this->maxniter; for (int iter = 1; iter <= this->maxniter; ++iter) { - writehead(GlobalV::ofs_running, istep, iter); + this->write_head(GlobalV::ofs_running, istep, iter); #ifdef __MPI auto iterstart = MPI_Wtime(); #else @@ -401,7 +411,7 @@ void ESolver_KS::Run(const int istep, UnitCell& ucell) #endif double diag_ethr = this->phsol->set_diagethr(istep, iter, drho); - eachiterinit(istep, iter); + this->iter_init(istep, iter); this->hamilt2density(istep, iter, diag_ethr); @@ -480,12 +490,14 @@ void ESolver_KS::Run(const int istep, UnitCell& ucell) // Hamilt should be used after it is constructed. // this->phamilt->update(conv_elec); - updatepot(istep, iter); - eachiterfinish(iter); + this->update_pot(istep, iter); + this->iter_finish(iter); #ifdef __MPI double duration = (double)(MPI_Wtime() - iterstart); #else - double duration = (std::chrono::duration_cast(std::chrono::system_clock::now() - iterstart)).count() / static_cast(1e6); + double duration = + (std::chrono::duration_cast(std::chrono::system_clock::now() + - iterstart)).count() / static_cast(1e6); #endif /* SCF print: G1 -3.435545e+03 0.000000e+00 3.607e-01 2.862e-01 @@ -496,7 +508,7 @@ void ESolver_KS::Run(const int istep, UnitCell& ucell) { dkin = p_chgmix->get_dkin(pelec->charge, GlobalV::nelec); } - printiter(iter, drho, dkin, duration, diag_ethr); + this->print_iter(iter, drho, dkin, duration, diag_ethr); #ifdef __RAPIDJSON //add Json of scf mag @@ -532,8 +544,8 @@ void ESolver_KS::Run(const int istep, UnitCell& ucell) this->conv_elec ); #endif //__RAPIDJSON - afterscf(istep); - ModuleBase::timer::tick(this->classname, "Run"); + after_scf(istep); + ModuleBase::timer::tick(this->classname, "run"); } #ifdef __RAPIDJSON @@ -545,27 +557,33 @@ void ESolver_KS::Run(const int istep, UnitCell& ucell) return; }; + template -void ESolver_KS::printhead() +void ESolver_KS::print_head(void) { std::cout << " " << std::setw(7) << "ITER"; + if (GlobalV::NSPIN == 2) { std::cout << std::setw(10) << "TMAG"; std::cout << std::setw(10) << "AMAG"; } + std::cout << std::setw(15) << "ETOT(eV)"; std::cout << std::setw(15) << "EDIFF(eV)"; std::cout << std::setw(11) << "DRHO"; + if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) { std::cout << std::setw(11) << "DKIN"; } + std::cout << std::setw(11) << "TIME(s)" << std::endl; } + template -void ESolver_KS::printiter( +void ESolver_KS::print_iter( const int iter, const double drho, const double dkin, @@ -575,8 +593,9 @@ void ESolver_KS::printiter( this->pelec->print_etot(this->conv_elec, iter, drho, dkin, duration, INPUT.printe, ethr); } + template -void ESolver_KS::writehead(std::ofstream& ofs_running, const int istep, const int iter) +void ESolver_KS::write_head(std::ofstream& ofs_running, const int istep, const int iter) { ofs_running << "\n " @@ -586,6 +605,7 @@ void ESolver_KS::writehead(std::ofstream& ofs_running, const int iste << "--------------------------------\n"; } + template int ESolver_KS::getniter() { @@ -599,7 +619,7 @@ ModuleIO::Output_Rho ESolver_KS::create_Output_Rho( int iter, const std::string& prefix) { - int precision = 3; + const int precision = 3; std::string tag = "CHG"; return ModuleIO::Output_Rho(this->pw_big, this->pw_rhod, @@ -618,7 +638,7 @@ ModuleIO::Output_Rho ESolver_KS::create_Output_Rho( template ModuleIO::Output_Rho ESolver_KS::create_Output_Kin(int is, int iter, const std::string& prefix) { - int precision = 11; + const int precision = 11; std::string tag = "TAU"; return ModuleIO::Output_Rho(this->pw_big, this->pw_rhod, diff --git a/source/module_esolver/esolver_ks.h b/source/module_esolver/esolver_ks.h index 9b2e898ce8..63f79a60d8 100644 --- a/source/module_esolver/esolver_ks.h +++ b/source/module_esolver/esolver_ks.h @@ -32,17 +32,22 @@ class ESolver_KS : public ESolver_FP virtual ~ESolver_KS(); double scf_thr; // scf threshold + double drho; // the difference between rho_in (before HSolver) and rho_out (After HSolver) + int maxniter; // maximum iter steps for scf + int niter; // iter steps actually used in scf + bool conv_elec; // If electron density is converged in scf. + int out_freq_elec;// frequency for output - virtual void Init(Input& inp, UnitCell& cell) override; + virtual void init(Input& inp, UnitCell& cell) override; virtual void init_after_vc(Input& inp, UnitCell& cell) override; // liuyu add 2023-03-09 - virtual void Run(const int istep, UnitCell& cell) override; + virtual void run(const int istep, UnitCell& cell) override; // calculate electron density from a specific Hamiltonian virtual void hamilt2density(const int istep, const int iter, const double ethr); @@ -55,19 +60,19 @@ class ESolver_KS : public ESolver_FP protected: //! Something to do before SCF iterations. - virtual void beforescf(int istep) {}; + virtual void before_scf(int istep) {}; //! Something to do before hamilt2density function in each iter loop. - virtual void eachiterinit(const int istep, const int iter) {}; + virtual void iter_init(const int istep, const int iter) {}; //! Something to do after hamilt2density function in each iter loop. - virtual void eachiterfinish(const int iter) {}; + virtual void iter_finish(const int iter) {}; //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. - virtual void afterscf(const int istep) {}; + virtual void after_scf(const int istep) {}; //! It should be replaced by a function in Hamilt Class - virtual void updatepot(const int istep, const int iter) {}; + virtual void update_pot(const int istep, const int iter) {}; //! choose strategy when charge density convergence achieved virtual bool do_after_converge(int& iter){return true;} @@ -76,19 +81,27 @@ class ESolver_KS : public ESolver_FP // Print the headline on the screen: // ITER ETOT(eV) EDIFF(eV) DRHO TIME(s) - void printhead(); + void print_head(); // Print inforamtion in each iter // G1 -3.435545e+03 0.000000e+00 3.607e-01 2.862e-01 // for metaGGA // ITER ETOT(eV) EDIFF(eV) DRHO DKIN TIME(s) // G1 -3.435545e+03 0.000000e+00 3.607e-01 3.522e-01 2.862e-01 - void printiter(const int iter, const double drho, const double dkin, const double duration, const double ethr); + void print_iter( + const int iter, + const double drho, + const double dkin, + const double duration, + const double ethr); // Write the headline in the running_log file // "PW/LCAO" ALGORITHM --------------- ION= 1 ELEC= 1-------------------------------- - void writehead(std::ofstream& ofs_running, const int istep, const int iter); + void write_head( + std::ofstream& ofs_running, + const int istep, + const int iter); /// @brief create a new ModuleIO::Output_Rho object to output charge density ModuleIO::Output_Rho create_Output_Rho(int is, int iter, const std::string& prefix="None"); @@ -99,8 +112,6 @@ class ESolver_KS : public ESolver_FP /// @brief create a new ModuleIO::Output_Potential object to print potential ModuleIO::Output_Potential create_Output_Potential(int iter, const std::string& prefix = "None"); - // TODO: control single precision at input files - //! Solve Hamitonian hsolver::HSolver* phsol = nullptr; diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index c3eb8931ba..03433f9198 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -78,9 +78,9 @@ ESolver_KS_LCAO::~ESolver_KS_LCAO() template -void ESolver_KS_LCAO::Init(Input& inp, UnitCell& ucell) +void ESolver_KS_LCAO::init(Input& inp, UnitCell& ucell) { - ModuleBase::TITLE("ESolver_KS_LCAO", "Init"); + ModuleBase::TITLE("ESolver_KS_LCAO", "init"); // if we are only calculating S, then there is no need // to prepare for potentials and so on @@ -102,7 +102,7 @@ void ESolver_KS_LCAO::Init(Input& inp, UnitCell& ucell) } else { - ESolver_KS::Init(inp, ucell); + ESolver_KS::init(inp, ucell); } // end ifnot get_S // init ElecState @@ -233,6 +233,7 @@ void ESolver_KS_LCAO::Init(Input& inp, UnitCell& ucell) } } + template void ESolver_KS_LCAO::init_after_vc(Input& inp, UnitCell& ucell) { @@ -273,14 +274,14 @@ void ESolver_KS_LCAO::init_after_vc(Input& inp, UnitCell& ucell) template -double ESolver_KS_LCAO::cal_Energy() +double ESolver_KS_LCAO::cal_energy() { return this->pelec->f_en.etot; } template -void ESolver_KS_LCAO::cal_Force(ModuleBase::matrix& force) +void ESolver_KS_LCAO::cal_force(ModuleBase::matrix& force) { Force_Stress_LCAO FSL(this->RA, GlobalC::ucell.nat); FSL.getForceStress(GlobalV::CAL_FORCE, @@ -302,7 +303,7 @@ void ESolver_KS_LCAO::cal_Force(ModuleBase::matrix& force) #endif & GlobalC::ucell.symm); - // delete RA after cal_Force + // delete RA after cal_force this->RA.delete_grid(); @@ -311,12 +312,12 @@ void ESolver_KS_LCAO::cal_Force(ModuleBase::matrix& force) template -void ESolver_KS_LCAO::cal_Stress(ModuleBase::matrix& stress) +void ESolver_KS_LCAO::cal_stress(ModuleBase::matrix& stress) { if (!this->have_force) { ModuleBase::matrix fcs; - this->cal_Force(fcs); + this->cal_force(fcs); } stress = this->scs; // copy the stress this->have_force = false; @@ -324,7 +325,7 @@ void ESolver_KS_LCAO::cal_Stress(ModuleBase::matrix& stress) template -void ESolver_KS_LCAO::postprocess() +void ESolver_KS_LCAO::post_process(void) { GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; GlobalV::ofs_running << std::setprecision(16); @@ -416,7 +417,10 @@ void ESolver_KS_LCAO::postprocess() template -void ESolver_KS_LCAO::Init_Basis_lcao(ORB_control& orb_con, Input& inp, UnitCell& ucell) +void ESolver_KS_LCAO::Init_Basis_lcao( + ORB_control& orb_con, + Input& inp, + UnitCell& ucell) { // autoset NB2D first if (GlobalV::NB2D == 0) @@ -457,22 +461,6 @@ void ESolver_KS_LCAO::Init_Basis_lcao(ORB_control& orb_con, Input& inp, two_center_bundle->build_orb_onsite(ucell.ntype, GlobalV::onsite_radius); // currently deepks only use one descriptor file, so cast bool to int is fine - //this->orb_con.read_orb_first(GlobalV::ofs_running, - // GlobalC::ORB, - // ucell.ntype, - // GlobalV::global_orbital_dir, - // ucell.orbital_fn, - // ucell.descriptor_file, - // ucell.lmax, - // inp.lcao_ecut, - // inp.lcao_dk, - // inp.lcao_dr, - // inp.lcao_rmax, - // GlobalV::deepks_setorb, - // inp.out_mat_r, - // GlobalV::CAL_FORCE, - // GlobalV::MY_RANK); - // TODO Due to the omnipresence of GlobalC::ORB, we still have to rely // on the old interface for now. two_center_bundle->to_LCAO_Orbitals(GlobalC::ORB, @@ -516,7 +504,7 @@ void ESolver_KS_LCAO::Init_Basis_lcao(ORB_control& orb_con, Input& inp, template -void ESolver_KS_LCAO::eachiterinit(const int istep, const int iter) +void ESolver_KS_LCAO::iter_init(const int istep, const int iter) { if (iter == 1) { @@ -524,7 +512,8 @@ void ESolver_KS_LCAO::eachiterinit(const int istep, const int iter) this->p_chgmix->mixing_restart = GlobalV::SCF_NMAX + 1; } // for mixing restart - if (iter == this->p_chgmix->mixing_restart && GlobalV::MIXING_RESTART > 0.0) + if (iter == this->p_chgmix->mixing_restart + && GlobalV::MIXING_RESTART > 0.0) { this->p_chgmix->init_mixing(); if (GlobalV::MIXING_DMR) // for mixing_dmr @@ -543,8 +532,9 @@ void ESolver_KS_LCAO::eachiterinit(const int istep, const int iter) // mohan move it outside 2011-01-13 // first need to calculate the weight according to // electrons number. - - if (istep == 0 && this->wf.init_wfc == "file" && this->LOWF.error == 0) + if (istep == 0 + && this->wf.init_wfc == "file" + && this->LOWF.error == 0) { if (iter == 1) { @@ -583,7 +573,11 @@ void ESolver_KS_LCAO::eachiterinit(const int istep, const int iter) { GlobalC::ucell.cal_ux(); } + + //! update the potentials by using new electron charge density this->pelec->pot->update_from_charge(this->pelec->charge, &GlobalC::ucell); + + //! compute the correction energy for metals this->pelec->f_en.descf = this->pelec->cal_delta_escf(); } } @@ -612,13 +606,15 @@ void ESolver_KS_LCAO::eachiterinit(const int istep, const int iter) dynamic_cast*>(this->pelec) ->get_DM(); } - GlobalC::dftu.cal_slater_UJ(this->pelec->charge->rho, this->pw_rho->nrxx); // Calculate U and J if Yukawa potential is used + // Calculate U and J if Yukawa potential is used + GlobalC::dftu.cal_slater_UJ(this->pelec->charge->rho, this->pw_rho->nrxx); } #ifdef __DEEPKS // the density matrixes of DeePKS have been updated in each iter GlobalC::ld.set_hr_cal(true); + // HR in HamiltLCAO should be recalculate if(GlobalV::deepks_scf) { @@ -758,7 +754,7 @@ void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) template -void ESolver_KS_LCAO::updatepot(const int istep, const int iter) +void ESolver_KS_LCAO::update_pot(const int istep, const int iter) { // print Hamiltonian and Overlap matrix if (this->conv_elec) @@ -854,7 +850,7 @@ void ESolver_KS_LCAO::updatepot(const int istep, const int iter) template -void ESolver_KS_LCAO::eachiterfinish(int iter) +void ESolver_KS_LCAO::iter_finish(int iter) { // mix density matrix if (GlobalV::MIXING_RESTART > 0 && iter >= this->p_chgmix->mixing_restart && GlobalV::MIXING_DMR ) @@ -884,11 +880,17 @@ void ESolver_KS_LCAO::eachiterfinish(int iter) for (int ik = 0;ik < this->kv.nks;++ik) { ModuleBase::GlobalFunc::ZEROS(Hexxk_save.data(), Hexxk_save.size()); + hamilt::OperatorEXX> opexx_save(&this->LM, nullptr, &Hexxk_save, this->kv); + opexx_save.contributeHk(ik); + GlobalC::restart.save_disk("Hexx", ik, this->LOWF.ParaV->get_local_size(), Hexxk_save.data()); } - if (GlobalV::MY_RANK == 0)GlobalC::restart.save_disk("Eexx", 0, 1, &this->pelec->f_en.exx); + if (GlobalV::MY_RANK == 0) + { + GlobalC::restart.save_disk("Eexx", 0, 1, &this->pelec->f_en.exx); + } } #endif //----------------------------------- @@ -926,7 +928,7 @@ void ESolver_KS_LCAO::eachiterfinish(int iter) template -void ESolver_KS_LCAO::afterscf(const int istep) +void ESolver_KS_LCAO::after_scf(const int istep) { // save charge difference into files for charge extrapolation if (GlobalV::CALCULATION != "scf") diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index a6e2a20baa..269bb2b3a1 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -16,6 +16,7 @@ #include "module_io/output_mat_sparse.h" #include "module_basis/module_nao/two_center_bundle.h" #include + namespace ModuleESolver { template @@ -25,32 +26,58 @@ namespace ModuleESolver ESolver_KS_LCAO(); ~ESolver_KS_LCAO(); - void Init(Input& inp, UnitCell& cell) override; + void init(Input& inp, UnitCell& cell) override; + void init_after_vc(Input& inp, UnitCell& cell) override; - double cal_Energy() override; - void cal_Force(ModuleBase::matrix& force) override; - void cal_Stress(ModuleBase::matrix& stress) override; - void postprocess() override; + double cal_energy() override; + + void cal_force(ModuleBase::matrix& force) override; + + void cal_stress(ModuleBase::matrix& stress) override; + + void post_process() override; + void nscf() override; + void get_S(); protected: - virtual void beforescf(const int istep) override; - virtual void eachiterinit(const int istep, const int iter) override; + + virtual void before_scf(const int istep) override; + + virtual void iter_init(const int istep, const int iter) override; + virtual void hamilt2density(const int istep, const int iter, const double ethr) override; - virtual void updatepot(const int istep, const int iter) override; - virtual void eachiterfinish(const int iter) override; - virtual void afterscf(const int istep) override; + + virtual void update_pot(const int istep, const int iter) override; + + virtual void iter_finish(const int iter) override; + + virtual void after_scf(const int istep) override; + virtual bool do_after_converge(int& iter) override; - virtual void othercalculation(const int istep)override; + virtual void others(const int istep)override; + + // we will get rid of this class soon, don't use it, mohan 2024-03-28 ORB_control orb_con; //Basis_LCAO + + // we will get rid of this class soon, don't use it, mohan 2024-03-28 Record_adj RA; + + // we will get rid of this class soon, don't use it, mohan 2024-03-28 Local_Orbital_wfc LOWF; + + // we will get rid of this class soon, don't use it, mohan 2024-03-28 Local_Orbital_Charge LOC; + + // we will get rid of this class soon, don't use it, mohan 2024-03-28 LCAO_Hamilt UHM; + + // we will get rid of this class soon, don't use it, mohan 2024-03-28 LCAO_Matrix LM; + Grid_Technique GridT; std::unique_ptr two_center_bundle; @@ -97,8 +124,5 @@ namespace ModuleESolver #endif }; - - - } #endif diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp index e16774c6f1..57a0bcd3ea 100644 --- a/source/module_esolver/esolver_ks_lcao_elec.cpp +++ b/source/module_esolver/esolver_ks_lcao_elec.cpp @@ -273,10 +273,10 @@ void ESolver_KS_LCAO::beforesolver(const int istep) } template -void ESolver_KS_LCAO::beforescf(int istep) +void ESolver_KS_LCAO::before_scf(int istep) { - ModuleBase::TITLE("ESolver_KS_LCAO", "beforescf"); - ModuleBase::timer::tick("ESolver_KS_LCAO", "beforescf"); + ModuleBase::TITLE("ESolver_KS_LCAO", "before_scf"); + ModuleBase::timer::tick("ESolver_KS_LCAO", "before_scf"); if (GlobalC::ucell.cell_parameter_updated) { @@ -335,15 +335,15 @@ void ESolver_KS_LCAO::beforescf(int istep) this->p_hamilt->non_first_scf = istep; - ModuleBase::timer::tick("ESolver_KS_LCAO", "beforescf"); + ModuleBase::timer::tick("ESolver_KS_LCAO", "before_scf"); return; } template -void ESolver_KS_LCAO::othercalculation(const int istep) +void ESolver_KS_LCAO::others(const int istep) { - ModuleBase::TITLE("ESolver_KS_LCAO", "othercalculation"); - ModuleBase::timer::tick("ESolver_KS_LCAO", "othercalculation"); + ModuleBase::TITLE("ESolver_KS_LCAO", "others"); + ModuleBase::timer::tick("ESolver_KS_LCAO", "others"); if (GlobalV::CALCULATION == "get_S") { this->get_S(); @@ -440,10 +440,10 @@ void ESolver_KS_LCAO::othercalculation(const int istep) } else { - ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::othercalculation", "CALCULATION type not supported"); + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::others", "CALCULATION type not supported"); } - ModuleBase::timer::tick("ESolver_KS_LCAO", "othercalculation"); + ModuleBase::timer::tick("ESolver_KS_LCAO", "others"); return; } template <> diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index ab2ad31acc..a4ff74d150 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -62,9 +62,9 @@ ESolver_KS_LCAO_TDDFT::~ESolver_KS_LCAO_TDDFT() } } -void ESolver_KS_LCAO_TDDFT::Init(Input& inp, UnitCell& ucell) +void ESolver_KS_LCAO_TDDFT::init(Input& inp, UnitCell& ucell) { - ESolver_KS::Init(inp, ucell); + ESolver_KS::init(inp, ucell); // Initialize the FFT. // this function belongs to cell LOOP @@ -257,7 +257,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density( this->pelec->f_en.deband = this->pelec->cal_delta_eband(); } -void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter) +void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) { // print Hamiltonian and Overlap matrix if (this->conv_elec) @@ -434,7 +434,7 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter) } -void ESolver_KS_LCAO_TDDFT::afterscf(const int istep) +void ESolver_KS_LCAO_TDDFT::after_scf(const int istep) { for (int is = 0; is < GlobalV::NSPIN; is++) { @@ -457,9 +457,10 @@ void ESolver_KS_LCAO_TDDFT::afterscf(const int istep) this->RA, this->UHM); } - ESolver_KS_LCAO, double>::afterscf(istep); + ESolver_KS_LCAO, double>::after_scf(istep); } + // use the original formula (Hamiltonian matrix) to calculate energy density matrix void ESolver_KS_LCAO_TDDFT::cal_edm_tddft(void) { diff --git a/source/module_esolver/esolver_ks_lcao_tddft.h b/source/module_esolver/esolver_ks_lcao_tddft.h index eaf785ae92..8e1f4a2c89 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.h +++ b/source/module_esolver/esolver_ks_lcao_tddft.h @@ -17,21 +17,32 @@ namespace ModuleESolver class ESolver_KS_LCAO_TDDFT : public ESolver_KS_LCAO, double> { public: + ESolver_KS_LCAO_TDDFT(); + ~ESolver_KS_LCAO_TDDFT(); - void Init(Input& inp, UnitCell& cell) override; + + void init(Input& inp, UnitCell& cell) override; psi::Psi>* psi_laststep = nullptr; + std::complex** Hk_laststep = nullptr; + std::complex** Sk_laststep = nullptr; + //same as pelec elecstate::ElecStateLCAO_TDDFT* pelec_td = nullptr; + int td_htype = 1; protected: + virtual void hamilt2density(const int istep, const int iter, const double ethr) override; - virtual void updatepot(const int istep, const int iter) override; - virtual void afterscf(const int istep) override; + + virtual void update_pot(const int istep, const int iter) override; + + virtual void after_scf(const int istep) override; + void cal_edm_tddft(); }; diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 0395ff660d..8e1054f6b0 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -105,6 +105,8 @@ ESolver_KS_PW::~ESolver_KS_PW() } delete this->psi; } + + template void ESolver_KS_PW::Init_GlobalC(Input& inp, UnitCell& cell) { @@ -217,9 +219,9 @@ void ESolver_KS_PW::Init_GlobalC(Input& inp, UnitCell& cell) template -void ESolver_KS_PW::Init(Input& inp, UnitCell& ucell) +void ESolver_KS_PW::init(Input& inp, UnitCell& ucell) { - ESolver_KS::Init(inp, ucell); + ESolver_KS::init(inp, ucell); // init HSolver if (this->phsol == nullptr) @@ -433,10 +435,11 @@ void ESolver_KS_PW::init_after_vc(Input& inp, UnitCell& ucell) ModuleBase::timer::tick("ESolver_KS_PW", "init_after_vc"); } + template -void ESolver_KS_PW::beforescf(int istep) +void ESolver_KS_PW::before_scf(int istep) { - ModuleBase::TITLE("ESolver_KS_PW", "beforescf"); + ModuleBase::TITLE("ESolver_KS_PW", "before_scf"); if (GlobalC::ucell.cell_parameter_updated) { @@ -523,7 +526,7 @@ void ESolver_KS_PW::beforescf(int istep) // before hamilt2density, we update Hk and initialize psi if(GlobalV::psi_initializer) { - // beforescf function will be called everytime before scf. However, once atomic coordinates changed, + // before_scf function will be called everytime before scf. However, once atomic coordinates changed, // structure factor will change, therefore all atomwise properties will change. So we need to reinitialize // psi every time before scf. But for random wavefunction, we dont, because random wavefunction is not // related to atomic coordinates. @@ -538,10 +541,10 @@ void ESolver_KS_PW::beforescf(int istep) template -void ESolver_KS_PW::othercalculation(const int istep) +void ESolver_KS_PW::others(const int istep) { - ModuleBase::TITLE("ESolver_KS_PW", "othercalculation"); - ModuleBase::timer::tick("ESolver_KS_PW", "othercalculation"); + ModuleBase::TITLE("ESolver_KS_PW", "others"); + ModuleBase::timer::tick("ESolver_KS_PW", "others"); if (GlobalV::CALCULATION == "test_memory") { Cal_Test::test_memory(this->pw_rho, @@ -571,15 +574,15 @@ void ESolver_KS_PW::othercalculation(const int istep) } else { - ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::othercalculation", "CALCULATION type not supported"); + ModuleBase::WARNING_QUIT("ESolver_KS_PW::others", "CALCULATION type not supported"); } - ModuleBase::timer::tick("ESolver_KS_PW", "othercalculation"); + ModuleBase::timer::tick("ESolver_KS_PW", "others"); return; } template -void ESolver_KS_PW::eachiterinit(const int istep, const int iter) +void ESolver_KS_PW::iter_init(const int istep, const int iter) { if (iter == 1) { @@ -804,7 +807,7 @@ void ESolver_KS_PW::hamilt2density( // before hamilt2density, we update Hk and initialize psi if(GlobalV::psi_initializer) { - // beforescf function will be called everytime before scf. However, once atomic coordinates changed, + // before_scf function will be called everytime before scf. However, once atomic coordinates changed, // structure factor will change, therefore all atomwise properties will change. So we need to reinitialize // psi every time before scf. But for random wavefunction, we dont, because random wavefunction is not // related to atomic coordinates. @@ -899,7 +902,7 @@ void ESolver_KS_PW::hamilt2density( // Temporary, it should be rewritten with Hamilt class. template -void ESolver_KS_PW::updatepot(const int istep, const int iter) +void ESolver_KS_PW::update_pot(const int istep, const int iter) { if (!this->conv_elec) { @@ -918,7 +921,7 @@ void ESolver_KS_PW::updatepot(const int istep, const int iter) template -void ESolver_KS_PW::eachiterfinish(const int iter) +void ESolver_KS_PW::iter_finish(const int iter) { // liuyu 2023-10-24 // D in uspp need vloc, thus needs update when veff updated @@ -966,7 +969,7 @@ void ESolver_KS_PW::eachiterfinish(const int iter) template -void ESolver_KS_PW::afterscf(const int istep) +void ESolver_KS_PW::after_scf(const int istep) { this->create_Output_Potential(istep).write(); @@ -1075,19 +1078,22 @@ void ESolver_KS_PW::afterscf(const int istep) template -double ESolver_KS_PW::cal_Energy() +double ESolver_KS_PW::cal_energy() { return this->pelec->f_en.etot; } template -void ESolver_KS_PW::cal_Force(ModuleBase::matrix& force) +void ESolver_KS_PW::cal_force(ModuleBase::matrix& force) { Forces ff(GlobalC::ucell.nat); - if (this->__kspw_psi != nullptr) - this->__kspw_psi = nullptr; - if (this->__kspw_psi == nullptr) + if (this->__kspw_psi != nullptr) + { + this->__kspw_psi = nullptr; + } + + if (this->__kspw_psi == nullptr) { this->__kspw_psi = GlobalV::precision_flag == "single" ? new psi::Psi, Device>(this->kspw_psi[0]) @@ -1107,11 +1113,14 @@ void ESolver_KS_PW::cal_Force(ModuleBase::matrix& force) template -void ESolver_KS_PW::cal_Stress(ModuleBase::matrix& stress) +void ESolver_KS_PW::cal_stress(ModuleBase::matrix& stress) { Stress_PW ss(this->pelec); if (this->__kspw_psi != nullptr) + { this->__kspw_psi = nullptr; + } + if (this->__kspw_psi == nullptr) { this->__kspw_psi = GlobalV::precision_flag == "single" @@ -1141,7 +1150,7 @@ void ESolver_KS_PW::cal_Stress(ModuleBase::matrix& stress) template -void ESolver_KS_PW::postprocess(void) +void ESolver_KS_PW::post_process(void) { GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; @@ -1361,7 +1370,9 @@ void ESolver_KS_PW::nscf(void) ModuleBase::TITLE("ESolver_KS_PW", "nscf"); ModuleBase::timer::tick("ESolver_KS_PW", "nscf"); - this->beforescf(0); + // mohan add istep_tmp 2024-03-31 + const int istep_tmp = 0; + this->before_scf(istep_tmp); //! Setup the parameters for diagonalization double diag_ethr = GlobalV::PW_DIAG_THR; diff --git a/source/module_esolver/esolver_ks_pw.h b/source/module_esolver/esolver_ks_pw.h index d08a24f552..41aab7892e 100644 --- a/source/module_esolver/esolver_ks_pw.h +++ b/source/module_esolver/esolver_ks_pw.h @@ -28,15 +28,15 @@ class ESolver_KS_PW : public ESolver_KS ~ESolver_KS_PW(); - void Init(Input& inp, UnitCell& cell) override; + void init(Input& inp, UnitCell& cell) override; void init_after_vc(Input& inp, UnitCell& cell) override; - double cal_Energy() override; + double cal_energy() override; - void cal_Force(ModuleBase::matrix& force) override; + void cal_force(ModuleBase::matrix& force) override; - void cal_Stress(ModuleBase::matrix& stress) override; + void cal_stress(ModuleBase::matrix& stress) override; virtual void hamilt2density(const int istep, const int iter, const double ethr) override; @@ -44,7 +44,7 @@ class ESolver_KS_PW : public ESolver_KS virtual void nscf() override; - void postprocess() override; + void post_process() override; /** * @brief calculate Onsager coefficients Lmn(\omega) and conductivities with Kubo-Greenwood formula @@ -88,17 +88,17 @@ class ESolver_KS_PW : public ESolver_KS protected: - virtual void beforescf(const int istep) override; + virtual void before_scf(const int istep) override; - virtual void eachiterinit(const int istep, const int iter) override; + virtual void iter_init(const int istep, const int iter) override; - virtual void updatepot(const int istep, const int iter) override; + virtual void update_pot(const int istep, const int iter) override; - virtual void eachiterfinish(const int iter) override; + virtual void iter_finish(const int iter) override; - virtual void afterscf(const int istep) override; + virtual void after_scf(const int istep) override; - virtual void othercalculation(const int istep)override; + virtual void others(const int istep)override; //temporary, this will be removed in the future; //Init Global class diff --git a/source/module_esolver/esolver_lj.cpp b/source/module_esolver/esolver_lj.cpp index f48ff290d8..0b5346a53f 100644 --- a/source/module_esolver/esolver_lj.cpp +++ b/source/module_esolver/esolver_lj.cpp @@ -7,7 +7,7 @@ namespace ModuleESolver { - void ESolver_LJ::Init(Input& inp, UnitCell& ucell) + void ESolver_LJ::init(Input& inp, UnitCell& ucell) { ucell_ = &ucell; lj_potential = 0; @@ -24,7 +24,7 @@ namespace ModuleESolver lj_sigma *= ModuleBase::ANGSTROM_AU; } - void ESolver_LJ::Run(const int istep, UnitCell& ucell) + void ESolver_LJ::run(const int istep, UnitCell& ucell) { Grid_Driver grid_neigh(GlobalV::test_deconstructor, GlobalV::test_grid_driver, GlobalV::test_grid); atom_arrange::search( @@ -93,18 +93,18 @@ namespace ModuleESolver #endif } - double ESolver_LJ::cal_Energy() + double ESolver_LJ::cal_energy() { return lj_potential; } - void ESolver_LJ::cal_Force(ModuleBase::matrix& force) + void ESolver_LJ::cal_force(ModuleBase::matrix& force) { force = lj_force; ModuleIO::print_force(GlobalV::ofs_running, *ucell_, "TOTAL-FORCE (eV/Angstrom)", force, false); } - void ESolver_LJ::cal_Stress(ModuleBase::matrix& stress) + void ESolver_LJ::cal_stress(ModuleBase::matrix& stress) { stress = lj_virial; @@ -119,7 +119,7 @@ namespace ModuleESolver ModuleIO::print_stress("TOTAL-STRESS", stress, true, false); } - void ESolver_LJ::postprocess() + void ESolver_LJ::post_process(void) { GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; GlobalV::ofs_running << std::setprecision(16); diff --git a/source/module_esolver/esolver_lj.h b/source/module_esolver/esolver_lj.h index 9d04243e1c..6963e4b854 100644 --- a/source/module_esolver/esolver_lj.h +++ b/source/module_esolver/esolver_lj.h @@ -14,17 +14,25 @@ namespace ModuleESolver classname = "ESolver_LJ"; } - void Init(Input& inp, UnitCell& cell) override; - void Run(const int istep, UnitCell& cell) override; - double cal_Energy() override; - void cal_Force(ModuleBase::matrix& force) override; - void cal_Stress(ModuleBase::matrix& stress) override; - void postprocess() override; + void init(Input& inp, UnitCell& cell) override; + + void run(const int istep, UnitCell& cell) override; + + double cal_energy() override; + + void cal_force(ModuleBase::matrix& force) override; + + void cal_stress(ModuleBase::matrix& stress) override; + + void post_process() override; private: + double LJ_energy(const double d); + ModuleBase::Vector3 LJ_force(const double d, const ModuleBase::Vector3 dr); + void LJ_virial(const ModuleBase::Vector3& force, const ModuleBase::Vector3& dtau); diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp index 53e1224e44..ccd6bc3b5b 100644 --- a/source/module_esolver/esolver_of.cpp +++ b/source/module_esolver/esolver_of.cpp @@ -57,9 +57,9 @@ ESolver_OF::~ESolver_OF() delete this->opt_cg_mag_; } -void ESolver_OF::Init(Input& inp, UnitCell& ucell) +void ESolver_OF::init(Input& inp, UnitCell& ucell) { - ESolver_FP::Init(inp, ucell); + ESolver_FP::init(inp, ucell); // save necessary parameters this->of_kinetic_ = inp.of_kinetic; @@ -217,9 +217,9 @@ void ESolver_OF::init_after_vc(Input& inp, UnitCell& ucell) } } -void ESolver_OF::Run(int istep, UnitCell& ucell) +void ESolver_OF::run(int istep, UnitCell& ucell) { - ModuleBase::timer::tick("ESolver_OF", "Run"); + ModuleBase::timer::tick("ESolver_OF", "run"); // get Ewald energy, initial rho and phi if necessary this->before_opt(istep, ucell); this->iter_ = 0; @@ -232,11 +232,13 @@ void ESolver_OF::Run(int istep, UnitCell& ucell) // calculate the energy of new rho and phi this->energy_llast_ = this->energy_last_; this->energy_last_ = this->energy_current_; - this->energy_current_ = this->cal_Energy(); + this->energy_current_ = this->cal_energy(); // check if the job is done - if (this->check_exit()) - break; + if (this->check_exit()) + { + break; + } // find the optimization direction and step lenghth theta according to the potential this->optimize(ucell); @@ -249,7 +251,7 @@ void ESolver_OF::Run(int istep, UnitCell& ucell) this->after_opt(istep, ucell); - ModuleBase::timer::tick("ESolver_OF", "Run"); + ModuleBase::timer::tick("ESolver_OF", "run"); } /** @@ -334,7 +336,10 @@ void ESolver_OF::update_potential(UnitCell& ucell) { // (1) get dL/dphi if (GlobalV::NSPIN == 4) + { ucell.cal_ux(); + } + this->pelec->pot->update_from_charge(pelec->charge, &ucell); // Hartree + XC + external this->kinetic_potential(pelec->charge->rho, this->pphi_, @@ -582,7 +587,7 @@ void ESolver_OF::after_opt(const int istep, UnitCell& ucell) /** * @brief Output the FINAL_ETOT */ -void ESolver_OF::postprocess() +void ESolver_OF::post_process() { GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; @@ -597,7 +602,7 @@ void ESolver_OF::postprocess() * * @return total energy */ -double ESolver_OF::cal_Energy() +double ESolver_OF::cal_energy() { this->pelec->cal_energies(2); double kinetic_energy = this->kinetic_energy(); // kinetic energy @@ -621,7 +626,7 @@ double ESolver_OF::cal_Energy() * * @param [out] force */ -void ESolver_OF::cal_Force(ModuleBase::matrix& force) +void ESolver_OF::cal_force(ModuleBase::matrix& force) { Forces ff(GlobalC::ucell.nat); ff.cal_force(force, *pelec, this->pw_rho, &GlobalC::ucell.symm, &sf); @@ -632,7 +637,7 @@ void ESolver_OF::cal_Force(ModuleBase::matrix& force) * * @param [out] stress */ -void ESolver_OF::cal_Stress(ModuleBase::matrix& stress) +void ESolver_OF::cal_stress(ModuleBase::matrix& stress) { ModuleBase::matrix kinetic_stress_; kinetic_stress_.create(3, 3); @@ -641,4 +646,4 @@ void ESolver_OF::cal_Stress(ModuleBase::matrix& stress) OF_Stress_PW ss(this->pelec, this->pw_rho); ss.cal_stress(stress, kinetic_stress_, GlobalC::ucell, &GlobalC::ucell.symm, &sf, &kv); } -} // namespace ModuleESolver \ No newline at end of file +} // namespace ModuleESolver diff --git a/source/module_esolver/esolver_of.h b/source/module_esolver/esolver_of.h index ab847321be..42722b27b1 100644 --- a/source/module_esolver/esolver_of.h +++ b/source/module_esolver/esolver_of.h @@ -20,14 +20,19 @@ class ESolver_OF : public ESolver_FP ESolver_OF(); ~ESolver_OF(); - virtual void Init(Input& inp, UnitCell& ucell) override; + virtual void init(Input& inp, UnitCell& ucell) override; + virtual void init_after_vc(Input& inp, UnitCell& ucell) override; - virtual void Run(int istep, UnitCell& ucell) override; - virtual void postprocess() override; - virtual double cal_Energy() override; - virtual void cal_Force(ModuleBase::matrix& force) override; - virtual void cal_Stress(ModuleBase::matrix& stress) override; + virtual void run(int istep, UnitCell& ucell) override; + + virtual void post_process() override; + + virtual double cal_energy() override; + + virtual void cal_force(ModuleBase::matrix& force) override; + + virtual void cal_stress(ModuleBase::matrix& stress) override; virtual int getniter() override { @@ -135,4 +140,4 @@ class ESolver_OF : public ESolver_FP }; } // namespace ModuleESolver -#endif \ No newline at end of file +#endif diff --git a/source/module_esolver/esolver_of_tool.cpp b/source/module_esolver/esolver_of_tool.cpp index 6eb398aaeb..df44a329c8 100644 --- a/source/module_esolver/esolver_of_tool.cpp +++ b/source/module_esolver/esolver_of_tool.cpp @@ -291,6 +291,7 @@ void ESolver_OF::adjust_direction() */ void ESolver_OF::check_direction(double* dEdtheta, double** ptemp_phi, UnitCell& ucell) { + assert(GlobalV::NSPIN>0); double* temp_theta = new double[GlobalV::NSPIN]; ModuleBase::GlobalFunc::ZEROS(temp_theta, GlobalV::NSPIN); @@ -508,4 +509,4 @@ void ESolver_OF::print_info() context.center_title(); GlobalV::ofs_running << context.str() << std::endl; } -} // namespace ModuleESolver \ No newline at end of file +} // namespace ModuleESolver diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp index e41e8629b1..342c513014 100644 --- a/source/module_esolver/esolver_sdft_pw.cpp +++ b/source/module_esolver/esolver_sdft_pw.cpp @@ -36,11 +36,12 @@ ESolver_SDFT_PW::~ESolver_SDFT_PW() { } -void ESolver_SDFT_PW::Init(Input& inp, UnitCell& ucell) +void ESolver_SDFT_PW::init(Input& inp, UnitCell& ucell) { this->nche_sto = inp.nche_sto; this->method_sto = inp.method_sto; - ESolver_KS::Init(inp, ucell); + + ESolver_KS::init(inp, ucell); this->pelec = new elecstate::ElecStatePW_SDFT(pw_wfc, &(chr), @@ -110,23 +111,23 @@ void ESolver_SDFT_PW::Init(Input& inp, UnitCell& ucell) } -void ESolver_SDFT_PW::beforescf(const int istep) +void ESolver_SDFT_PW::before_scf(const int istep) { - ESolver_KS_PW::beforescf(istep); + ESolver_KS_PW::before_scf(istep); if (istep > 0 && INPUT.nbands_sto != 0 && INPUT.initsto_freq > 0 && istep % INPUT.initsto_freq == 0) { Update_Sto_Orbitals(this->stowf, INPUT.seed_sto); } } -void ESolver_SDFT_PW::eachiterfinish(int iter) +void ESolver_SDFT_PW::iter_finish(int iter) { // this->pelec->print_eigenvalue(GlobalV::ofs_running); this->pelec->cal_energies(2); } -void ESolver_SDFT_PW::afterscf(const int istep) +void ESolver_SDFT_PW::after_scf(const int istep) { // save charge difference into files for charge extrapolation if (GlobalV::CALCULATION != "scf") @@ -220,13 +221,13 @@ void ESolver_SDFT_PW::hamilt2density(int istep, int iter, double ethr) } } -double ESolver_SDFT_PW::cal_Energy() +double ESolver_SDFT_PW::cal_energy() { return this->pelec->f_en.etot; } -void ESolver_SDFT_PW::cal_Force(ModuleBase::matrix& force) +void ESolver_SDFT_PW::cal_force(ModuleBase::matrix& force) { Sto_Forces ff(GlobalC::ucell.nat); @@ -242,7 +243,7 @@ void ESolver_SDFT_PW::cal_Force(ModuleBase::matrix& force) } -void ESolver_SDFT_PW::cal_Stress(ModuleBase::matrix& stress) +void ESolver_SDFT_PW::cal_stress(ModuleBase::matrix& stress) { Sto_Stress_PW ss; ss.cal_stress( @@ -259,7 +260,7 @@ void ESolver_SDFT_PW::cal_Stress(ModuleBase::matrix& stress) } -void ESolver_SDFT_PW::postprocess(void) +void ESolver_SDFT_PW::post_process(void) { GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; @@ -338,10 +339,10 @@ void ESolver_SDFT_PW::postprocess(void) } -void ESolver_SDFT_PW::othercalculation(const int istep) +void ESolver_SDFT_PW::others(const int istep) { - ModuleBase::TITLE("ESolver_SDFT_PW", "othercalculation"); - ModuleBase::timer::tick("ESolver_SDFT_PW", "othercalculation"); + ModuleBase::TITLE("ESolver_SDFT_PW", "others"); + ModuleBase::timer::tick("ESolver_SDFT_PW", "others"); if (GlobalV::CALCULATION == "nscf") { @@ -349,9 +350,9 @@ void ESolver_SDFT_PW::othercalculation(const int istep) } else { - ModuleBase::WARNING_QUIT("ESolver_SDFT_PW::othercalculation", "CALCULATION type not supported"); + ModuleBase::WARNING_QUIT("ESolver_SDFT_PW::others", "CALCULATION type not supported"); } - ModuleBase::timer::tick("ESolver_SDFT_PW", "othercalculation"); + ModuleBase::timer::tick("ESolver_SDFT_PW", "others"); return; } @@ -369,7 +370,7 @@ void ESolver_SDFT_PW::nscf(void) std::cout << " DIGA_THR : " << diag_thr << std::endl; - this->beforescf(istep); + this->before_scf(istep); this->hamilt2density(istep, iter, diag_thr); diff --git a/source/module_esolver/esolver_sdft_pw.h b/source/module_esolver/esolver_sdft_pw.h index 6b3698b3c6..cc9b2aa8b4 100644 --- a/source/module_esolver/esolver_sdft_pw.h +++ b/source/module_esolver/esolver_sdft_pw.h @@ -15,23 +15,33 @@ class ESolver_SDFT_PW : public ESolver_KS_PW> public: ESolver_SDFT_PW(); ~ESolver_SDFT_PW(); - void Init(Input& inp, UnitCell& cell) override; - double cal_Energy() override; - void cal_Force(ModuleBase::matrix& force) override; - void cal_Stress(ModuleBase::matrix& stress) override; + + void init(Input& inp, UnitCell& cell) override; + + double cal_energy() override; + + void cal_force(ModuleBase::matrix& force) override; + + void cal_stress(ModuleBase::matrix& stress) override; public: + Stochastic_WF stowf; protected: - virtual void beforescf(const int istep) override; - // virtual void eachiterinit(int iter) override; + virtual void before_scf(const int istep) override; + virtual void hamilt2density(const int istep, const int iter, const double ethr) override; + virtual void nscf() override; - virtual void othercalculation(const int istep) override; - virtual void eachiterfinish(const int iter) override; - virtual void afterscf(const int istep) override; - virtual void postprocess() override; + + virtual void others(const int istep) override; + + virtual void iter_finish(const int iter) override; + + virtual void after_scf(const int istep) override; + + virtual void post_process() override; public: /** @@ -133,4 +143,4 @@ extern const ModuleBase::matrix* veff; } -#endif \ No newline at end of file +#endif diff --git a/source/module_esolver/test/esolver_dp_test.cpp b/source/module_esolver/test/esolver_dp_test.cpp index 8902c0e10a..8f9fe878a2 100644 --- a/source/module_esolver/test/esolver_dp_test.cpp +++ b/source/module_esolver/test/esolver_dp_test.cpp @@ -10,12 +10,12 @@ /** * - Tested Functions: - * - ESolver_DP::Init() - * - ESolver_DP::Run() - * - ESolver_DP::cal_Energy() - * - ESolver_DP::cal_Force() - * - ESolver_DP::cal_Stress() - * - ESolver_DP::postprocess() + * - ESolver_DP::init() + * - ESolver_DP::run() + * - ESolver_DP::cal_energy() + * - ESolver_DP::cal_force() + * - ESolver_DP::cal_stress() + * - ESolver_DP::post_process() * - ESolver_DP::type_map() */ namespace ModuleIO @@ -39,7 +39,7 @@ class ESolverDPTest : public ::testing::Test { // Initialize variables before each test esolver = new ModuleESolver::ESolver_DP("./support/case_1.pb"); - esolver->Init(inp, ucell); + esolver->init(inp, ucell); } void TearDown() override @@ -85,7 +85,7 @@ TEST_F(ESolverDPTest, InitCase2) esolver->dp_type[0] = 0; esolver->dp_type[1] = 0; esolver->dp_file = "./support/case_2.pb"; - esolver->Init(inp, ucell); + esolver->init(inp, ucell); // Check the initialized variables EXPECT_EQ(esolver->dp_type[0], 0); @@ -100,17 +100,17 @@ TEST_F(ESolverDPTest, RunWarningQuit) int istep = 0; testing::internal::CaptureStdout(); - EXPECT_EXIT(esolver->Run(istep, ucell), ::testing::ExitedWithCode(0), ""); + EXPECT_EXIT(esolver->run(istep, ucell), ::testing::ExitedWithCode(0), ""); std::string output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("Please recompile with -D__DPMD")); } -// Test the cal_Energy() funciton +// Test the cal_energy() funciton TEST_F(ESolverDPTest, CalEnergy) { double etot = 0.0; esolver->dp_potential = 9.8; - etot = esolver->cal_Energy(); + etot = esolver->cal_energy(); // Check the results EXPECT_DOUBLE_EQ(etot, 9.8); @@ -128,7 +128,7 @@ TEST_F(ESolverDPTest, CalForce) } } - esolver->cal_Force(force); + esolver->cal_force(force); // Check the results for (int i = 0; i < ucell.nat; ++i) @@ -152,7 +152,7 @@ TEST_F(ESolverDPTest, CalStress) } } - esolver->cal_Stress(stress); + esolver->cal_stress(stress); // Check the results for (int i = 0; i < 3; ++i) @@ -171,7 +171,7 @@ TEST_F(ESolverDPTest, Postprocess) // Check the results GlobalV::ofs_running.open("log"); - esolver->postprocess(); + esolver->post_process(); GlobalV::ofs_running.close(); std::string expected_output = "\n\n --------------------------------------------\n !FINAL_ETOT_IS 133.3358404 eV\n " @@ -230,4 +230,4 @@ TEST_F(ESolverDPTest, TypeMapWarningQuit) EXPECT_EXIT(esolver->type_map(ucell), ::testing::ExitedWithCode(0), ""); std::string output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("can not find the DP model")); -} \ No newline at end of file +} diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h index 53002629da..90d23c4be3 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h @@ -18,25 +18,34 @@ class LCAO_Hamilt public: LCAO_Hamilt(); + ~LCAO_Hamilt(); void grid_prepare(const Grid_Technique& gt, const ModulePW::PW_Basis& rhopw, const ModulePW::PW_Basis_Big& bigpw); // jingan add 2021-6-4 void set_R_range_sparse(); + void calculate_HContainer_sparse_d(const int ¤t_spin, const double &sparse_threshold, const hamilt::HContainer& hR, std::map, std::map>>& target); + void calculate_HContainer_sparse_cd(const int ¤t_spin, const double &sparse_threshold, const hamilt::HContainer>& hR, std::map, std::map>>>& target); + void calculate_dSTN_R_sparse(const int ¤t_spin, const double &sparse_threshold); + void calculate_STN_R_sparse_for_S(const double &sparse_threshold); + void calculate_STN_R_sparse_for_T(const double &sparse_threshold); + void calculat_HR_dftu_sparse(const int ¤t_spin, const double &sparse_threshold); + void calculat_HR_dftu_soc_sparse(const int ¤t_spin, const double &sparse_threshold); + #ifdef __EXX template void calculate_HR_exx_sparse( const int ¤t_spin, @@ -44,13 +53,25 @@ class LCAO_Hamilt const int (&nmp)[3], const std::vector< std::map >, RI::Tensor > >>& Hexxs); #endif - void calculate_HSR_sparse(const int ¤t_spin, const double &sparse_threshold, const int (&nmp)[3], hamilt::Hamilt>* p_ham); + + void calculate_HSR_sparse( + const int ¤t_spin, + const double &sparse_threshold, + const int (&nmp)[3], + hamilt::Hamilt>* p_ham); + void calculate_SR_sparse(const double &sparse_threshold, hamilt::Hamilt>* p_ham); + void clear_zero_elements(const int ¤t_spin, const double &sparse_threshold); + void destroy_all_HSR_sparse(void); + void calculate_TR_sparse(const double &sparse_threshold); + void destroy_TR_sparse(void); + void calculate_dH_sparse(const int ¤t_spin, const double &sparse_threshold); + void destroy_dH_R_sparse(void); // used for gamma only algorithms. diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h index 6f9d48faf6..f395d051b8 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h @@ -192,41 +192,96 @@ class LCAO_Matrix double* DHloc_fixed_23; double* DHloc_fixed_33; - template - static void set_mat2d(const int& global_ir, const int& global_ic, const T& v, const Parallel_Orbitals& pv, T* mat); - void set_HSgamma(const int& iw1_all, const int& iw2_all, const double& v, double* HSloc); - void set_HSk(const int &iw1_all, const int &iw2_all, const std::complex &v, const char &dtype, const int spin = 0); - - void set_force (const int& iw1_all, const int& iw2_all, const double& vx, const double& vy, - const double& vz, const char &dtype); - void set_stress (const int& iw1_all, const int& iw2_all, const double& vx, const double& vy, - const double& vz, const char &dtype, const ModuleBase::Vector3 &dtau); - - void set_HR_tr(const int &Rx, const int &Ry, const int &Rz, const int &iw1_all, const int &iw2_all, const double &v); - void set_HR_tr_soc(const int &Rx, const int &Ry, const int &Rz, - const int &iw1_all, const int &iw2_all, const std::complex &v); //LiuXh add 2019-07-16 + template + static void set_mat2d( + const int& global_ir, + const int& global_ic, + const T& v, + const Parallel_Orbitals& pv, + T* mat); + + void set_HSgamma( + const int& iw1_all, + const int& iw2_all, + const double& v, + double* HSloc); + + void set_HSk( + const int &iw1_all, + const int &iw2_all, + const std::complex &v, + const char &dtype, + const int spin = 0); + + void set_force ( + const int& iw1_all, + const int& iw2_all, + const double& vx, + const double& vy, + const double& vz, + const char &dtype); + + void set_stress ( + const int& iw1_all, + const int& iw2_all, + const double& vx, + const double& vy, + const double& vz, + const char &dtype, + const ModuleBase::Vector3 &dtau); + + void set_HR_tr( + const int &Rx, + const int &Ry, + const int &Rz, + const int &iw1_all, + const int &iw2_all, + const double &v); + + void set_HR_tr_soc( + const int &Rx, + const int &Ry, + const int &Rz, + const int &iw1_all, + const int &iw2_all, + const std::complex &v); //LiuXh add 2019-07-16 void zeros_HSgamma(const char &mtype); + void zeros_HSk(const char &mtype); + void zeros_HSR(const char &mtype); void print_HSgamma(const char &mtype, std::ostream &os=std::cout); - void print_HSk(const char &mtype, const char &vtype = 'C', const double &accuracy = 1.0e-5, std::ostream &os=std::cout); + + void print_HSk( + const char &mtype, + const char &vtype = 'C', + const double &accuracy = 1.0e-5, + std::ostream &os=std::cout); + void update_Hloc(void); + void update_Hloc2(const int &ik); void allocate_HS_R(const int &nnr); void output_HSk(const char &mtype, std::string &fn); + //LiuXh add 2019-07-15 void allocate_Hloc_fixedR_tr(void); + void allocate_HR_tr(void); + void allocate_SlocR_tr(void); + void destroy_Hloc_fixedR_tr(void); // jingan add 2021-6-4, modify 2021-12-2 void destroy_HS_R_sparse(void); + void destroy_T_R_sparse(void); + void destroy_dH_R_sparse(void); }; diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h b/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h index 1ce3787ae3..50049ee970 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h @@ -13,6 +13,10 @@ #include "module_psi/psi.h" #include "module_elecstate/elecstate.h" #include "module_hamilt_lcao/hamilt_lcaodft/DM_gamma_2d_to_grid.h" + + +// try to get rid of this class, so please do not use it +// in new code anymoer, mohan 2024-03-28 class Local_Orbital_Charge { public: @@ -27,21 +31,24 @@ class Local_Orbital_Charge psi::Psi* psi, const K_Vectors& kv, const int& istep); + void allocate_dm_wfc(const Grid_Technique& gt, elecstate::ElecState* pelec, Local_Orbital_wfc& lowf, psi::Psi>* psi, const K_Vectors& kv, const int& istep); - //----------------- + // in DM_gamma.cpp - //----------------- - void allocate_gamma(const int& lgd, psi::Psi* psid, elecstate::ElecState* pelec, const int& nks, const int& istep); + void allocate_gamma(const int& lgd, + psi::Psi* psid, + elecstate::ElecState* pelec, + const int& nks, + const int& istep); + void cal_dk_gamma_from_2D_pub(void); - //----------------- // in DM_k.cpp - //----------------- void allocate_DM_k(const int& nks, const int& nnrg); // liaochen modify on 2010-3-23 @@ -53,9 +60,7 @@ class Local_Orbital_Charge static int out_dm; // output density matrix or not. static int out_dm1; - //----------------- // in dm_2d.cpp - //----------------- // dm stands for density matrix std::vector dm_gamma; // dm_gamma[is](iw1,iw2); std::vector dm_k; // dm_k[ik](iw1,iw2); @@ -72,13 +77,10 @@ class Local_Orbital_Charge double** dm2d, const K_Vectors& kv); //output, dm2d[NSPIN][LNNR] - //----------------- // wavefunctions' pointer - //----------------- Local_Orbital_wfc* LOWF; - //----------------- + // Parallel Variables' pointer - //----------------- const Parallel_Orbitals* ParaV; //temporary set it to public for ElecStateLCAO class, would be refactor later @@ -87,12 +89,14 @@ class Local_Orbital_Charge std::map, std::map>> DMR_sparse; void set_dm_k(int ik, std::complex* dm_k_in); // set dm_k from a pointer + void set_dm_gamma(int is, double* dm_gamma_in); // set dm_gamma from a pointer private: // whether the DM array has been allocated bool init_DM; + // whether the DM(R) array has been allocated bool init_DM_R; diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.cpp index a242e20a77..cbf4ab56e7 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.cpp @@ -40,7 +40,11 @@ void Local_Orbital_wfc::gamma_file(psi::Psi* psid, elecstate::ElecState* //allocate psi int ncol = this->ParaV->ncol_bands; - if (GlobalV::KS_SOLVER == "genelpa" || GlobalV::KS_SOLVER == "lapack_gvx" || GlobalV::KS_SOLVER == "scalapack_gvx" || GlobalV::KS_SOLVER == "cg_in_lcao" + + if (GlobalV::KS_SOLVER == "genelpa" + || GlobalV::KS_SOLVER == "lapack_gvx" + || GlobalV::KS_SOLVER == "scalapack_gvx" + || GlobalV::KS_SOLVER == "cg_in_lcao" #ifdef __CUSOLVER_LCAO || GlobalV::KS_SOLVER == "cusolver" #endif @@ -48,6 +52,7 @@ void Local_Orbital_wfc::gamma_file(psi::Psi* psid, elecstate::ElecState* { ncol = this->ParaV->ncol; } + if (psid == nullptr) { ModuleBase::WARNING_QUIT("gamma_file", "psid should be allocated first!"); @@ -125,10 +130,8 @@ void Local_Orbital_wfc::allocate_k(const int& lgd, } // allocate the second part. //if(lgd != 0) xiaohui modify 2015-02-04, fixed memory bug - //if(lgd != 0 && this->complex_flag == false) if(lgd != 0) { - //std::cout<<"lgd="<) * GlobalV::NBANDS*GlobalV::NLOCAL); - //ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"MemoryForWaveFunctions (MB)",mem); - //std::cout<<"wfc_k_grid["<complex_flag = true; } } @@ -152,7 +151,10 @@ void Local_Orbital_wfc::allocate_k(const int& lgd, } else if (INPUT.init_wfc == "file") { - if (istep > 0)return; + if (istep > 0) + { + return; + } std::cout << " Read in wave functions files: " << nkstot << std::endl; if (psi == nullptr) { @@ -201,6 +203,8 @@ void Local_Orbital_wfc::allocate_k(const int& lgd, return; } + + int Local_Orbital_wfc::globalIndex(int localindex, int nblk, int nprocs, int myproc) { int iblock, gIndex; @@ -229,9 +233,9 @@ void Local_Orbital_wfc::wfc_2d_to_grid(const int istep, const Parallel_Orbitals* pv = this->ParaV; const int inc = 1; - int myid; + int myid=0; MPI_Comm_rank(pv->comm_2D, &myid); - int info; + int info=0; //calculate maxnloc for bcasting 2d-wfc long maxnloc; // maximum number of elements in local matrix @@ -268,14 +272,18 @@ void Local_Orbital_wfc::wfc_2d_to_grid(const int istep, info=MPI_Bcast(naroc, 2, MPI_INT, src_rank, pv->comm_2D); info=MPI_Bcast(work, maxnloc, MPI_DOUBLE, src_rank, pv->comm_2D); - if (out_wfc_lcao) - info = this->set_wfc_grid(naroc, pv->nb, - pv->dim0, pv->dim1, iprow, ipcol, - work, wfc_grid, myid, ctot); - else - info = this->set_wfc_grid(naroc, pv->nb, - pv->dim0, pv->dim1, iprow, ipcol, - work, wfc_grid); + if (out_wfc_lcao) + { + info = this->set_wfc_grid(naroc, pv->nb, + pv->dim0, pv->dim1, iprow, ipcol, + work, wfc_grid, myid, ctot); + } + else + { + info = this->set_wfc_grid(naroc, pv->nb, + pv->dim0, pv->dim1, iprow, ipcol, + work, wfc_grid); + } }//loop ipcol }//loop iprow @@ -338,12 +346,12 @@ void Local_Orbital_wfc::wfc_2d_to_grid(const int istep, const Parallel_Orbitals* pv = this->ParaV; const int inc = 1; - int myid; + int myid=0; MPI_Comm_rank(pv->comm_2D, &myid); - int info; + int info=0; //calculate maxnloc for bcasting 2d-wfc - long maxnloc; // maximum number of elements in local matrix + long maxnloc=0; // maximum number of elements in local matrix info=MPI_Reduce(&pv->nloc_wfc, &maxnloc, 1, MPI_LONG, MPI_MAX, 0, pv->comm_2D); info=MPI_Bcast(&maxnloc, 1, MPI_LONG, 0, pv->comm_2D); std::complex *work=new std::complex[maxnloc]; // work/buffer matrix @@ -360,7 +368,7 @@ void Local_Orbital_wfc::wfc_2d_to_grid(const int istep, ModuleBase::Memory::record("LOWF::ctot", sizeof(std::complex) * GlobalV::NBANDS * GlobalV::NLOCAL); } - int naroc[2]; // maximum number of row or column + int naroc[2] = {0}; // maximum number of row or column for(int iprow=0; iprowdim0; ++iprow) { for(int ipcol=0; ipcoldim1; ++ipcol) @@ -377,15 +385,19 @@ void Local_Orbital_wfc::wfc_2d_to_grid(const int istep, info=MPI_Bcast(naroc, 2, MPI_INT, src_rank, pv->comm_2D); info = MPI_Bcast(work, maxnloc, MPI_DOUBLE_COMPLEX, src_rank, pv->comm_2D); - if (out_wfc_lcao) - info = this->set_wfc_grid(naroc, pv->nb, - pv->dim0, pv->dim1, iprow, ipcol, - work, wfc_grid, myid, ctot); - else - // mohan update 2021-02-12, delte BFIELD option - info = this->set_wfc_grid(naroc, pv->nb, - pv->dim0, pv->dim1, iprow, ipcol, - work, wfc_grid); + if (out_wfc_lcao) + { + info = this->set_wfc_grid(naroc, pv->nb, + pv->dim0, pv->dim1, iprow, ipcol, + work, wfc_grid, myid, ctot); + } + else + { + // mohan update 2021-02-12, delte BFIELD option + info = this->set_wfc_grid(naroc, pv->nb, + pv->dim0, pv->dim1, iprow, ipcol, + work, wfc_grid); + } }//loop ipcol }//loop iprow @@ -434,4 +446,4 @@ void Local_Orbital_wfc::wfc_2d_to_grid(const int istep, delete[] work; ModuleBase::timer::tick("Local_Orbital_wfc","wfc_2d_to_grid"); } -#endif \ No newline at end of file +#endif diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.h b/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.h index c071d68ccf..d994289029 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.h @@ -81,7 +81,9 @@ class Local_Orbital_wfc #endif int error = 0; + private: + template int set_wfc_grid(int naroc[2], int nb, @@ -101,6 +103,7 @@ class Local_Orbital_wfc }; #ifdef __MPI +// the function should not be defined here!! mohan 2024-03-28 template int Local_Orbital_wfc::set_wfc_grid(int naroc[2], int nb, @@ -115,22 +118,28 @@ int Local_Orbital_wfc::set_wfc_grid(int naroc[2], { ModuleBase::TITLE(" Local_Orbital_wfc", "set_wfc_grid"); if (!wfc && !ctot) + { return 0; + } for (int j = 0; j < naroc[1]; ++j) { int igcol = globalIndex(j, nb, dim1, ipcol); - if (igcol >= GlobalV::NBANDS) - continue; - for (int i = 0; i < naroc[0]; ++i) - { - int igrow = globalIndex(i, nb, dim0, iprow); - int mu_local = this->gridt->trace_lo[igrow]; - if (wfc && mu_local >= 0) - { - wfc[igcol][mu_local] = work[j * naroc[0] + i]; - } - if (ctot && myid == 0) - ctot[igcol][igrow] = work[j * naroc[0] + i]; + if (igcol >= GlobalV::NBANDS) + { + continue; + } + for (int i = 0; i < naroc[0]; ++i) + { + int igrow = globalIndex(i, nb, dim0, iprow); + int mu_local = this->gridt->trace_lo[igrow]; + if (wfc && mu_local >= 0) + { + wfc[igcol][mu_local] = work[j * naroc[0] + i]; + } + if (ctot && myid == 0) + { + ctot[igcol][igrow] = work[j * naroc[0] + i]; + } } } return 0; diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp index 1d0f480a32..7c7739dac9 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp @@ -400,8 +400,10 @@ void Wavefunc_in_pw::produce_local_basis_in_pw(const int& ik, for(int m = 0;m<2*L+1;m++) { const int lm = L*L +m; - if (iwall + 2 * L + 1 > GlobalC::ucell.natomwfc) - ModuleBase::WARNING_QUIT("this->wf.atomic_wfc()", "error: too many wfcs"); + if (iwall + 2 * L + 1 > GlobalC::ucell.natomwfc) + { + ModuleBase::WARNING_QUIT("this->wf.atomic_wfc()", "error: too many wfcs"); + } for (int ig = 0; ig < npw; ig++) { aux[ig] = sk[ig] * ylm(lm,ig) * flq[ig]; @@ -430,7 +432,8 @@ void Wavefunc_in_pw::produce_local_basis_in_pw(const int& ik, } // end else GlobalC::ucell.atoms[it].has_so } // end for is_N } // end if GlobalV::NONCOLIN - else{//LSDA and nomagnet case + else + {//LSDA and nomagnet case for(int m=0; m<2*L+1; m++) { const int lm = L*L+m; @@ -455,241 +458,3 @@ void Wavefunc_in_pw::produce_local_basis_in_pw(const int& ik, delete[] gk; } -// void Wavefunc_in_pw::produce_local_basis_q_in_pw(const int &ik, -// ModuleBase::ComplexMatrix &psi, -// ModulePW::PW_Basis_K *wfc_basis, -// const ModuleBase::realArray &table_local, -// ModuleBase::Vector3 q) // pengfei 2016-11-23 -// { -// ModuleBase::TITLE("Wavefunc_in_pw","produce_local_basis_in_pw"); -// assert(ik>=0); -// const int npw = kv.ngk[ik]; -// const int total_lm = ( GlobalC::ucell.lmax + 1) * ( GlobalC::ucell.lmax + 1); -// ModuleBase::matrix ylm(total_lm, npw); -// std::complex *aux = new std::complex[npw]; -// double *chiaux = nullptr; - -// ModuleBase::Vector3 *gkq = new ModuleBase::Vector3[npw]; - -// for(int ig=0;iggetgpluskcar(ik, ig) + q; -// } - -// ModuleBase::YlmReal::Ylm_Real(total_lm, npw, gkq, ylm); - -// //int index = 0; -// double *flq = new double[npw]; -// int iwall=0; -// for (int it = 0;it < GlobalC::ucell.ntype;it++) -// { -// for (int ia = 0;ia < GlobalC::ucell.atoms[it].na;ia++) -// { -// std::complex *skq = GlobalC::sf.get_skq(ik, it, ia, wfc_basis, q); -// int ic=0; -// for(int L = 0; L < GlobalC::ucell.atoms[it].nwl+1; L++) -// { -// std::complex lphase = pow(ModuleBase::NEG_IMAG_UNIT, L); //mohan 2010-04-19 -// for(int N=0; N < GlobalC::ucell.atoms[it].l_nchi[L]; N++) -// { -// for(int ig=0; ig ((GlobalV::NQX-4) * GlobalV::DQ) ) -// { -// flq[ig] = 0.0; -// } -// else -// { -// flq[ig] = ModuleBase::PolyInt::Polynomial_Interpolation(table_local, it, ic, GlobalV::NQX, -// GlobalV::DQ, gkq[ig].norm() * GlobalC::ucell.tpiba ); -// } -// } - -// if(GlobalV::NSPIN==4) -// { -// Soc soc; -// soc.rot_ylm(GlobalC::ucell.atoms[it].nwl+1); -// for(int is_N = 0; is_N < 2; is_N++) -// { -// if(L==0&&is_N==1) continue; -// if(GlobalC::ucell.atoms[it].ncpp.has_so) -// { -// const double j = double(L+is_N) - 0.5; -// if ( !(GlobalV::DOMAG||GlobalV::DOMAG_Z)) -// {//atomic_wfc_so -// double fact[2]; -// for(int m=-L-1;m1e-8||fabs(fact[1])>1e-8) -// { -// for(int is=0;is<2;is++) -// { -// if(fabs(fact[is])>1e-8) -// { -// const int ind = GlobalC::ppcell.lmaxkb + soc.sph_ind(L,j,m,is); -// ModuleBase::GlobalFunc::ZEROS(aux, npw); -// for(int n1=0;n1<2*L+1;n1++){ -// const int lm = L*L +n1; -// if(std::abs(soc.rotylm(n1,ind))>1e-8) -// for(int ig=0; ignpwk_max*is ) = lphase * fact[is] * skq[ig] -// * aux[ig] -// * flq[ig]; -// } -// else -// for(int ig=0; ignpwk_max*is) = -// std::complex(0.0 , 0.0); -// }//is -// iwall++; -// }//if -// }//m -// }//if -// else -// {//atomic_wfc_so_mag - -// double alpha, gamma; -// std::complex fup,fdown; -// //int nc; -// //This routine creates two functions only in the case j=l+1/2 or exit in the other -// case if(fabs(j-L+0.5<1e-4)) continue; delete[] chiaux; chiaux = new double [npw]; -// //Find the functions j= l- 1/2 -// if(L==0) -// for(int ig=0;igGlobalC::ucell.natomwfc) -// ModuleBase::WARNING_QUIT("this->wf.atomic_wfc()","error: too many wfcs"); for(int ig = 0;ignpwk_max) = (cos(0.5 * gamma) - ModuleBase::IMAG_UNIT * -// sin(0.5*gamma)) -// * fdown; -// //second rotation with angle gamma around(OZ) -// fup = cos(0.5 * (alpha + ModuleBase::PI))*aux[ig]; -// fdown = ModuleBase::IMAG_UNIT * sin(0.5 * (alpha + ModuleBase::PI))*aux[ig]; -// psi(iwall+2*L+1,ig) = (cos(0.5*gamma) + -// ModuleBase::IMAG_UNIT*sin(0.5*gamma))*fup; psi(iwall+2*L+1,ig+ wfc_basis->npwk_max) -// = (cos(0.5*gamma) -// - ModuleBase::IMAG_UNIT*sin(0.5*gamma))*fdown; -// } -// iwall++; -// } -// iwall += 2*L +1; -// } -// } -// else -// {//atomic_wfc_nc -// double alpha, gamman; -// std::complex fup, fdown; -// //alpha = GlobalC::ucell.magnet.angle1_[it]; -// //gamman = -GlobalC::ucell.magnet.angle2_[it] + 0.5*ModuleBase::PI; -// alpha = GlobalC::ucell.atoms[it].angle1[ia]; -// gamman = -GlobalC::ucell.atoms[it].angle2[ia]+ 0.5*ModuleBase::PI; -// for(int m = 0;m<2*L+1;m++) -// { -// const int lm = L*L +m; -// // if(iwall+2*l+1>GlobalC::ucell.natomwfc) -// ModuleBase::WARNING_QUIT("this->wf.atomic_wfc()","error: too many wfcs"); for(int -// ig = 0;ignpwk_max) = (cos(0.5 * gamman) - ModuleBase::IMAG_UNIT * -// sin(0.5*gamman)) -// * fdown; -// //second rotation with angle gamma around(OZ) -// fup = cos(0.5 * (alpha + ModuleBase::PI)) * aux[ig]; -// fdown = ModuleBase::IMAG_UNIT * sin(0.5 * (alpha + ModuleBase::PI)) * aux[ig]; -// psi(iwall+2*L+1,ig) = (cos(0.5*gamman) + -// ModuleBase::IMAG_UNIT*sin(0.5*gamman))*fup; psi(iwall+2*L+1,ig+ wfc_basis->npwk_max) -// = (cos(0.5*gamman) -// - ModuleBase::IMAG_UNIT*sin(0.5*gamman))*fdown; -// } -// iwall++; -// } -// iwall += 2*L+1; -// } -// // iwall++; -// }//end is_N -// }//end if -// else{//LSDA and nomagnet case -// for(int m=0; m<2*L+1; m++) -// { -// const int lm = L*L+m; -// for(int ig=0; ig q); // pengfei 2016-11-23 } #endif diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces.cpp index e1a81d9217..7c2c15321f 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/forces.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/forces.cpp @@ -420,16 +420,22 @@ void Forces::cal_force(ModuleBase::matrix& force, { ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "PAW FORCE (eV/Angstrom)", forcepaw, 0); } - if (GlobalV::EFIELD_FLAG) - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "EFIELD FORCE (eV/Angstrom)", force_e, 0); - if (GlobalV::GATE_FLAG) - ModuleIO::print_force(GlobalV::ofs_running, - GlobalC::ucell, - "GATEFIELD FORCE (eV/Angstrom)", - force_gate, - 0); - if (GlobalV::imp_sol) - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "IMP_SOL FORCE (eV/Angstrom)", forcesol, 0); + if (GlobalV::EFIELD_FLAG) + { + ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "EFIELD FORCE (eV/Angstrom)", force_e, 0); + } + if (GlobalV::GATE_FLAG) + { + ModuleIO::print_force(GlobalV::ofs_running, + GlobalC::ucell, + "GATEFIELD FORCE (eV/Angstrom)", + force_gate, + 0); + } + if (GlobalV::imp_sol) + { + ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "IMP_SOL FORCE (eV/Angstrom)", forcesol, 0); + } } ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "TOTAL-FORCE (eV/Angstrom)", force, 0); @@ -1304,4 +1310,4 @@ void Forces::cal_force_scc(ModuleBase::matrix& forcescc, template class Forces; #if ((defined __CUDA) || (defined __ROCM)) template class Forces; -#endif \ No newline at end of file +#endif diff --git a/source/module_io/output_dm1.cpp b/source/module_io/output_dm1.cpp index 9087431bf2..21f78295e9 100644 --- a/source/module_io/output_dm1.cpp +++ b/source/module_io/output_dm1.cpp @@ -5,32 +5,39 @@ namespace ModuleIO { -Output_DM1::Output_DM1(int nspin, int istep, Local_Orbital_Charge& LOC, Record_adj& RA, - K_Vectors& kv, const elecstate::DensityMatrix,double>* DM) - : _nspin(nspin), _istep(istep), _LOC(LOC), _RA(RA), _kv(kv), _DM(DM) +Output_DM1::Output_DM1( + int nspin, + int istep, + Local_Orbital_Charge& LOC, + Record_adj &RA, + K_Vectors& kv, + const elecstate::DensityMatrix,double>* DM) + : _nspin(nspin), _istep(istep), loc(LOC), _RA(RA), _kv(kv), _DM(DM) { } -void Output_DM1::write() +void Output_DM1::write(void) { double** dm2d; dm2d = new double*[_nspin]; for (int is = 0; is < _nspin; is++) { - dm2d[is] = new double[_LOC.ParaV->nnr]; - ModuleBase::GlobalFunc::ZEROS(dm2d[is], _LOC.ParaV->nnr); + dm2d[is] = new double[loc.ParaV->nnr]; + ModuleBase::GlobalFunc::ZEROS(dm2d[is], loc.ParaV->nnr); } + // get DMK from DM - _LOC.dm_k.resize(_kv.nks); + loc.dm_k.resize(_kv.nks); for (int ik = 0; ik < _kv.nks; ++ik) { - _LOC.set_dm_k(ik,_DM->get_DMK_pointer(ik)); + loc.set_dm_k(ik,_DM->get_DMK_pointer(ik)); } + // cal DMR in LOC - _LOC.cal_dm_R(_LOC.dm_k, _RA, dm2d, _kv); + loc.cal_dm_R(loc.dm_k, _RA, dm2d, _kv); for (int is = 0; is < _nspin; is++) { - write_dm1(is, _istep, dm2d, _LOC.ParaV, _LOC.DMR_sparse); + write_dm1(is, _istep, dm2d, loc.ParaV, loc.DMR_sparse); } for (int is = 0; is < _nspin; is++) @@ -40,4 +47,4 @@ void Output_DM1::write() delete[] dm2d; } -} // namespace ModuleIO \ No newline at end of file +} // namespace ModuleIO diff --git a/source/module_io/output_dm1.h b/source/module_io/output_dm1.h index fd0e39723f..1e835077d8 100644 --- a/source/module_io/output_dm1.h +++ b/source/module_io/output_dm1.h @@ -14,18 +14,27 @@ namespace ModuleIO class Output_DM1 : public Output_Interface { public: - Output_DM1(int nspin, int istep, Local_Orbital_Charge& LOC, Record_adj& RA, K_Vectors& kv, const elecstate::DensityMatrix,double>* DM); - void write() override; + Output_DM1( + int nspin, + int istep, + Local_Orbital_Charge& LOC, + Record_adj& RA, + K_Vectors& kv, + const elecstate::DensityMatrix,double>* DM); + + void write() override; private: - int _nspin; - int _istep; - Local_Orbital_Charge& _LOC; - Record_adj& _RA; - K_Vectors& _kv; - const elecstate::DensityMatrix,double>* _DM; + + int _nspin; + int _istep; + Local_Orbital_Charge& loc; + Record_adj& _RA; + K_Vectors& _kv; + const elecstate::DensityMatrix,double>* _DM; + }; } // namespace ModuleIO -#endif \ No newline at end of file +#endif diff --git a/source/module_io/write_dm_sparse.cpp b/source/module_io/write_dm_sparse.cpp index f6851a9cf8..a67f08167d 100644 --- a/source/module_io/write_dm_sparse.cpp +++ b/source/module_io/write_dm_sparse.cpp @@ -6,8 +6,13 @@ #include "module_cell/module_neighbor/sltk_grid_driver.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" -void ModuleIO::write_dm1(const int &is, const int &istep, double** dm2d, const Parallel_Orbitals* ParaV, - std::map, std::map>> &DMR_sparse) +void ModuleIO::write_dm1( + const int &is, + const int &istep, + double** dm2d, + const Parallel_Orbitals* ParaV, + std::map, + std::map>> &DMR_sparse) { ModuleBase::timer::tick("ModuleIO","write_dm"); ModuleBase::TITLE("ModuleIO","write_dm"); @@ -19,8 +24,11 @@ void ModuleIO::write_dm1(const int &is, const int &istep, double** dm2d, const P ModuleBase::timer::tick("ModuleIO","write_dm"); } -void ModuleIO::get_dm_sparse(const int &is, double** dm2d, const Parallel_Orbitals* ParaV, - std::map, std::map>> &DMR_sparse) +void ModuleIO::get_dm_sparse( + const int &is, + double** dm2d, + const Parallel_Orbitals* ParaV, + std::map, std::map>> &DMR_sparse) { ModuleBase::timer::tick("ModuleIO","get_dm_sparse"); ModuleBase::TITLE("ModuleIO","get_dm_sparse"); @@ -31,7 +39,7 @@ void ModuleIO::get_dm_sparse(const int &is, double** dm2d, const Parallel_Orbita ModuleBase::Vector3 dtau, tau1, tau2; ModuleBase::Vector3 dtau1, dtau2, tau0; - double temp_value_double; + double temp_value_double = 0.0; for(int T1 = 0; T1 < GlobalC::ucell.ntype; ++T1) { @@ -56,7 +64,10 @@ void ModuleIO::get_dm_sparse(const int &is, double** dm2d, const Parallel_Orbita bool adj = false; - if(distance < rcut) adj = true; + if(distance < rcut) + { + adj = true; + } else if(distance >= rcut) { for(int ad0 = 0; ad0 < GlobalC::GridD.getAdjacentNum()+1; ++ad0) @@ -85,21 +96,24 @@ void ModuleIO::get_dm_sparse(const int &is, double** dm2d, const Parallel_Orbita { const int start2 = GlobalC::ucell.itiaiw2iwt(T2,I2,0); - Abfs::Vector3_Order dR(GlobalC::GridD.getBox(ad).x, GlobalC::GridD.getBox(ad).y, GlobalC::GridD.getBox(ad).z); + Abfs::Vector3_Order dR( + GlobalC::GridD.getBox(ad).x, + GlobalC::GridD.getBox(ad).y, + GlobalC::GridD.getBox(ad).z); for(int ii=0; iinw; ii++) { const int iw1_all = start + ii; const int mu = ParaV->global2local_row(iw1_all); - if(mu<0)continue; + if(mu<0) continue; for(int jj=0; jjnw; jj++) { int iw2_all = start2 + jj; const int nu = ParaV->global2local_col(iw2_all); - if(nu<0)continue; + if(nu<0) continue; temp_value_double = dm2d[is][index]; if (std::abs(temp_value_double) > sparse_threshold) @@ -118,8 +132,12 @@ void ModuleIO::get_dm_sparse(const int &is, double** dm2d, const Parallel_Orbita ModuleBase::timer::tick("ModuleIO","get_dm_sparse"); } -void ModuleIO::write_dm_sparse(const int &is, const int &istep, const Parallel_Orbitals* ParaV, - std::map, std::map>> &DMR_sparse) + +void ModuleIO::write_dm_sparse( + const int &is, + const int &istep, + const Parallel_Orbitals* ParaV, + std::map, std::map>> &DMR_sparse) { ModuleBase::timer::tick("ModuleIO","write_dm_sparse"); ModuleBase::TITLE("ModuleIO","write_dm_sparse"); @@ -245,6 +263,7 @@ void ModuleIO::write_dm_sparse(const int &is, const int &istep, const Parallel_O ModuleBase::timer::tick("ModuleIO","write_dm_sparse"); } + void ModuleIO::destroy_dm_sparse(std::map, std::map>> &DMR_sparse) { std::map, std::map>> empty_DMR_sparse; diff --git a/source/module_io/write_dm_sparse.h b/source/module_io/write_dm_sparse.h index fae6d43ad9..0a5ad3e498 100644 --- a/source/module_io/write_dm_sparse.h +++ b/source/module_io/write_dm_sparse.h @@ -10,10 +10,13 @@ namespace ModuleIO { void write_dm1(const int &is, const int &istep, double** dm2d, const Parallel_Orbitals* ParaV, std::map, std::map>> &DMR_sparse); + void get_dm_sparse(const int &is, double** dm2d, const Parallel_Orbitals* ParaV, std::map, std::map>> &DMR_sparse); + void write_dm_sparse(const int &is, const int &istep, const Parallel_Orbitals* ParaV, std::map, std::map>> &DMR_sparse); + void destroy_dm_sparse( std::map, std::map>> &DMR_sparse); } diff --git a/source/module_md/md_func.cpp b/source/module_md/md_func.cpp index b5e1d80b74..05432471ce 100644 --- a/source/module_md/md_func.cpp +++ b/source/module_md/md_func.cpp @@ -6,11 +6,11 @@ namespace MD_func { -double gaussrand() +double gaussrand(void) { static double V1, V2, S; static int phase = 0; - double X; + double X=0.0; if (phase == 0) { @@ -228,16 +228,16 @@ void force_virial(ModuleESolver::ESolver* p_esolver, ModuleBase::TITLE("MD_func", "force_virial"); ModuleBase::timer::tick("MD_func", "force_virial"); - p_esolver->Run(istep, unit_in); + p_esolver->run(istep, unit_in); - potential = p_esolver->cal_Energy(); + potential = p_esolver->cal_energy(); ModuleBase::matrix force_temp(unit_in.nat, 3); - p_esolver->cal_Force(force_temp); + p_esolver->cal_force(force_temp); if (cal_stress) { - p_esolver->cal_Stress(virial); + p_esolver->cal_stress(virial); } /// convert Rydberg to Hartree @@ -472,4 +472,4 @@ void current_md_info(const int& my_rank, const std::string& file_dir, int& md_st #endif } -} // namespace MD_func \ No newline at end of file +} // namespace MD_func diff --git a/source/module_md/test/fire_test.cpp b/source/module_md/test/fire_test.cpp index c275ec2093..96a8e5612e 100644 --- a/source/module_md/test/fire_test.cpp +++ b/source/module_md/test/fire_test.cpp @@ -46,7 +46,7 @@ class FIREtest : public testing::Test Setcell::parameters(); ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); - p_esolver->Init(INPUT, ucell); + p_esolver->init(INPUT, ucell); mdrun = new FIRE(INPUT.mdp, ucell); mdrun->setup(p_esolver, GlobalV::global_readin_dir); @@ -208,4 +208,4 @@ TEST_F(FIREtest, PrintMD) EXPECT_THAT(output_str, testing::HasSubstr(" LARGEST GRAD (eV/A) : 0.049479926")); ifs.close(); remove("running.log"); -} \ No newline at end of file +} diff --git a/source/module_md/test/langevin_test.cpp b/source/module_md/test/langevin_test.cpp index 77bac815ce..8fc152c651 100644 --- a/source/module_md/test/langevin_test.cpp +++ b/source/module_md/test/langevin_test.cpp @@ -46,7 +46,7 @@ class Langevin_test : public testing::Test Setcell::parameters(); ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); - p_esolver->Init(INPUT, ucell); + p_esolver->init(INPUT, ucell); mdrun = new Langevin(INPUT.mdp, ucell); mdrun->setup(p_esolver, GlobalV::global_readin_dir); @@ -189,4 +189,4 @@ TEST_F(Langevin_test, print_md) " ------------------------------------------------------------------------------------------------")); ifs.close(); remove("running.log"); -} \ No newline at end of file +} diff --git a/source/module_md/test/lj_pot_test.cpp b/source/module_md/test/lj_pot_test.cpp index a6cbced392..89a062a1b0 100644 --- a/source/module_md/test/lj_pot_test.cpp +++ b/source/module_md/test/lj_pot_test.cpp @@ -35,7 +35,7 @@ class LJ_pot_test : public testing::Test Setcell::parameters(); ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); - p_esolver->Init(INPUT, ucell); + p_esolver->init(INPUT, ucell); MD_func::force_virial(p_esolver, 0, ucell, potential, force, true, stress); } @@ -77,4 +77,4 @@ TEST_F(LJ_pot_test, stress) EXPECT_NEAR(stress(2, 0), 0, doublethreshold); EXPECT_NEAR(stress(2, 1), -1.1858461261560206e-22, doublethreshold); EXPECT_NEAR(stress(2, 2), 6.4275429572682057e-07, doublethreshold); -} \ No newline at end of file +} diff --git a/source/module_md/test/msst_test.cpp b/source/module_md/test/msst_test.cpp index dcf3364873..a5d3ac97e4 100644 --- a/source/module_md/test/msst_test.cpp +++ b/source/module_md/test/msst_test.cpp @@ -46,7 +46,7 @@ class MSST_test : public testing::Test Setcell::parameters(); ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); - p_esolver->Init(INPUT, ucell); + p_esolver->init(INPUT, ucell); mdrun = new MSST(INPUT.mdp, ucell); mdrun->setup(p_esolver, GlobalV::global_readin_dir); diff --git a/source/module_md/test/nhchain_test.cpp b/source/module_md/test/nhchain_test.cpp index 93ea37b8eb..49ec231081 100644 --- a/source/module_md/test/nhchain_test.cpp +++ b/source/module_md/test/nhchain_test.cpp @@ -46,7 +46,7 @@ class NHC_test : public testing::Test Setcell::parameters(); ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); - p_esolver->Init(INPUT, ucell); + p_esolver->init(INPUT, ucell); INPUT.mdp.md_type = "npt"; INPUT.mdp.md_pmode = "tri"; @@ -237,4 +237,4 @@ TEST_F(NHC_test, print_md) " ------------------------------------------------------------------------------------------------")); ifs.close(); remove("running.log"); -} \ No newline at end of file +} diff --git a/source/module_md/test/verlet_test.cpp b/source/module_md/test/verlet_test.cpp index 9b2f138d20..8852743edc 100644 --- a/source/module_md/test/verlet_test.cpp +++ b/source/module_md/test/verlet_test.cpp @@ -46,7 +46,7 @@ class Verlet_test : public testing::Test Setcell::parameters(); ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); - p_esolver->Init(INPUT, ucell); + p_esolver->init(INPUT, ucell); mdrun = new Verlet(INPUT.mdp, ucell); mdrun->setup(p_esolver, GlobalV::global_readin_dir); @@ -330,4 +330,4 @@ TEST_F(Verlet_test, print_md) " ------------------------------------------------------------------------------------------------")); ifs.close(); remove("running.log"); -} \ No newline at end of file +} diff --git a/source/module_relax/relax_driver.cpp b/source/module_relax/relax_driver.cpp index e9d8a092a2..28f2485bc9 100644 --- a/source/module_relax/relax_driver.cpp +++ b/source/module_relax/relax_driver.cpp @@ -48,7 +48,7 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver *p_esolve #endif //__RAPIDJSON // mohan added eiter to count for the electron iteration number, 2021-01-28 - p_esolver->Run(istep - 1, GlobalC::ucell); + p_esolver->run(istep - 1, GlobalC::ucell); time_t eend = time(NULL); time_t fstart = time(NULL); @@ -62,17 +62,17 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver *p_esolve // but I'll use force and stress explicitly here for now // calculate the total energy - this->etot = p_esolver->cal_Energy(); + this->etot = p_esolver->cal_energy(); // calculate and gather all parts of total ionic forces if (GlobalV::CAL_FORCE) { - p_esolver->cal_Force(force); + p_esolver->cal_force(force); } // calculate and gather all parts of stress if (GlobalV::CAL_STRESS) { - p_esolver->cal_Stress(stress); + p_esolver->cal_stress(stress); } if (GlobalV::CALCULATION == "relax" || GlobalV::CALCULATION == "cell-relax") @@ -105,8 +105,12 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver *p_esolve GlobalC::ucell.print_cell_cif("STRU_NOW.cif"); } - ModuleESolver::ESolver_KS* p_esolver_ks = dynamic_cast*>(p_esolver); - if (p_esolver_ks && stop && p_esolver_ks->maxniter == p_esolver_ks->niter && !(p_esolver_ks->conv_elec)) + ModuleESolver::ESolver_KS* p_esolver_ks + = dynamic_cast*>(p_esolver); + if (p_esolver_ks + && stop + && p_esolver_ks->maxniter == p_esolver_ks->niter + && !(p_esolver_ks->conv_elec)) { std::cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; std::cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; @@ -153,4 +157,4 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver *p_esolve } template class Relax_Driver; -template class Relax_Driver; \ No newline at end of file +template class Relax_Driver;