diff --git a/source/driver.cpp b/source/driver.cpp index c4490da862..40b314b543 100644 --- a/source/driver.cpp +++ b/source/driver.cpp @@ -183,7 +183,7 @@ void Driver::atomic_world() //-------------------------------------------------- // where the actual stuff is done - this->driver_run(); + this->driver_run(GlobalC::ucell); ModuleBase::timer::finish(GlobalV::ofs_running); ModuleBase::Memory::print_all(GlobalV::ofs_running); diff --git a/source/driver.h b/source/driver.h index 00137821d5..fdcb83fade 100644 --- a/source/driver.h +++ b/source/driver.h @@ -1,6 +1,8 @@ #ifndef DRIVER_H #define DRIVER_H +#include "module_cell/unitcell.h" + class Driver { public: @@ -34,7 +36,7 @@ class Driver void atomic_world(); // the actual calculations - void driver_run(); + void driver_run(UnitCell& ucell); }; #endif diff --git a/source/driver_run.cpp b/source/driver_run.cpp index f4a055ca15..90d2c505f5 100644 --- a/source/driver_run.cpp +++ b/source/driver_run.cpp @@ -24,7 +24,8 @@ * the configuration-changing subroutine takes force and stress and updates the * configuration */ -void Driver::driver_run() { +void Driver::driver_run(UnitCell& ucell) +{ ModuleBase::TITLE("Driver", "driver_line"); ModuleBase::timer::tick("Driver", "driver_line"); @@ -39,20 +40,18 @@ void Driver::driver_run() { #endif // the life of ucell should begin here, mohan 2024-05-12 - // delete ucell as a GlobalC in near future - GlobalC::ucell.setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running); - Check_Atomic_Stru::check_atomic_stru(GlobalC::ucell, - PARAM.inp.min_dist_coef); + ucell.setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running); + Check_Atomic_Stru::check_atomic_stru(ucell, PARAM.inp.min_dist_coef); //! 2: initialize the ESolver (depends on a set-up ucell after `setup_cell`) - ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(PARAM.inp, GlobalC::ucell); + ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(PARAM.inp, ucell); //! 3: initialize Esolver and fill json-structure - p_esolver->before_all_runners(PARAM.inp, GlobalC::ucell); + p_esolver->before_all_runners(ucell, PARAM.inp); // this Json part should be moved to before_all_runners, mohan 2024-05-12 #ifdef __RAPIDJSON - Json::gen_stru_wrapper(&GlobalC::ucell); + Json::gen_stru_wrapper(&ucell); #endif const std::string cal_type = PARAM.inp.calculation; @@ -60,16 +59,16 @@ void Driver::driver_run() { //! 4: different types of calculations if (cal_type == "md") { - Run_MD::md_line(GlobalC::ucell, p_esolver, PARAM); + Run_MD::md_line(ucell, p_esolver, PARAM); } else if (cal_type == "scf" || cal_type == "relax" || cal_type == "cell-relax" || cal_type == "nscf") { Relax_Driver rl_driver; - rl_driver.relax_driver(p_esolver); + rl_driver.relax_driver(p_esolver, ucell); } else if (cal_type == "get_S") { - p_esolver->runner(0, GlobalC::ucell); + p_esolver->runner(ucell, 0); } else { @@ -79,11 +78,11 @@ void Driver::driver_run() { //! test_neighbour(LCAO), //! gen_bessel(PW), et al. const int istep = 0; - p_esolver->others(istep); + p_esolver->others(ucell, istep); } //! 5: clean up esolver - p_esolver->after_all_runners(); + p_esolver->after_all_runners(ucell); ModuleESolver::clean_esolver(p_esolver); diff --git a/source/module_base/opt_TN.hpp b/source/module_base/opt_TN.hpp index 5abcc53be8..aad30b0f60 100644 --- a/source/module_base/opt_TN.hpp +++ b/source/module_base/opt_TN.hpp @@ -1,9 +1,9 @@ #ifndef OPT_TN_H #define OPT_TN_H -#include +#include "opt_CG.h" -#include "./opt_CG.h" +#include namespace ModuleBase { @@ -25,7 +25,7 @@ class Opt_TN { this->mach_prec_ = std::numeric_limits::epsilon(); // get machine precise } - ~Opt_TN(){}; + ~Opt_TN() {}; /** * @brief Allocate the space for the arrays in cg_. @@ -54,7 +54,9 @@ class Opt_TN { this->iter_ = 0; if (nx_new != 0) + { this->nx_ = nx_new; + } this->cg_.refresh(nx_new); } @@ -167,17 +169,23 @@ void Opt_TN::next_direct(double* px, epsilon = this->get_epsilon(px, cg_direct); // epsilon = 1e-9; for (int i = 0; i < this->nx_; ++i) + { temp_x[i] = px[i] + epsilon * cg_direct[i]; + } (t->*p_calGradient)(temp_x, temp_gradient); for (int i = 0; i < this->nx_; ++i) + { temp_Hcgd[i] = (temp_gradient[i] - pgradient[i]) / epsilon; + } // get CG step length and update rdirect cg_alpha = cg_.step_length(temp_Hcgd, cg_direct, cg_ifPD); if (cg_ifPD == -1) // Hessian is not positive definite, and cgiter = 1. { for (int i = 0; i < this->nx_; ++i) + { rdirect[i] += cg_alpha * cg_direct[i]; + } flag = -1; break; } @@ -188,14 +196,18 @@ void Opt_TN::next_direct(double* px, } for (int i = 0; i < this->nx_; ++i) + { rdirect[i] += cg_alpha * cg_direct[i]; + } // store residuals used in truncated conditions last_residual = curr_residual; curr_residual = cg_.get_residual(); cg_iter = cg_.get_iter(); if (cg_iter == 1) + { init_residual = curr_residual; + } // check truncated conditions // if (curr_residual < 1e-12) diff --git a/source/module_cell/unitcell.cpp b/source/module_cell/unitcell.cpp index 6526b0fa75..653e632094 100755 --- a/source/module_cell/unitcell.cpp +++ b/source/module_cell/unitcell.cpp @@ -551,6 +551,9 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) { this->atoms = new Atom[this->ntype]; // atom species. this->set_atom_flag = true; + this->symm.epsilon = PARAM.inp.symmetry_prec; + this->symm.epsilon_input = PARAM.inp.symmetry_prec; + bool ok = true; bool ok2 = true; diff --git a/source/module_esolver/esolver.cpp b/source/module_esolver/esolver.cpp index e2a24148e8..86cb82c38b 100644 --- a/source/module_esolver/esolver.cpp +++ b/source/module_esolver/esolver.cpp @@ -247,8 +247,8 @@ ESolver* init_esolver(const Input_para& inp, UnitCell& ucell) { p_esolver = new ESolver_KS_LCAO, std::complex>(); } - p_esolver->before_all_runners(inp, ucell); - p_esolver->runner(0, ucell); // scf-only + p_esolver->before_all_runners(ucell, inp); + p_esolver->runner(ucell, 0); // scf-only // force and stress is not needed currently, // they will be supported after the analytical gradient // of LR-TDDFT is implemented. diff --git a/source/module_esolver/esolver.h b/source/module_esolver/esolver.h index b4be298ed6..800d40aeb9 100644 --- a/source/module_esolver/esolver.h +++ b/source/module_esolver/esolver.h @@ -20,26 +20,26 @@ class ESolver } //! initialize the energy solver by using input parameters and cell modules - virtual void before_all_runners(const Input_para& inp, UnitCell& cell) = 0; + virtual void before_all_runners(UnitCell& ucell, const Input_para& inp) = 0; //! run energy solver - virtual void runner(const int istep, UnitCell& cell) = 0; + virtual void runner(UnitCell& cell, const int istep) = 0; //! perform post processing calculations - virtual void after_all_runners(){}; + virtual void after_all_runners(UnitCell& ucell){}; //! deal with exx and other calculation than scf/md/relax/cell-relax: //! such as nscf, get_wf and get_pchg - virtual void others(const int istep){}; + virtual void others(UnitCell& ucell, const int istep) {}; //! calculate total energy of a given system virtual double cal_energy() = 0; //! calcualte forces for the atoms in the given cell - virtual void cal_force(ModuleBase::matrix& force) = 0; + virtual void cal_force(UnitCell& ucell, ModuleBase::matrix& force) = 0; //! calcualte stress of given cell - virtual void cal_stress(ModuleBase::matrix& stress) = 0; + virtual void cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) = 0; bool conv_esolver = true; // whether esolver is converged diff --git a/source/module_esolver/esolver_dp.cpp b/source/module_esolver/esolver_dp.cpp index ea2da70c98..9a1e4d52d0 100644 --- a/source/module_esolver/esolver_dp.cpp +++ b/source/module_esolver/esolver_dp.cpp @@ -31,7 +31,7 @@ namespace ModuleESolver { -void ESolver_DP::before_all_runners(const Input_para& inp, UnitCell& ucell) +void ESolver_DP::before_all_runners(UnitCell& ucell, const Input_para& inp) { ucell_ = &ucell; dp_potential = 0; @@ -57,7 +57,7 @@ void ESolver_DP::before_all_runners(const Input_para& inp, UnitCell& ucell) #endif } -void ESolver_DP::runner(const int istep, UnitCell& ucell) +void ESolver_DP::runner(UnitCell& ucell, const int istep) { ModuleBase::TITLE("ESolver_DP", "runner"); ModuleBase::timer::tick("ESolver_DP", "runner"); @@ -127,13 +127,13 @@ double ESolver_DP::cal_energy() return dp_potential; } -void ESolver_DP::cal_force(ModuleBase::matrix& force) +void ESolver_DP::cal_force(UnitCell& ucell, 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(UnitCell& ucell, ModuleBase::matrix& stress) { stress = dp_virial; @@ -148,7 +148,7 @@ void ESolver_DP::cal_stress(ModuleBase::matrix& stress) ModuleIO::print_stress("TOTAL-STRESS", stress, true, false); } -void ESolver_DP::after_all_runners() +void ESolver_DP::after_all_runners(UnitCell& ucell) { 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 96738cfea4..5b59e8f442 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 before_all_runners(const Input_para& inp, UnitCell& cell) override; + void before_all_runners(UnitCell& ucell, const Input_para& inp) 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 runner(const int istep, UnitCell& cell) override; + void runner(UnitCell& cell, const int istep) override; /** * @brief get the total energy without ion kinetic energy @@ -59,21 +59,21 @@ class ESolver_DP : public ESolver * * @param force the computed atomic forces */ - void cal_force(ModuleBase::matrix& force) override; + void cal_force(UnitCell& ucell, 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(UnitCell& ucell, 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 after_all_runners() override; + void after_all_runners(UnitCell& ucell) override; private: /** diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index e67fe546ad..a336f0b864 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -40,8 +40,6 @@ ESolver_FP::ESolver_FP() pw_big = static_cast(pw_rhod); pw_big->setbxyz(PARAM.inp.bx, PARAM.inp.by, PARAM.inp.bz); sf.set(pw_rhod, PARAM.inp.nbspline); - - GlobalC::ucell.symm.epsilon = GlobalC::ucell.symm.epsilon_input = PARAM.inp.symmetry_prec; } ESolver_FP::~ESolver_FP() @@ -54,14 +52,14 @@ ESolver_FP::~ESolver_FP() delete this->pelec; } -void ESolver_FP::before_all_runners(const Input_para& inp, UnitCell& cell) +void ESolver_FP::before_all_runners(UnitCell& ucell, const Input_para& inp) { ModuleBase::TITLE("ESolver_FP", "before_all_runners"); //! 1) read pseudopotentials if (!PARAM.inp.use_paw) { - cell.read_pseudo(GlobalV::ofs_running); + ucell.read_pseudo(GlobalV::ofs_running); } //! 2) initialie the plane wave basis for rho @@ -75,11 +73,11 @@ void ESolver_FP::before_all_runners(const Input_para& inp, UnitCell& cell) if (inp.nx * inp.ny * inp.nz == 0) { - this->pw_rho->initgrids(inp.ref_cell_factor * cell.lat0, cell.latvec, 4.0 * inp.ecutwfc); + this->pw_rho->initgrids(inp.ref_cell_factor * ucell.lat0, ucell.latvec, 4.0 * inp.ecutwfc); } else { - this->pw_rho->initgrids(inp.ref_cell_factor * cell.lat0, cell.latvec, inp.nx, inp.ny, inp.nz); + this->pw_rho->initgrids(inp.ref_cell_factor * ucell.lat0, ucell.latvec, inp.nx, inp.ny, inp.nz); } this->pw_rho->initparameters(false, 4.0 * inp.ecutwfc); @@ -101,11 +99,11 @@ void ESolver_FP::before_all_runners(const Input_para& inp, UnitCell& cell) } if (inp.ndx * inp.ndy * inp.ndz == 0) { - this->pw_rhod->initgrids(inp.ref_cell_factor * cell.lat0, cell.latvec, inp.ecutrho); + this->pw_rhod->initgrids(inp.ref_cell_factor * ucell.lat0, ucell.latvec, inp.ecutrho); } else { - this->pw_rhod->initgrids(inp.ref_cell_factor * cell.lat0, cell.latvec, inp.ndx, inp.ndy, inp.ndz); + this->pw_rhod->initgrids(inp.ref_cell_factor * ucell.lat0, ucell.latvec, inp.ndx, inp.ndy, inp.ndz); } this->pw_rhod->initparameters(false, inp.ecutrho); this->pw_rhod->fft_bundle.initfftmode(inp.fft_mode); @@ -113,8 +111,8 @@ void ESolver_FP::before_all_runners(const Input_para& inp, UnitCell& cell) this->pw_rhod->collect_local_pw(); this->pw_rhod->collect_uniqgg(); } - ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU.cif", - cell, + ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU.cif", + ucell, "# Generated by ABACUS ModuleIO::CifParser", "data_?"); @@ -122,13 +120,13 @@ void ESolver_FP::before_all_runners(const Input_para& inp, UnitCell& cell) ModuleIO::print_rhofft(this->pw_rhod, this->pw_rho, this->pw_big, GlobalV::ofs_running); //! 5) initialize the charge extrapolation method if necessary - this->CE.Init_CE(PARAM.inp.nspin, GlobalC::ucell.nat, this->pw_rhod->nrxx, inp.chg_extrap); + this->CE.Init_CE(PARAM.inp.nspin, ucell.nat, this->pw_rhod->nrxx, inp.chg_extrap); return; } //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. -void ESolver_FP::after_scf(const int istep) +void ESolver_FP::after_scf(UnitCell& ucell, const int istep) { // 0) output convergence information ModuleIO::output_convergence_after_scf(this->conv_esolver, this->pelec->f_en.etot); @@ -137,7 +135,7 @@ void ESolver_FP::after_scf(const int istep) ModuleIO::output_efermi(this->conv_esolver, this->pelec->eferm.ef); // 2) update delta rho for charge extrapolation - CE.update_delta_rho(GlobalC::ucell, &(this->chr), &(this->sf)); + CE.update_delta_rho(ucell, &(this->chr), &(this->sf)); if (istep % PARAM.inp.out_interval == 0) { @@ -159,26 +157,26 @@ void ESolver_FP::after_scf(const int istep) } std::string fn =PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_CHG.cube"; ModuleIO::write_vdata_palgrid(GlobalC::Pgrid, - data, - is, - PARAM.inp.nspin, - istep, - fn, - this->pelec->eferm.get_efval(is), - &(GlobalC::ucell), - PARAM.inp.out_chg[1], - 1); + data, + is, + PARAM.inp.nspin, + istep, + fn, + this->pelec->eferm.get_efval(is), + &(ucell), + PARAM.inp.out_chg[1], + 1); if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) { fn =PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_TAU.cube"; ModuleIO::write_vdata_palgrid(GlobalC::Pgrid, - this->pelec->charge->kin_r_save[is], - is, - PARAM.inp.nspin, - istep, - fn, - this->pelec->eferm.get_efval(is), - &(GlobalC::ucell)); + this->pelec->charge->kin_r_save[is], + is, + PARAM.inp.nspin, + istep, + fn, + this->pelec->eferm.get_efval(is), + &(ucell)); } } } @@ -194,7 +192,7 @@ void ESolver_FP::after_scf(const int istep) PARAM.globalv.gamma_only_pw || PARAM.globalv.gamma_only_local, this->pw_rhod, PARAM.inp.nspin, - GlobalC::ucell.GT, + ucell.GT, rhog_tot, GlobalV::MY_POOL, GlobalV::RANK_IN_POOL, @@ -209,15 +207,15 @@ void ESolver_FP::after_scf(const int istep) std::string fn =PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_POT.cube"; ModuleIO::write_vdata_palgrid(GlobalC::Pgrid, - this->pelec->pot->get_effective_v(is), - is, - PARAM.inp.nspin, - istep, - fn, - 0.0, // efermi - &(GlobalC::ucell), - 3, // precision - 0); // out_fermi + this->pelec->pot->get_effective_v(is), + is, + PARAM.inp.nspin, + istep, + fn, + 0.0, // efermi + &(ucell), + 3, // precision + 0); // out_fermi } } else if (PARAM.inp.out_pot == 2) @@ -232,7 +230,7 @@ void ESolver_FP::after_scf(const int istep) istep, this->pw_rhod, this->pelec->charge, - &(GlobalC::ucell), + &(ucell), this->pelec->pot->get_fixed_v()); } @@ -243,7 +241,7 @@ void ESolver_FP::after_scf(const int istep) Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) { - srho.begin(is, *(this->pelec->charge), this->pw_rhod, GlobalC::ucell.symm); + srho.begin(is, *(this->pelec->charge), this->pw_rhod, ucell.symm); } std::string out_dir =PARAM.globalv.global_out_dir; @@ -258,27 +256,27 @@ void ESolver_FP::after_scf(const int istep) this->pelec->charge->rho, this->pelec->charge->kin_r, this->pw_rhod, - &(GlobalC::ucell), + &(ucell), PARAM.inp.out_elf[1]); } } } -void ESolver_FP::before_scf(const int istep) +void ESolver_FP::before_scf(UnitCell& ucell, const int istep) { ModuleBase::TITLE("ESolver_FP", "before_scf"); - if (GlobalC::ucell.cell_parameter_updated) + if (ucell.cell_parameter_updated) { // only G-vector and K-vector are changed due to the change of lattice // vector FFT grids do not change!! - pw_rho->initgrids(GlobalC::ucell.lat0, GlobalC::ucell.latvec, pw_rho->nx, pw_rho->ny, pw_rho->nz); + pw_rho->initgrids(ucell.lat0, ucell.latvec, pw_rho->nx, pw_rho->ny, pw_rho->nz); pw_rho->collect_local_pw(); pw_rho->collect_uniqgg(); if (PARAM.globalv.double_grid) { - this->pw_rhod->initgrids(GlobalC::ucell.lat0, GlobalC::ucell.latvec, pw_rhod->nx, pw_rhod->ny, pw_rhod->nz); + this->pw_rhod->initgrids(ucell.lat0, ucell.latvec, pw_rhod->nx, pw_rhod->ny, pw_rhod->nz); this->pw_rhod->collect_local_pw(); this->pw_rhod->collect_uniqgg(); } @@ -286,18 +284,15 @@ void ESolver_FP::before_scf(const int istep) GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, pw_rhod); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); - this->pelec->omega = GlobalC::ucell.omega; + this->pelec->omega = ucell.omega; if (ModuleSymmetry::Symmetry::symm_flag == 1) { - GlobalC::ucell.symm.analy_sys(GlobalC::ucell.lat, - GlobalC::ucell.st, - GlobalC::ucell.atoms, - GlobalV::ofs_running); + ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY"); } - kv.set_after_vc(PARAM.inp.nspin, GlobalC::ucell.G, GlobalC::ucell.latvec); + kv.set_after_vc(PARAM.inp.nspin, ucell.G, ucell.latvec); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS"); } diff --git a/source/module_esolver/esolver_fp.h b/source/module_esolver/esolver_fp.h index e78cb3bd6d..708fda3c35 100644 --- a/source/module_esolver/esolver_fp.h +++ b/source/module_esolver/esolver_fp.h @@ -30,14 +30,14 @@ namespace ModuleESolver virtual ~ESolver_FP(); //! Initialize of the first-principels energy solver - virtual void before_all_runners(const Input_para& inp, UnitCell& cell) override; + virtual void before_all_runners(UnitCell& ucell, const Input_para& inp) override; protected: //! Something to do before SCF iterations. - virtual void before_scf(const int istep); + virtual void before_scf(UnitCell& ucell, const int istep); //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. - virtual void after_scf(const int istep); + virtual void after_scf(UnitCell& ucell, const int istep); //! Electronic states elecstate::ElecState* pelec = nullptr; diff --git a/source/module_esolver/esolver_gets.cpp b/source/module_esolver/esolver_gets.cpp index 097dbdd32d..60bfb983d3 100644 --- a/source/module_esolver/esolver_gets.cpp +++ b/source/module_esolver/esolver_gets.cpp @@ -25,7 +25,7 @@ ESolver_GetS::~ESolver_GetS() } template -void ESolver_GetS::before_all_runners(const Input_para& inp, UnitCell& ucell) +void ESolver_GetS::before_all_runners(UnitCell& ucell, const Input_para& inp) { ModuleBase::TITLE("ESolver_GetS", "before_all_runners"); ModuleBase::timer::tick("ESolver_GetS", "before_all_runners"); @@ -82,14 +82,14 @@ void ESolver_GetS::before_all_runners(const Input_para& inp, UnitCell& u } template <> -void ESolver_GetS::runner(const int istep, UnitCell& ucell) +void ESolver_GetS::runner(UnitCell& ucell, const int istep) { ModuleBase::TITLE("ESolver_GetS", "runner"); ModuleBase::WARNING_QUIT("ESolver_GetS::runner", "not implemented"); } template <> -void ESolver_GetS, std::complex>::runner(const int istep, UnitCell& ucell) +void ESolver_GetS, std::complex>::runner(UnitCell& ucell, const int istep) { ModuleBase::TITLE("ESolver_GetS", "runner"); ModuleBase::timer::tick("ESolver_GetS", "runner"); @@ -130,7 +130,7 @@ void ESolver_GetS, std::complex>::runner(const int } template <> -void ESolver_GetS, double>::runner(const int istep, UnitCell& ucell) +void ESolver_GetS, double>::runner(UnitCell& ucell, const int istep) { ModuleBase::TITLE("ESolver_GetS", "runner"); ModuleBase::timer::tick("ESolver_GetS", "runner"); diff --git a/source/module_esolver/esolver_gets.h b/source/module_esolver/esolver_gets.h index f675d66051..79cc012f29 100644 --- a/source/module_esolver/esolver_gets.h +++ b/source/module_esolver/esolver_gets.h @@ -18,20 +18,20 @@ class ESolver_GetS : public ESolver_KS ESolver_GetS(); ~ESolver_GetS(); - void before_all_runners(const Input_para& inp, UnitCell& ucell) override; + void before_all_runners(UnitCell& ucell, const Input_para& inp) override; - void after_all_runners() {}; + void after_all_runners(UnitCell& ucell){}; - void runner(const int istep, UnitCell& ucell) override; + void runner(UnitCell& ucell, const int istep) override; //! calculate total energy of a given system double cal_energy() {}; //! calcualte forces for the atoms in the given cell - void cal_force(ModuleBase::matrix& force) {}; + void cal_force(UnitCell& ucell, ModuleBase::matrix& force) {}; //! calcualte stress of given cell - void cal_stress(ModuleBase::matrix& stress) {}; + void cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) {}; protected: // we will get rid of this class soon, don't use it, mohan 2024-03-28 diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index 7981a5f433..a3f9d91710 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -87,12 +87,12 @@ ESolver_KS::~ESolver_KS() //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS::before_all_runners(const Input_para& inp, UnitCell& ucell) +void ESolver_KS::before_all_runners(UnitCell& ucell, const Input_para& inp) { ModuleBase::TITLE("ESolver_KS", "before_all_runners"); //! 1) initialize "before_all_runniers" in ESolver_FP - ESolver_FP::before_all_runners(inp, ucell); + ESolver_FP::before_all_runners(ucell, inp); //! 2) setup the charge mixing parameters p_chgmix->set_mixing(PARAM.inp.mixing_mode, @@ -332,7 +332,7 @@ void ESolver_KS::before_all_runners(const Input_para& inp, UnitCell& //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS::hamilt2density_single(const int istep, const int iter, const double ethr) +void ESolver_KS::hamilt2density_single(UnitCell& ucell, const int istep, const int iter, const double ethr) { ModuleBase::timer::tick(this->classname, "hamilt2density_single"); // Temporarily, before HSolver is constructed, it should be overrided by @@ -343,10 +343,10 @@ void ESolver_KS::hamilt2density_single(const int istep, const int ite } template -void ESolver_KS::hamilt2density(const int istep, const int iter, const double ethr) +void ESolver_KS::hamilt2density(UnitCell& ucell, const int istep, const int iter, const double ethr) { // 7) use Hamiltonian to obtain charge density - this->hamilt2density_single(istep, iter, diag_ethr); + this->hamilt2density_single(ucell, istep, iter, diag_ethr); // 8) for MPI: STOGROUP? need to rewrite // It may be changed when more clever parallel algorithm is @@ -383,7 +383,7 @@ void ESolver_KS::hamilt2density(const int istep, const int iter, cons diag_ethr, PARAM.inp.nelec); - this->hamilt2density_single(istep, iter, diag_ethr); + this->hamilt2density_single(ucell, istep, iter, diag_ethr); drho = p_chgmix->get_drho(pelec->charge, PARAM.inp.nelec); @@ -416,14 +416,14 @@ void ESolver_KS::hamilt2density(const int istep, const int iter, cons //! 16) Json again //------------------------------------------------------------------------------ template -void ESolver_KS::runner(const int istep, UnitCell& ucell) +void ESolver_KS::runner(UnitCell& ucell, const int istep) { ModuleBase::TITLE("ESolver_KS", "runner"); ModuleBase::timer::tick(this->classname, "runner"); // 2) before_scf (electronic iteration loops) ModuleBase::timer::tick(this->classname, "before_scf"); - this->before_scf(istep); + this->before_scf(ucell, istep); ModuleBase::timer::tick(this->classname, "before_scf"); // 3) write charge density @@ -442,12 +442,12 @@ void ESolver_KS::runner(const int istep, UnitCell& ucell) for (int iter = 1; iter <= this->maxniter; ++iter) { // 6) initialization of SCF iterations - this->iter_init(istep, iter); + this->iter_init(ucell, istep, iter); - this->hamilt2density(istep, iter, diag_ethr); + this->hamilt2density(ucell, istep, iter, diag_ethr); // 10) finish scf iterations - this->iter_finish(istep, iter); + this->iter_finish(ucell, istep, iter); // 13) check convergence if (this->conv_esolver || this->oscillate_esolver) @@ -463,7 +463,7 @@ void ESolver_KS::runner(const int istep, UnitCell& ucell) // 15) after scf ModuleBase::timer::tick(this->classname, "after_scf"); - this->after_scf(istep); + this->after_scf(ucell, istep); ModuleBase::timer::tick(this->classname, "after_scf"); ModuleBase::timer::tick(this->classname, "runner"); @@ -471,7 +471,7 @@ void ESolver_KS::runner(const int istep, UnitCell& ucell) }; template -void ESolver_KS::iter_init(const int istep, const int iter) +void ESolver_KS::iter_init(UnitCell& ucell, const int istep, const int iter) { ModuleIO::write_head(GlobalV::ofs_running, istep, iter, this->basisname); @@ -515,7 +515,7 @@ void ESolver_KS::iter_init(const int istep, const int iter) } template -void ESolver_KS::iter_finish(const int istep, int& iter) +void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& iter) { if (PARAM.inp.out_bandgap) { @@ -535,10 +535,10 @@ void ESolver_KS::iter_finish(const int istep, int& iter) } // compute magnetization, only for LSDA(spin==2) - GlobalC::ucell.magnet.compute_magnetization(this->pelec->charge->nrxx, - this->pelec->charge->nxyz, - this->pelec->charge->rho, - this->pelec->nelec_spin.data()); + ucell.magnet.compute_magnetization(this->pelec->charge->nrxx, + this->pelec->charge->nxyz, + this->pelec->charge->rho, + this->pelec->nelec_spin.data()); if (GlobalV::MY_STOGROUP == 0) { @@ -574,7 +574,7 @@ void ESolver_KS::iter_finish(const int istep, int& iter) if (this->scf_ene_thr > 0.0) { // calculate energy of output charge density - this->update_pot(istep, iter); + this->update_pot(ucell, istep, iter); this->pelec->cal_energies(2); // 2 means Kohn-Sham functional // now, etot_old is the energy of input density, while etot is the energy of output density this->pelec->f_en.etot_delta = this->pelec->f_en.etot - this->pelec->f_en.etot_old; @@ -632,7 +632,7 @@ void ESolver_KS::iter_finish(const int istep, int& iter) // update potential // Hamilt should be used after it is constructed. // this->phamilt->update(conv_esolver); - this->update_pot(istep, iter); + this->update_pot(ucell, istep, iter); // 1 means Harris-Foulkes functional // 2 means Kohn-Sham functional @@ -665,8 +665,8 @@ void ESolver_KS::iter_finish(const int istep, int& iter) // Json, need to be moved to somewhere else #ifdef __RAPIDJSON // add Json of scf mag - Json::add_output_scf_mag(GlobalC::ucell.magnet.tot_magnetization, - GlobalC::ucell.magnet.abs_magnetization, + Json::add_output_scf_mag(ucell.magnet.tot_magnetization, + ucell.magnet.abs_magnetization, this->pelec->f_en.etot * ModuleBase::Ry_to_eV, this->pelec->f_en.etot_delta * ModuleBase::Ry_to_eV, drho, @@ -683,10 +683,10 @@ void ESolver_KS::iter_finish(const int istep, int& iter) //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. template -void ESolver_KS::after_scf(const int istep) +void ESolver_KS::after_scf(UnitCell& ucell, const int istep) { // 1) call after_scf() of ESolver_FP - ESolver_FP::after_scf(istep); + ESolver_FP::after_scf(ucell, istep); // 2) write eigenvalues if (istep % PARAM.inp.out_interval == 0) diff --git a/source/module_esolver/esolver_ks.h b/source/module_esolver/esolver_ks.h index d47826739c..323b60316b 100644 --- a/source/module_esolver/esolver_ks.h +++ b/source/module_esolver/esolver_ks.h @@ -30,31 +30,31 @@ class ESolver_KS : public ESolver_FP //! Deconstructor virtual ~ESolver_KS(); - virtual void before_all_runners(const Input_para& inp, UnitCell& cell) override; + virtual void before_all_runners(UnitCell& ucell, const Input_para& inp) override; - virtual void runner(const int istep, UnitCell& cell) override; + virtual void runner(UnitCell& ucell, const int istep) override; protected: //! Something to do before SCF iterations. - virtual void before_scf(const int istep) {}; + virtual void before_scf(UnitCell& ucell, const int istep) {}; //! Something to do before hamilt2density function in each iter loop. - virtual void iter_init(const int istep, const int iter); + virtual void iter_init(UnitCell& ucell, const int istep, const int iter); //! Something to do after hamilt2density function in each iter loop. - virtual void iter_finish(const int istep, int& iter); + virtual void iter_finish(UnitCell& ucell, const int istep, int& iter); // calculate electron density from a specific Hamiltonian with ethr - virtual void hamilt2density_single(const int istep, const int iter, const double ethr); + virtual void hamilt2density_single(UnitCell& ucell, const int istep, const int iter, const double ethr); // calculate electron density from a specific Hamiltonian - void hamilt2density(const int istep, const int iter, const double ethr); + void hamilt2density(UnitCell& ucell, const int istep, const int iter, const double ethr); //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. - virtual void after_scf(const int istep) override; + virtual void after_scf(UnitCell& ucell, const int istep) override; //! It should be replaced by a function in Hamilt Class - virtual void update_pot(const int istep, const int iter) {}; + virtual void update_pot(UnitCell& ucell, const int istep, const int iter) {}; //! Hamiltonian hamilt::Hamilt* p_hamilt = nullptr; diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index e9901bedc6..ef884e1c5b 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -117,12 +117,12 @@ ESolver_KS_LCAO::~ESolver_KS_LCAO() //! 13) print a warning if needed //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::before_all_runners(const Input_para& inp, UnitCell& ucell) +void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_para& inp) { ModuleBase::TITLE("ESolver_KS_LCAO", "before_all_runners"); ModuleBase::timer::tick("ESolver_KS_LCAO", "before_all_runners"); - ESolver_KS::before_all_runners(inp, ucell); + ESolver_KS::before_all_runners(ucell, inp); // 2) init ElecState // autoset nbands in ElecState, it should before basis_init (for Psi 2d division) @@ -174,12 +174,12 @@ void ESolver_KS_LCAO::before_all_runners(const Input_para& inp, UnitCell if (GlobalC::exx_info.info_ri.real_number) { this->exx_lri_double->init(MPI_COMM_WORLD, this->kv, orb_); - this->exd->exx_before_all_runners(this->kv, GlobalC::ucell, this->pv); + this->exd->exx_before_all_runners(this->kv, ucell, this->pv); } else { this->exx_lri_complex->init(MPI_COMM_WORLD, this->kv, orb_); - this->exc->exx_before_all_runners(this->kv, GlobalC::ucell, this->pv); + this->exc->exx_before_all_runners(this->kv, ucell, this->pv); } } } @@ -196,14 +196,14 @@ void ESolver_KS_LCAO::before_all_runners(const Input_para& inp, UnitCell // 9) inititlize the charge density this->pelec->charge->allocate(PARAM.inp.nspin); - this->pelec->omega = GlobalC::ucell.omega; + this->pelec->omega = ucell.omega; // 10) initialize the potential if (this->pelec->pot == nullptr) { this->pelec->pot = new elecstate::Potential(this->pw_rhod, this->pw_rho, - &GlobalC::ucell, + &ucell, &(GlobalC::ppcell.vloc), &(this->sf), &(this->pelec->f_en.etxc), @@ -270,12 +270,12 @@ double ESolver_KS_LCAO::cal_energy() //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::cal_force(ModuleBase::matrix& force) +void ESolver_KS_LCAO::cal_force(UnitCell& ucell, ModuleBase::matrix& force) { ModuleBase::TITLE("ESolver_KS_LCAO", "cal_force"); ModuleBase::timer::tick("ESolver_KS_LCAO", "cal_force"); - Force_Stress_LCAO fsl(this->RA, GlobalC::ucell.nat); + Force_Stress_LCAO fsl(this->RA, ucell.nat); fsl.getForceStress(PARAM.inp.cal_force, PARAM.inp.cal_stress, @@ -297,7 +297,7 @@ void ESolver_KS_LCAO::cal_force(ModuleBase::matrix& force) *this->exx_lri_double, *this->exx_lri_complex, #endif - &GlobalC::ucell.symm); + &ucell.symm); // delete RA after cal_force @@ -313,7 +313,7 @@ void ESolver_KS_LCAO::cal_force(ModuleBase::matrix& force) //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::cal_stress(ModuleBase::matrix& stress) +void ESolver_KS_LCAO::cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) { ModuleBase::TITLE("ESolver_KS_LCAO", "cal_stress"); ModuleBase::timer::tick("ESolver_KS_LCAO", "cal_stress"); @@ -321,7 +321,7 @@ void ESolver_KS_LCAO::cal_stress(ModuleBase::matrix& stress) if (!this->have_force) { ModuleBase::matrix fcs; - this->cal_force(fcs); + this->cal_force(ucell, fcs); } stress = this->scs; // copy the stress this->have_force = false; @@ -334,7 +334,7 @@ void ESolver_KS_LCAO::cal_stress(ModuleBase::matrix& stress) //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::after_all_runners() +void ESolver_KS_LCAO::after_all_runners(UnitCell& ucell) { ModuleBase::TITLE("ESolver_KS_LCAO", "after_all_runners"); ModuleBase::timer::tick("ESolver_KS_LCAO", "after_all_runners"); @@ -404,7 +404,7 @@ void ESolver_KS_LCAO::after_all_runners() if (PARAM.inp.out_proj_band) // Projeced band structure added by jiyy-2022-4-20 { - ModuleIO::write_proj_band_lcao(this->psi, this->pv, this->pelec, this->kv, GlobalC::ucell, this->p_hamilt); + ModuleIO::write_proj_band_lcao(this->psi, this->pv, this->pelec, this->kv, ucell, this->p_hamilt); } if (PARAM.inp.out_dos) @@ -418,7 +418,7 @@ void ESolver_KS_LCAO::after_all_runners() PARAM.inp.dos_sigma, *(this->pelec->klist), GlobalC::Pkpoints, - GlobalC::ucell, + ucell, this->pelec->eferm, PARAM.inp.nbands, this->p_hamilt); @@ -431,7 +431,7 @@ void ESolver_KS_LCAO::after_all_runners() GlobalV::DRANK, &this->pv, *this->psi, - GlobalC::ucell, + ucell, this->sf, *this->pw_rho, *this->pw_rhod, @@ -458,7 +458,7 @@ void ESolver_KS_LCAO::after_all_runners() GlobalV::DRANK, &this->pv, *this->psi, - GlobalC::ucell, + ucell, this->sf, *this->pw_rho, *this->pw_rhod, @@ -487,12 +487,12 @@ void ESolver_KS_LCAO::after_all_runners() //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::iter_init(const int istep, const int iter) +void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const int iter) { ModuleBase::TITLE("ESolver_KS_LCAO", "iter_init"); // call iter_init() of ESolver_KS - ESolver_KS::iter_init(istep, iter); + ESolver_KS::iter_init(ucell, istep, iter); if (iter == 1) { @@ -593,11 +593,11 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) if (PARAM.inp.nspin == 4) { - GlobalC::ucell.cal_ux(); + ucell.cal_ux(); } //! update the potentials by using new electron charge density - this->pelec->pot->update_from_charge(this->pelec->charge, &GlobalC::ucell); + this->pelec->pot->update_from_charge(this->pelec->charge, &ucell); //! compute the correction energy for metals this->pelec->f_en.descf = this->pelec->cal_delta_escf(); @@ -682,7 +682,7 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) //! 12) calculate delta energy //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::hamilt2density_single(int istep, int iter, double ethr) +void ESolver_KS_LCAO::hamilt2density_single(UnitCell& ucell, int istep, int iter, double ethr) { ModuleBase::TITLE("ESolver_KS_LCAO", "hamilt2density_single"); @@ -735,7 +735,7 @@ void ESolver_KS_LCAO::hamilt2density_single(int istep, int iter, double Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) { - srho.begin(is, *(this->pelec->charge), this->pw_rho, GlobalC::ucell.symm); + srho.begin(is, *(this->pelec->charge), this->pw_rho, ucell.symm); } // 12) calculate delta energy @@ -750,7 +750,7 @@ void ESolver_KS_LCAO::hamilt2density_single(int istep, int iter, double //! 3) print potential //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::update_pot(const int istep, const int iter) +void ESolver_KS_LCAO::update_pot(UnitCell& ucell, const int istep, const int iter) { ModuleBase::TITLE("ESolver_KS_LCAO", "update_pot"); @@ -829,9 +829,9 @@ void ESolver_KS_LCAO::update_pot(const int istep, const int iter) { if (PARAM.inp.nspin == 4) { - GlobalC::ucell.cal_ux(); + ucell.cal_ux(); } - this->pelec->pot->update_from_charge(this->pelec->charge, &GlobalC::ucell); + this->pelec->pot->update_from_charge(this->pelec->charge, &ucell); this->pelec->f_en.descf = this->pelec->cal_delta_escf(); } else @@ -849,7 +849,7 @@ void ESolver_KS_LCAO::update_pot(const int istep, const int iter) //! 4) output charge density and density matrix //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::iter_finish(const int istep, int& iter) +void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& iter) { ModuleBase::TITLE("ESolver_KS_LCAO", "iter_finish"); @@ -891,7 +891,7 @@ void ESolver_KS_LCAO::iter_finish(const int istep, int& iter) } // call iter_finish() of ESolver_KS - ESolver_KS::iter_finish(istep, iter); + ESolver_KS::iter_finish(ucell, istep, iter); // 1) mix density matrix if mixing_restart + mixing_dmr + not first // mixing_restart at every iter @@ -918,7 +918,7 @@ void ESolver_KS_LCAO::iter_finish(const int istep, int& iter) if (GlobalC::exx_info.info_global.cal_exx) { GlobalC::exx_info.info_ri.real_number ? this->exd->exx_iter_finish(this->kv, - GlobalC::ucell, + ucell, *this->p_hamilt, *this->pelec, *this->p_chgmix, @@ -927,7 +927,7 @@ void ESolver_KS_LCAO::iter_finish(const int istep, int& iter) istep, this->conv_esolver) : this->exc->exx_iter_finish(this->kv, - GlobalC::ucell, + ucell, *this->p_hamilt, *this->pelec, *this->p_chgmix, @@ -955,26 +955,26 @@ void ESolver_KS_LCAO::iter_finish(const int istep, int& iter) } std::string fn = PARAM.globalv.global_out_dir + "/tmp_SPIN" + std::to_string(is + 1) + "_CHG.cube"; ModuleIO::write_vdata_palgrid(GlobalC::Pgrid, - data, - is, - PARAM.inp.nspin, - 0, - fn, - this->pelec->eferm.get_efval(is), - &(GlobalC::ucell), - 3, - 1); + data, + is, + PARAM.inp.nspin, + 0, + fn, + this->pelec->eferm.get_efval(is), + &(ucell), + 3, + 1); if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) { fn = PARAM.globalv.global_out_dir + "/tmp_SPIN" + std::to_string(is + 1) + "_TAU.cube"; ModuleIO::write_vdata_palgrid(GlobalC::Pgrid, - this->pelec->charge->kin_r_save[is], - is, - PARAM.inp.nspin, - 0, - fn, - this->pelec->eferm.get_efval(is), - &(GlobalC::ucell)); + this->pelec->charge->kin_r_save[is], + is, + PARAM.inp.nspin, + 0, + fn, + this->pelec->eferm.get_efval(is), + &(ucell)); } } } @@ -1003,7 +1003,7 @@ void ESolver_KS_LCAO::iter_finish(const int istep, int& iter) //! 18) write quasi-orbitals //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::after_scf(const int istep) +void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) { ModuleBase::TITLE("ESolver_KS_LCAO", "after_scf"); // 1) calculate the kinetic energy density tau, sunliang 2024-09-18 @@ -1013,7 +1013,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) } // 2) call after_scf() of ESolver_KS - ESolver_KS::after_scf(istep); + ESolver_KS::after_scf(ucell, istep); // 3) write density matrix for sparse matrix ModuleIO::write_dmr(dynamic_cast*>(this->pelec)->get_DM()->get_DMR_vector(), @@ -1035,7 +1035,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) ModuleIO::write_dmk(dynamic_cast*>(this->pelec)->get_DM()->get_DMK_vector(), precision, efermis, - &(GlobalC::ucell), + &(ucell), this->pv); } @@ -1049,11 +1049,11 @@ void ESolver_KS_LCAO::after_scf(const int istep) const std::string file_name_exx = PARAM.globalv.global_out_dir + "HexxR" + std::to_string(GlobalV::MY_RANK); if (GlobalC::exx_info.info_ri.real_number) { - ModuleIO::write_Hexxs_csr(file_name_exx, GlobalC::ucell, this->exd->get_Hexxs()); + ModuleIO::write_Hexxs_csr(file_name_exx, ucell, this->exd->get_Hexxs()); } else { - ModuleIO::write_Hexxs_csr(file_name_exx, GlobalC::ucell, this->exc->get_Hexxs()); + ModuleIO::write_Hexxs_csr(file_name_exx, ucell, this->exc->get_Hexxs()); } } } @@ -1066,11 +1066,11 @@ void ESolver_KS_LCAO::after_scf(const int istep) ModuleBase::timer::tick("ESolver_KS_LCAO", "out_deepks_labels"); LDI.out_deepks_labels(this->pelec->f_en.etot, this->pelec->klist->get_nks(), - GlobalC::ucell.nat, + ucell.nat, PARAM.globalv.nlocal, this->pelec->ekb, this->pelec->klist->kvec_d, - GlobalC::ucell, + ucell, orb_, GlobalC::GridD, &(this->pv), @@ -1104,7 +1104,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) hamilt::HamiltLCAO, double>* p_ham_lcao = dynamic_cast, double>*>(this->p_hamilt); std::string zipname = "output_HR0.npz"; - ModuleIO::output_mat_npz(GlobalC::ucell, zipname, *(p_ham_lcao->getHR())); + ModuleIO::output_mat_npz(ucell, zipname, *(p_ham_lcao->getHR())); if (PARAM.inp.nspin == 2) { @@ -1112,7 +1112,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) hamilt::HamiltLCAO, double>* p_ham_lcao = dynamic_cast, double>*>(this->p_hamilt); zipname = "output_HR1.npz"; - ModuleIO::output_mat_npz(GlobalC::ucell, zipname, *(p_ham_lcao->getHR())); + ModuleIO::output_mat_npz(ucell, zipname, *(p_ham_lcao->getHR())); } } @@ -1122,12 +1122,12 @@ void ESolver_KS_LCAO::after_scf(const int istep) const elecstate::DensityMatrix* dm = dynamic_cast*>(this->pelec)->get_DM(); std::string zipname = "output_DM0.npz"; - ModuleIO::output_mat_npz(GlobalC::ucell, zipname, *(dm->get_DMR_pointer(1))); + ModuleIO::output_mat_npz(ucell, zipname, *(dm->get_DMR_pointer(1))); if (PARAM.inp.nspin == 2) { zipname = "output_DM1.npz"; - ModuleIO::output_mat_npz(GlobalC::ucell, zipname, *(dm->get_DMR_pointer(2))); + ModuleIO::output_mat_npz(ucell, zipname, *(dm->get_DMR_pointer(2))); } } @@ -1144,14 +1144,14 @@ void ESolver_KS_LCAO::after_scf(const int istep) this->GK, // mohan add 2024-04-01 two_center_bundle_, orb_, - GlobalC::ucell, + ucell, GlobalC::GridD, // mohan add 2024-04-06 this->kv, this->p_hamilt); // mulliken charge analysis if (PARAM.inp.out_mul) { - ModuleIO::cal_mag(&(this->pv), this->p_hamilt, this->kv, this->pelec, GlobalC::ucell, istep, true); + ModuleIO::cal_mag(&(this->pv), this->p_hamilt, this->kv, this->pelec, ucell, istep, true); } } @@ -1178,7 +1178,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) tqo.initialize(PARAM.globalv.global_out_dir, PARAM.inp.pseudo_dir, PARAM.inp.orbital_dir, - &GlobalC::ucell, + &ucell, this->kv.kvec_d, GlobalV::ofs_running, GlobalV::MY_RANK, @@ -1194,7 +1194,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) = new hamilt::EkineticNew>(&hsk, this->kv.kvec_d, &hR, - &GlobalC::ucell, + &ucell, orb_.cutoffs(), &GlobalC::GridD, two_center_bundle_.kinetic_orb.get()); diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index 95cdd9f77f..20d11a68a7 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -27,30 +27,30 @@ class ESolver_KS_LCAO : public ESolver_KS { ESolver_KS_LCAO(); ~ESolver_KS_LCAO(); - void before_all_runners(const Input_para& inp, UnitCell& cell) override; + void before_all_runners(UnitCell& ucell, const Input_para& inp) override; double cal_energy() override; - void cal_force(ModuleBase::matrix& force) override; + void cal_force(UnitCell& ucell, ModuleBase::matrix& force) override; - void cal_stress(ModuleBase::matrix& stress) override; + void cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) override; - void after_all_runners() override; + void after_all_runners(UnitCell& ucell) override; protected: - virtual void before_scf(const int istep) override; + virtual void before_scf(UnitCell& ucell, const int istep) override; - virtual void iter_init(const int istep, const int iter) override; + virtual void iter_init(UnitCell& ucell, const int istep, const int iter) override; - virtual void hamilt2density_single(const int istep, const int iter, const double ethr) override; + virtual void hamilt2density_single(UnitCell& ucell, const int istep, const int iter, const double ethr) override; - virtual void update_pot(const int istep, const int iter) override; + virtual void update_pot(UnitCell& ucell, const int istep, const int iter) override; - virtual void iter_finish(const int istep, int& iter) override; + virtual void iter_finish(UnitCell& ucell, const int istep, int& iter) override; - virtual void after_scf(const int istep) override; + virtual void after_scf(UnitCell& ucell, const int istep) override; - virtual void others(const int istep) override; + virtual void others(UnitCell& ucell, const int istep) override; // we will get rid of this class soon, don't use it, mohan 2024-03-28 Record_adj RA; diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index 731016e603..92be005770 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -65,16 +65,16 @@ ESolver_KS_LCAO_TDDFT::~ESolver_KS_LCAO_TDDFT() } } -void ESolver_KS_LCAO_TDDFT::before_all_runners(const Input_para& inp, UnitCell& ucell) +void ESolver_KS_LCAO_TDDFT::before_all_runners(UnitCell& ucell, const Input_para& inp) { // 1) run before_all_runners in ESolver_KS_LCAO - ESolver_KS_LCAO, double>::before_all_runners(inp, ucell); + ESolver_KS_LCAO, double>::before_all_runners(ucell, inp); // this line should be optimized // this->pelec = dynamic_cast(this->pelec); } -void ESolver_KS_LCAO_TDDFT::hamilt2density_single(const int istep, const int iter, const double ethr) +void ESolver_KS_LCAO_TDDFT::hamilt2density_single(UnitCell& ucell, const int istep, const int iter, const double ethr) { if (PARAM.inp.init_wfc == "file") { @@ -133,7 +133,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density_single(const int istep, const int ite Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) { - srho.begin(is, *(pelec->charge), pw_rho, GlobalC::ucell.symm); + srho.begin(is, *(pelec->charge), pw_rho, ucell.symm); } } @@ -141,7 +141,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density_single(const int istep, const int ite this->pelec->f_en.deband = this->pelec->cal_delta_eband(); } -void ESolver_KS_LCAO_TDDFT::iter_finish(const int istep, int& iter) +void ESolver_KS_LCAO_TDDFT::iter_finish(UnitCell& ucell, const int istep, int& iter) { // print occupation of each band if (iter == 1 && istep <= 2) @@ -167,10 +167,10 @@ void ESolver_KS_LCAO_TDDFT::iter_finish(const int istep, int& iter) << std::endl; } - ESolver_KS_LCAO, double>::iter_finish(istep, iter); + ESolver_KS_LCAO, double>::iter_finish(ucell, istep, iter); } -void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) +void ESolver_KS_LCAO_TDDFT::update_pot(UnitCell& ucell, const int istep, const int iter) { // print Hamiltonian and Overlap matrix if (this->conv_esolver) @@ -239,9 +239,9 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) { if (PARAM.inp.nspin == 4) { - GlobalC::ucell.cal_ux(); + ucell.cal_ux(); } - this->pelec->pot->update_from_charge(this->pelec->charge, &GlobalC::ucell); + this->pelec->pot->update_from_charge(this->pelec->charge, &ucell); this->pelec->f_en.descf = this->pelec->cal_delta_escf(); } else @@ -343,7 +343,7 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) } } -void ESolver_KS_LCAO_TDDFT::after_scf(const int istep) +void ESolver_KS_LCAO_TDDFT::after_scf(UnitCell& ucell, const int istep) { for (int is = 0; is < PARAM.inp.nspin; is++) { @@ -368,7 +368,7 @@ void ESolver_KS_LCAO_TDDFT::after_scf(const int istep) orb_, this->RA); } - ESolver_KS_LCAO, double>::after_scf(istep); + ESolver_KS_LCAO, double>::after_scf(ucell, istep); } void ESolver_KS_LCAO_TDDFT::weight_dm_rho() diff --git a/source/module_esolver/esolver_ks_lcao_tddft.h b/source/module_esolver/esolver_ks_lcao_tddft.h index c36e41cb9d..40f4486475 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.h +++ b/source/module_esolver/esolver_ks_lcao_tddft.h @@ -15,16 +15,16 @@ class ESolver_KS_LCAO_TDDFT : public ESolver_KS_LCAO, doubl ~ESolver_KS_LCAO_TDDFT(); - void before_all_runners(const Input_para& inp, UnitCell& cell) override; + void before_all_runners(UnitCell& ucell, const Input_para& inp) override; protected: - virtual void hamilt2density_single(const int istep, const int iter, const double ethr) override; + virtual void hamilt2density_single(UnitCell& ucell, const int istep, const int iter, const double ethr) override; - virtual void update_pot(const int istep, const int iter) override; + virtual void update_pot(UnitCell& ucell, const int istep, const int iter) override; - virtual void iter_finish(const int istep, int& iter) override; + virtual void iter_finish(UnitCell& ucell, const int istep, int& iter) override; - virtual void after_scf(const int istep) override; + virtual void after_scf(UnitCell& ucell, const int istep) override; //! wave functions of last time step psi::Psi>* psi_laststep = nullptr; diff --git a/source/module_esolver/esolver_ks_lcaopw.cpp b/source/module_esolver/esolver_ks_lcaopw.cpp index 0eb22b4ebf..eb54bd9b66 100644 --- a/source/module_esolver/esolver_ks_lcaopw.cpp +++ b/source/module_esolver/esolver_ks_lcaopw.cpp @@ -81,18 +81,26 @@ namespace ModuleESolver } template - void ESolver_KS_LIP::before_all_runners(const Input_para& inp, UnitCell& cell) + void ESolver_KS_LIP::before_all_runners(UnitCell& ucell, const Input_para& inp) { - ESolver_KS_PW::before_all_runners(inp, cell); + ESolver_KS_PW::before_all_runners(ucell, inp); #ifdef __EXX if (PARAM.inp.calculation == "scf" || PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax" || PARAM.inp.calculation == "md") { if (GlobalC::exx_info.info_global.cal_exx) { - XC_Functional::set_xc_first_loop(cell); + XC_Functional::set_xc_first_loop(ucell); this->exx_lip = std::unique_ptr>(new Exx_Lip(GlobalC::exx_info.info_lip, - cell.symm, &this->kv, this->p_wf_init, this->kspw_psi, this->pw_wfc, this->pw_rho, this->sf, &cell, this->pelec)); + ucell.symm, + &this->kv, + this->p_wf_init, + this->kspw_psi, + this->pw_wfc, + this->pw_rho, + this->sf, + &ucell, + this->pelec)); // this->exx_lip.init(GlobalC::exx_info.info_lip, cell.symm, &this->kv, this->p_wf_init, this->kspw_psi, this->pw_wfc, this->pw_rho, this->sf, &cell, this->pelec); } } @@ -100,9 +108,9 @@ namespace ModuleESolver } template - void ESolver_KS_LIP::iter_init(const int istep, const int iter) + void ESolver_KS_LIP::iter_init(UnitCell& ucell, const int istep, const int iter) { - ESolver_KS_PW::iter_init(istep, iter); + ESolver_KS_PW::iter_init(ucell, istep, iter); #ifdef __EXX if (GlobalC::exx_info.info_global.cal_exx && !GlobalC::exx_info.info_global.separate_loop && this->two_level_step) { this->exx_lip->cal_exx(); @@ -111,7 +119,7 @@ namespace ModuleESolver } template - void ESolver_KS_LIP::hamilt2density_single(const int istep, const int iter, const double ethr) + void ESolver_KS_LIP::hamilt2density_single(UnitCell& ucell, const int istep, const int iter, const double ethr) { ModuleBase::TITLE("ESolver_KS_LIP", "hamilt2density_single"); ModuleBase::timer::tick("ESolver_KS_LIP", "hamilt2density_single"); @@ -152,7 +160,7 @@ namespace ModuleESolver Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) { - srho.begin(is, *(this->pelec->charge), this->pw_rhod, GlobalC::ucell.symm); + srho.begin(is, *(this->pelec->charge), this->pw_rhod, ucell.symm); } // deband is calculated from "output" charge density calculated @@ -164,9 +172,9 @@ namespace ModuleESolver } template - void ESolver_KS_LIP::iter_finish(const int istep, int& iter) + void ESolver_KS_LIP::iter_finish(UnitCell& ucell, const int istep, int& iter) { - ESolver_KS_PW::iter_finish(istep, iter); + ESolver_KS_PW::iter_finish(ucell, istep, iter); #ifdef __EXX if (GlobalC::exx_info.info_global.cal_exx && this->conv_esolver) @@ -183,7 +191,7 @@ namespace ModuleESolver if (!this->two_level_step) { // update exx and redo scf - XC_Functional::set_xc_type(GlobalC::ucell.atoms[0].ncpp.xc_func); + XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); iter = 0; std::cout << " Entering 2nd SCF, where EXX is updated" << std::endl; this->two_level_step++; @@ -202,7 +210,7 @@ namespace ModuleESolver // update exx and redo scf if (this->two_level_step == 0) { - XC_Functional::set_xc_type(GlobalC::ucell.atoms[0].ncpp.xc_func); + XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); } std::cout << " Updating EXX " << std::flush; @@ -226,19 +234,29 @@ namespace ModuleESolver } template - void ESolver_KS_LIP::after_all_runners() + void ESolver_KS_LIP::after_all_runners(UnitCell& ucell) { - ESolver_KS_PW::after_all_runners(); - + ESolver_KS_PW::after_all_runners(ucell); + #ifdef __LCAO if (PARAM.inp.out_mat_xc) { - ModuleIO::write_Vxc(PARAM.inp.nspin, PARAM.globalv.nlocal, - GlobalV::DRANK, *this->kspw_psi, GlobalC::ucell, this->sf, - *this->pw_wfc, *this->pw_rho, *this->pw_rhod, - GlobalC::ppcell.vloc, *this->pelec->charge, this->kv, this->pelec->wg + ModuleIO::write_Vxc(PARAM.inp.nspin, + PARAM.globalv.nlocal, + GlobalV::DRANK, + *this->kspw_psi, + ucell, + this->sf, + *this->pw_wfc, + *this->pw_rho, + *this->pw_rhod, + GlobalC::ppcell.vloc, + *this->pelec->charge, + this->kv, + this->pelec->wg #ifdef __EXX - , *this->exx_lip + , + *this->exx_lip #endif ); } diff --git a/source/module_esolver/esolver_ks_lcaopw.h b/source/module_esolver/esolver_ks_lcaopw.h index 68fc0a905b..6684281689 100644 --- a/source/module_esolver/esolver_ks_lcaopw.h +++ b/source/module_esolver/esolver_ks_lcaopw.h @@ -20,18 +20,21 @@ namespace ModuleESolver ~ESolver_KS_LIP(); - void before_all_runners(const Input_para& inp, UnitCell& cell) override; - void after_all_runners()override; + void before_all_runners(UnitCell& ucell, const Input_para& inp) override; + void after_all_runners(UnitCell& ucell) override; - protected: - virtual void iter_init(const int istep, const int iter) override; - virtual void iter_finish(const int istep, int& iter) override; + protected: + virtual void iter_init(UnitCell& ucell, const int istep, const int iter) override; + virtual void iter_finish(UnitCell& ucell, const int istep, int& iter) override; - /// All the other interfaces except this one are the same as ESolver_KS_PW. - virtual void hamilt2density_single(const int istep, const int iter, const double ethr) override; + /// All the other interfaces except this one are the same as ESolver_KS_PW. + virtual void hamilt2density_single(UnitCell& ucell, + const int istep, + const int iter, + const double ethr) override; - virtual void allocate_hamilt() override; - virtual void deallocate_hamilt() override; + virtual void allocate_hamilt() override; + virtual void deallocate_hamilt() override; #ifdef __EXX std::unique_ptr> exx_lip; diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index a053f2de95..30b1507265 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -129,10 +129,10 @@ void ESolver_KS_PW::deallocate_hamilt() } template -void ESolver_KS_PW::before_all_runners(const Input_para& inp, UnitCell& ucell) +void ESolver_KS_PW::before_all_runners(UnitCell& ucell, const Input_para& inp) { // 1) call before_all_runners() of ESolver_KS - ESolver_KS::before_all_runners(inp, ucell); + ESolver_KS::before_all_runners(ucell, inp); // 3) initialize ElecState, if (this->pelec == nullptr) @@ -237,23 +237,19 @@ void ESolver_KS_PW::before_all_runners(const Input_para& inp, UnitCel } template -void ESolver_KS_PW::before_scf(const int istep) +void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) { ModuleBase::TITLE("ESolver_KS_PW", "before_scf"); //! 1) call before_scf() of ESolver_FP - ESolver_FP::before_scf(istep); + ESolver_FP::before_scf(ucell, istep); - if (GlobalC::ucell.cell_parameter_updated) + if (ucell.cell_parameter_updated) { - GlobalC::ppcell.init_vnl(GlobalC::ucell, this->pw_rhod); + GlobalC::ppcell.init_vnl(ucell, this->pw_rhod); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "NON-LOCAL POTENTIAL"); - this->pw_wfc->initgrids(GlobalC::ucell.lat0, - GlobalC::ucell.latvec, - this->pw_wfc->nx, - this->pw_wfc->ny, - this->pw_wfc->nz); + this->pw_wfc->initgrids(ucell.lat0, ucell.latvec, this->pw_wfc->nx, this->pw_wfc->ny, this->pw_wfc->nz); this->pw_wfc->initparameters(false, PARAM.inp.ecutwfc, this->kv.get_nks(), this->kv.kvec_d.data()); @@ -261,14 +257,14 @@ void ESolver_KS_PW::before_scf(const int istep) this->p_wf_init->make_table(this->kv.get_nks(), &this->sf); } - if (GlobalC::ucell.ionic_position_updated) + if (ucell.ionic_position_updated) { - this->CE.update_all_dis(GlobalC::ucell); + this->CE.update_all_dis(ucell); this->CE.extrapolate_charge( #ifdef __MPI &(GlobalC::Pgrid), #endif - GlobalC::ucell, + ucell, this->pelec->charge, &this->sf, GlobalV::ofs_running, @@ -286,7 +282,7 @@ void ESolver_KS_PW::before_scf(const int istep) //---------------------------------------------------------- // about vdw, jiyy add vdwd3 and linpz add vdwd2 //---------------------------------------------------------- - auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp, &(GlobalV::ofs_running)); + auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp, &(GlobalV::ofs_running)); if (vdw_solver != nullptr) { this->pelec->f_en.evdw = vdw_solver->get_energy(); @@ -295,18 +291,18 @@ void ESolver_KS_PW::before_scf(const int istep) // calculate ewald energy if (!PARAM.inp.test_skip_ewald) { - this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(GlobalC::ucell, this->pw_rhod, this->sf.strucFac); + this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(ucell, this->pw_rhod, this->sf.strucFac); } //! cal_ux should be called before init_scf because //! the direction of ux is used in noncoline_rho if (PARAM.inp.nspin == 4) { - GlobalC::ucell.cal_ux(); + ucell.cal_ux(); } //! calculate the total local pseudopotential in real space - this->pelec->init_scf(istep, this->sf.strucFac, GlobalC::ucell.symm, (void*)this->pw_wfc); + this->pelec->init_scf(istep, this->sf.strucFac, ucell.symm, (void*)this->pw_wfc); //! output the initial charge density if (PARAM.inp.out_chg[0] == 2) @@ -316,13 +312,13 @@ void ESolver_KS_PW::before_scf(const int istep) std::stringstream ss; ss << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_CHG_INI.cube"; ModuleIO::write_vdata_palgrid(GlobalC::Pgrid, - this->pelec->charge->rho[is], - is, - PARAM.inp.nspin, - istep, - ss.str(), - this->pelec->eferm.ef, - &(GlobalC::ucell)); + this->pelec->charge->rho[is], + is, + PARAM.inp.nspin, + istep, + ss.str(), + this->pelec->eferm.ef, + &(ucell)); } } @@ -334,15 +330,15 @@ void ESolver_KS_PW::before_scf(const int istep) std::stringstream ss; ss << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_POT_INI.cube"; ModuleIO::write_vdata_palgrid(GlobalC::Pgrid, - this->pelec->pot->get_effective_v(is), - is, - PARAM.inp.nspin, - istep, - ss.str(), - 0.0, // efermi - &(GlobalC::ucell), - 11, // precsion - 0); // out_fermi + this->pelec->pot->get_effective_v(is), + is, + PARAM.inp.nspin, + istep, + ss.str(), + 0.0, // efermi + &(ucell), + 11, // precsion + 0); // out_fermi } } @@ -352,7 +348,7 @@ void ESolver_KS_PW::before_scf(const int istep) Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) { - srho.begin(is, *(this->pelec->charge), this->pw_rhod, GlobalC::ucell.symm); + srho.begin(is, *(this->pelec->charge), this->pw_rhod, ucell.symm); } // liuyu move here 2023-10-09 @@ -361,7 +357,7 @@ void ESolver_KS_PW::before_scf(const int istep) // projectors ModuleBase::matrix veff = this->pelec->pot->get_effective_v(); - GlobalC::ppcell.cal_effective_D(veff, this->pw_rhod, GlobalC::ucell); + GlobalC::ppcell.cal_effective_D(veff, this->pw_rhod, ucell); // after init_rho (in pelec->init_scf), we have rho now. // before hamilt2density, we update Hk and initialize psi @@ -388,10 +384,10 @@ void ESolver_KS_PW::before_scf(const int istep) } template -void ESolver_KS_PW::iter_init(const int istep, const int iter) +void ESolver_KS_PW::iter_init(UnitCell& ucell, const int istep, const int iter) { // call iter_init() of ESolver_KS - ESolver_KS::iter_init(istep, iter); + ESolver_KS::iter_init(ucell, istep, iter); if (iter == 1) { @@ -410,7 +406,10 @@ void ESolver_KS_PW::iter_init(const int istep, const int iter) // Temporary, it should be replaced by hsolver later. template -void ESolver_KS_PW::hamilt2density_single(const int istep, const int iter, const double ethr) +void ESolver_KS_PW::hamilt2density_single(UnitCell& ucell, + const int istep, + const int iter, + const double ethr) { ModuleBase::timer::tick("ESolver_KS_PW", "hamilt2density_single"); @@ -453,7 +452,7 @@ void ESolver_KS_PW::hamilt2density_single(const int istep, const int Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) { - srho.begin(is, *(this->pelec->charge), this->pw_rhod, GlobalC::ucell.symm); + srho.begin(is, *(this->pelec->charge), this->pw_rhod, ucell.symm); } // deband is calculated from "output" charge density calculated @@ -466,15 +465,15 @@ void ESolver_KS_PW::hamilt2density_single(const int istep, const int // Temporary, it should be rewritten with Hamilt class. template -void ESolver_KS_PW::update_pot(const int istep, const int iter) +void ESolver_KS_PW::update_pot(UnitCell& ucell, const int istep, const int iter) { if (!this->conv_esolver) { if (PARAM.inp.nspin == 4) { - GlobalC::ucell.cal_ux(); + ucell.cal_ux(); } - this->pelec->pot->update_from_charge(this->pelec->charge, &GlobalC::ucell); + this->pelec->pot->update_from_charge(this->pelec->charge, &ucell); this->pelec->f_en.descf = this->pelec->cal_delta_escf(); #ifdef __MPI MPI_Bcast(&(this->pelec->f_en.descf), 1, MPI_DOUBLE, 0, PARAPW_WORLD); @@ -487,10 +486,10 @@ void ESolver_KS_PW::update_pot(const int istep, const int iter) } template -void ESolver_KS_PW::iter_finish(const int istep, int& iter) +void ESolver_KS_PW::iter_finish(UnitCell& ucell, const int istep, int& iter) { // 1) Call iter_finish() of ESolver_KS - ESolver_KS::iter_finish(istep, iter); + ESolver_KS::iter_finish(ucell, istep, iter); // 2) Update USPP-related quantities // D in uspp need vloc, thus needs update when veff updated @@ -500,7 +499,7 @@ void ESolver_KS_PW::iter_finish(const int istep, int& iter) if (PARAM.globalv.use_uspp) { ModuleBase::matrix veff = this->pelec->pot->get_effective_v(); - GlobalC::ppcell.cal_effective_D(veff, this->pw_rhod, GlobalC::ucell); + GlobalC::ppcell.cal_effective_D(veff, this->pw_rhod, ucell); } // 3) Print out charge density @@ -521,26 +520,26 @@ void ESolver_KS_PW::iter_finish(const int istep, int& iter) } std::string fn = PARAM.globalv.global_out_dir + "/tmp_SPIN" + std::to_string(is + 1) + "_CHG.cube"; ModuleIO::write_vdata_palgrid(GlobalC::Pgrid, - data, - is, - PARAM.inp.nspin, - 0, - fn, - this->pelec->eferm.get_efval(is), - &(GlobalC::ucell), - 3, - 1); + data, + is, + PARAM.inp.nspin, + 0, + fn, + this->pelec->eferm.get_efval(is), + &(ucell), + 3, + 1); if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) { fn = PARAM.globalv.global_out_dir + "/tmp_SPIN" + std::to_string(is + 1) + "_TAU.cube"; ModuleIO::write_vdata_palgrid(GlobalC::Pgrid, - this->pelec->charge->kin_r_save[is], - is, - PARAM.inp.nspin, - 0, - fn, - this->pelec->eferm.get_efval(is), - &(GlobalC::ucell)); + this->pelec->charge->kin_r_save[is], + is, + PARAM.inp.nspin, + 0, + fn, + this->pelec->eferm.get_efval(is), + &(ucell)); } } } @@ -560,7 +559,7 @@ void ESolver_KS_PW::iter_finish(const int istep, int& iter) } template -void ESolver_KS_PW::after_scf(const int istep) +void ESolver_KS_PW::after_scf(UnitCell& ucell, const int istep) { // 1) calculate the kinetic energy density tau, sunliang 2024-09-18 if (PARAM.inp.out_elf[0] > 0) @@ -569,7 +568,7 @@ void ESolver_KS_PW::after_scf(const int istep) } // 2) call after_scf() of ESolver_KS - ESolver_KS::after_scf(istep); + ESolver_KS::after_scf(ucell, istep); // 3) output wavefunctions if (PARAM.inp.out_wfc_pw == 1 || PARAM.inp.out_wfc_pw == 2) @@ -606,7 +605,7 @@ void ESolver_KS_PW::after_scf(const int istep) this->pw_big->bz, this->pw_big->nbz, this->pelec->charge->ngmc, - &GlobalC::ucell, + &ucell, this->psi, this->pw_rhod, this->pw_wfc, @@ -649,9 +648,9 @@ double ESolver_KS_PW::cal_energy() } template -void ESolver_KS_PW::cal_force(ModuleBase::matrix& force) +void ESolver_KS_PW::cal_force(UnitCell& ucell, ModuleBase::matrix& force) { - Forces ff(GlobalC::ucell.nat); + Forces ff(ucell.nat); if (this->__kspw_psi != nullptr && PARAM.inp.precision == "single") { delete reinterpret_cast, Device>*>(this->__kspw_psi); @@ -663,18 +662,11 @@ void ESolver_KS_PW::cal_force(ModuleBase::matrix& force) : reinterpret_cast, Device>*>(this->kspw_psi); // Calculate forces - ff.cal_force(force, - *this->pelec, - this->pw_rhod, - &GlobalC::ucell.symm, - &this->sf, - &this->kv, - this->pw_wfc, - this->__kspw_psi); + ff.cal_force(force, *this->pelec, this->pw_rhod, &ucell.symm, &this->sf, &this->kv, this->pw_wfc, this->__kspw_psi); } template -void ESolver_KS_PW::cal_stress(ModuleBase::matrix& stress) +void ESolver_KS_PW::cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) { Stress_PW ss(this->pelec); if (this->__kspw_psi != nullptr && PARAM.inp.precision == "single") @@ -687,10 +679,10 @@ void ESolver_KS_PW::cal_stress(ModuleBase::matrix& stress) ? new psi::Psi, Device>(this->kspw_psi[0]) : reinterpret_cast, Device>*>(this->kspw_psi); ss.cal_stress(stress, - GlobalC::ucell, + ucell, &GlobalC::ppcell, this->pw_rhod, - &GlobalC::ucell.symm, + &ucell.symm, &this->sf, &this->kv, this->pw_wfc, @@ -707,7 +699,7 @@ void ESolver_KS_PW::cal_stress(ModuleBase::matrix& stress) } template -void ESolver_KS_PW::after_all_runners() +void ESolver_KS_PW::after_all_runners(UnitCell& ucell) { //! 1) Output information to screen GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; @@ -813,7 +805,7 @@ void ESolver_KS_PW::after_all_runners() << " a.u." << std::endl; } Numerical_Basis numerical_basis; - numerical_basis.output_overlap(this->psi[0], this->sf, this->kv, this->pw_wfc, GlobalC::ucell, i); + numerical_basis.output_overlap(this->psi[0], this->sf, this->kv, this->pw_wfc, ucell, i); } ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "BASIS OVERLAP (Q and S) GENERATION."); } @@ -828,7 +820,7 @@ void ESolver_KS_PW::after_all_runners() //! 7) Use Kubo-Greenwood method to compute conductivities if (PARAM.inp.cal_cond) { - EleCond elec_cond(&GlobalC::ucell, &this->kv, this->pelec, this->pw_wfc, this->psi, &GlobalC::ppcell); + EleCond elec_cond(&ucell, &this->kv, this->pelec, this->pw_wfc, this->psi, &GlobalC::ppcell); elec_cond.KG(PARAM.inp.cond_smear, PARAM.inp.cond_fwhm, PARAM.inp.cond_wcut, diff --git a/source/module_esolver/esolver_ks_pw.h b/source/module_esolver/esolver_ks_pw.h index e84b193719..e3806ab5f3 100644 --- a/source/module_esolver/esolver_ks_pw.h +++ b/source/module_esolver/esolver_ks_pw.h @@ -21,30 +21,30 @@ class ESolver_KS_PW : public ESolver_KS ~ESolver_KS_PW(); - void before_all_runners(const Input_para& inp, UnitCell& cell) override; + void before_all_runners(UnitCell& ucell, const Input_para& inp) override; double cal_energy() override; - void cal_force(ModuleBase::matrix& force) override; + void cal_force(UnitCell& ucell, ModuleBase::matrix& force) override; - void cal_stress(ModuleBase::matrix& stress) override; + void cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) override; - void after_all_runners() override; + void after_all_runners(UnitCell& ucell) override; protected: - virtual void before_scf(const int istep) override; + virtual void before_scf(UnitCell& ucell, const int istep) override; - virtual void iter_init(const int istep, const int iter) override; + virtual void iter_init(UnitCell& ucell, const int istep, const int iter) override; - virtual void update_pot(const int istep, const int iter) override; + virtual void update_pot(UnitCell& ucell, const int istep, const int iter) override; - virtual void iter_finish(const int istep, int& iter) override; + virtual void iter_finish(UnitCell& ucell, const int istep, int& iter) override; - virtual void after_scf(const int istep) override; + virtual void after_scf(UnitCell& ucell, const int istep) override; - virtual void others(const int istep) override; + virtual void others(UnitCell& ucell, const int istep) override; - virtual void hamilt2density_single(const int istep, const int iter, const double ethr) override; + virtual void hamilt2density_single(UnitCell& ucell, const int istep, const int iter, const double ethr) override; virtual void allocate_hamilt(); virtual void deallocate_hamilt(); diff --git a/source/module_esolver/esolver_lj.cpp b/source/module_esolver/esolver_lj.cpp index 6bfb943fa4..379aaced17 100644 --- a/source/module_esolver/esolver_lj.cpp +++ b/source/module_esolver/esolver_lj.cpp @@ -8,87 +8,86 @@ namespace ModuleESolver { - void ESolver_LJ::before_all_runners(const Input_para& inp, UnitCell& ucell) - { - ucell_ = &ucell; - lj_potential = 0; - lj_force.create(ucell.nat, 3); - lj_virial.create(3, 3); +void ESolver_LJ::before_all_runners(UnitCell& ucell, const Input_para& inp) +{ + ucell_ = &ucell; + lj_potential = 0; + lj_force.create(ucell.nat, 3); + lj_virial.create(3, 3); - ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU.cif", - ucell, - "# Generated by ABACUS ModuleIO::CifParser", - "data_?"); + ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU.cif", + ucell, + "# Generated by ABACUS ModuleIO::CifParser", + "data_?"); - // determine the maximum rcut and lj_rcut - rcut_search_radius(inp.mdp.lj_rcut); + // determine the maximum rcut and lj_rcut + rcut_search_radius(inp.mdp.lj_rcut); - // determine the LJ parameters - set_c6_c12(inp.mdp.lj_rule, inp.mdp.lj_epsilon, inp.mdp.lj_sigma); + // determine the LJ parameters + set_c6_c12(inp.mdp.lj_rule, inp.mdp.lj_epsilon, inp.mdp.lj_sigma); - // calculate the energy shift so that LJ energy is zero at rcut - cal_en_shift(inp.mdp.lj_eshift); - } + // calculate the energy shift so that LJ energy is zero at rcut + cal_en_shift(inp.mdp.lj_eshift); +} - void ESolver_LJ::runner(const int istep, UnitCell& ucell) +void ESolver_LJ::runner(UnitCell& ucell, const int istep) +{ + Grid_Driver grid_neigh(PARAM.inp.test_deconstructor, PARAM.inp.test_grid); + atom_arrange::search(PARAM.inp.search_pbc, + GlobalV::ofs_running, + grid_neigh, + ucell, + search_radius, + PARAM.inp.test_atom_input); + + double distance = 0.0; + int index = 0; + + // Important! potential, force, virial must be zero per step + lj_potential = 0; + lj_force.zero_out(); + lj_virial.zero_out(); + + ModuleBase::Vector3 tau1, tau2, dtau; + for (int it = 0; it < ucell.ntype; ++it) { - Grid_Driver grid_neigh(PARAM.inp.test_deconstructor, PARAM.inp.test_grid); - atom_arrange::search( - PARAM.inp.search_pbc, - GlobalV::ofs_running, - grid_neigh, - ucell, - search_radius, - PARAM.inp.test_atom_input); - - double distance=0.0; - int index = 0; - - // Important! potential, force, virial must be zero per step - lj_potential = 0; - lj_force.zero_out(); - lj_virial.zero_out(); - - ModuleBase::Vector3 tau1, tau2, dtau; - for (int it = 0; it < ucell.ntype; ++it) + Atom* atom1 = &ucell.atoms[it]; + for (int ia = 0; ia < atom1->na; ++ia) { - Atom* atom1 = &ucell.atoms[it]; - for (int ia = 0; ia < atom1->na; ++ia) + tau1 = atom1->tau[ia]; + grid_neigh.Find_atom(ucell, tau1, it, ia); + for (int ad = 0; ad < grid_neigh.getAdjacentNum(); ++ad) { - tau1 = atom1->tau[ia]; - grid_neigh.Find_atom(ucell, tau1, it, ia); - for (int ad = 0; ad < grid_neigh.getAdjacentNum(); ++ad) + tau2 = grid_neigh.getAdjacentTau(ad); + int it2 = grid_neigh.getType(ad); + dtau = (tau1 - tau2) * ucell.lat0; + distance = dtau.norm(); + if (distance < lj_rcut(it, it2)) { - tau2 = grid_neigh.getAdjacentTau(ad); - int it2 = grid_neigh.getType(ad); - dtau = (tau1 - tau2) * ucell.lat0; - distance = dtau.norm(); - if (distance < lj_rcut(it, it2)) - { - lj_potential += LJ_energy(distance, it, it2) - en_shift(it, it2); - ModuleBase::Vector3 f_ij = LJ_force(dtau, it, it2); - lj_force(index, 0) += f_ij.x; - lj_force(index, 1) += f_ij.y; - lj_force(index, 2) += f_ij.z; - LJ_virial(f_ij, dtau); - } + lj_potential += LJ_energy(distance, it, it2) - en_shift(it, it2); + ModuleBase::Vector3 f_ij = LJ_force(dtau, it, it2); + lj_force(index, 0) += f_ij.x; + lj_force(index, 1) += f_ij.y; + lj_force(index, 2) += f_ij.z; + LJ_virial(f_ij, dtau); } - index++; } + index++; } + } - lj_potential /= 2.0; - GlobalV::ofs_running << " final etot is " << std::setprecision(11) << lj_potential * ModuleBase::Ry_to_eV - << " eV" << std::endl; + lj_potential /= 2.0; + GlobalV::ofs_running << " final etot is " << std::setprecision(11) << lj_potential * ModuleBase::Ry_to_eV << " eV" + << std::endl; - // Post treatment for virial - for (int i = 0; i < 3; ++i) + // Post treatment for virial + for (int i = 0; i < 3; ++i) + { + for (int j = 0; j < 3; ++j) { - for (int j = 0; j < 3; ++j) - { - lj_virial(i, j) /= (2.0 * ucell.omega); - } + lj_virial(i, j) /= (2.0 * ucell.omega); } + } #ifdef __MPI atom_arrange::delete_vector( GlobalV::ofs_running, @@ -105,13 +104,13 @@ namespace ModuleESolver return lj_potential; } - void ESolver_LJ::cal_force(ModuleBase::matrix& force) + void ESolver_LJ::cal_force(UnitCell& ucell, 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(UnitCell& ucell, ModuleBase::matrix& stress) { stress = lj_virial; @@ -126,7 +125,7 @@ namespace ModuleESolver ModuleIO::print_stress("TOTAL-STRESS", stress, true, false); } - void ESolver_LJ::after_all_runners() + void ESolver_LJ::after_all_runners(UnitCell& ucell) { 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 668e39c287..8f50a1be35 100644 --- a/source/module_esolver/esolver_lj.h +++ b/source/module_esolver/esolver_lj.h @@ -14,17 +14,17 @@ namespace ModuleESolver classname = "ESolver_LJ"; } - void before_all_runners(const Input_para& inp, UnitCell& cell) override; + void before_all_runners(UnitCell& ucell, const Input_para& inp) override; - void runner(const int istep, UnitCell& cell) override; + void runner(UnitCell& cell, const int istep) override; double cal_energy() override; - void cal_force(ModuleBase::matrix& force) override; + void cal_force(UnitCell& ucell, ModuleBase::matrix& force) override; - void cal_stress(ModuleBase::matrix& stress) override; + void cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) override; - void after_all_runners() override; + void after_all_runners(UnitCell& ucell) override; private: diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp index 7b6fece3b8..8b9019e40c 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::before_all_runners(const Input_para& inp, UnitCell& ucell) +void ESolver_OF::before_all_runners(UnitCell& ucell, const Input_para& inp) { - ESolver_FP::before_all_runners(inp, ucell); + ESolver_FP::before_all_runners(ucell, inp); // save necessary parameters this->of_kinetic_ = inp.of_kinetic; @@ -69,6 +69,8 @@ void ESolver_OF::before_all_runners(const Input_para& inp, UnitCell& ucell) this->of_tolp_ = inp.of_tolp; this->max_iter_ = inp.scf_nmax; this->dV_ = ucell.omega / this->pw_rho->nxyz; + this->bound_cal_potential_ + = std::bind(&ESolver_OF::cal_potential, this, std::placeholders::_1, std::placeholders::_2, std::ref(ucell)); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SETUP UNITCELL"); @@ -119,7 +121,7 @@ void ESolver_OF::before_all_runners(const Input_para& inp, UnitCell& ucell) this->init_elecstate(ucell); // calculate the total local pseudopotential in real space - this->pelec->init_scf(0, sf.strucFac, GlobalC::ucell.symm); // atomic_rho, v_of_rho, set_vrs + this->pelec->init_scf(0, sf.strucFac, ucell.symm); // atomic_rho, v_of_rho, set_vrs // liuyu move here 2023-10-09 // D in uspp need vloc, thus behind init_scf() @@ -153,7 +155,7 @@ void ESolver_OF::before_all_runners(const Input_para& inp, UnitCell& ucell) this->allocate_array(); } -void ESolver_OF::runner(int istep, UnitCell& ucell) +void ESolver_OF::runner(UnitCell& ucell, const int istep) { ModuleBase::timer::tick("ESolver_OF", "runner"); // get Ewald energy, initial rho and phi if necessary @@ -201,7 +203,7 @@ void ESolver_OF::runner(int istep, UnitCell& ucell) void ESolver_OF::before_opt(const int istep, UnitCell& ucell) { //! 1) call before_scf() of ESolver_FP - ESolver_FP::before_scf(istep); + ESolver_FP::before_scf(ucell, istep); if (ucell.cell_parameter_updated) { @@ -256,7 +258,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell) GlobalV::ofs_warning); } - this->pelec->init_scf(istep, sf.strucFac, GlobalC::ucell.symm); + this->pelec->init_scf(istep, sf.strucFac, ucell.symm); // calculate ewald energy this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(ucell, this->pw_rho, sf.strucFac); @@ -264,7 +266,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell) Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) { - srho.begin(is, *(pelec->charge), this->pw_rho, GlobalC::ucell.symm); + srho.begin(is, *(pelec->charge), this->pw_rho, ucell.symm); } for (int is = 0; is < PARAM.inp.nspin; ++is) @@ -359,7 +361,7 @@ void ESolver_OF::update_potential(UnitCell& ucell) void ESolver_OF::optimize(UnitCell& ucell) { // (1) get |d0> with optimization algorithm - this->get_direction(); + this->get_direction(ucell); // initialize temp_phi and temp_rho used in line search double** ptemp_phi = new double*[PARAM.inp.nspin]; for (int is = 0; is < PARAM.inp.nspin; ++is) @@ -415,7 +417,7 @@ void ESolver_OF::update_rho() // Symmetry_rho srho; // for (int is = 0; is < PARAM.inp.nspin; is++) // { - // srho.begin(is, *(pelec->charge), this->pw_rho, GlobalC::Pgrid, GlobalC::ucell.symm); + // srho.begin(is, *(pelec->charge), this->pw_rho, GlobalC::Pgrid, ucell.symm); // for (int ibs = 0; ibs < this->pw_rho->nrxx; ++ibs) // { // this->pphi_[is][ibs] = sqrt(pelec->charge->rho[is][ibs]); @@ -494,13 +496,13 @@ void ESolver_OF::after_opt(const int istep, UnitCell& ucell) } // 2) call after_scf() of ESolver_FP - ESolver_FP::after_scf(istep); + ESolver_FP::after_scf(ucell, istep); } /** * @brief Output the FINAL_ETOT */ -void ESolver_OF::after_all_runners() +void ESolver_OF::after_all_runners(UnitCell& ucell) { GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; @@ -539,10 +541,10 @@ double ESolver_OF::cal_energy() * * @param [out] force */ -void ESolver_OF::cal_force(ModuleBase::matrix& force) +void ESolver_OF::cal_force(UnitCell& ucell, ModuleBase::matrix& force) { - Forces ff(GlobalC::ucell.nat); - ff.cal_force(force, *pelec, this->pw_rho, &GlobalC::ucell.symm, &sf); + Forces ff(ucell.nat); + ff.cal_force(force, *pelec, this->pw_rho, &ucell.symm, &sf); } /** @@ -550,13 +552,13 @@ void ESolver_OF::cal_force(ModuleBase::matrix& force) * * @param [out] stress */ -void ESolver_OF::cal_stress(ModuleBase::matrix& stress) +void ESolver_OF::cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) { ModuleBase::matrix kinetic_stress_; kinetic_stress_.create(3, 3); this->kinetic_stress(kinetic_stress_); OF_Stress_PW ss(this->pelec, this->pw_rho); - ss.cal_stress(stress, kinetic_stress_, GlobalC::ucell, &GlobalC::ucell.symm, &sf, &kv); + ss.cal_stress(stress, kinetic_stress_, ucell, &ucell.symm, &sf, &kv); } } // namespace ModuleESolver diff --git a/source/module_esolver/esolver_of.h b/source/module_esolver/esolver_of.h index 183976df8f..081a077606 100644 --- a/source/module_esolver/esolver_of.h +++ b/source/module_esolver/esolver_of.h @@ -19,17 +19,17 @@ class ESolver_OF : public ESolver_FP ESolver_OF(); ~ESolver_OF(); - virtual void before_all_runners(const Input_para& inp, UnitCell& ucell) override; + virtual void before_all_runners(UnitCell& ucell, const Input_para& inp) override; - virtual void runner(const int istep, UnitCell& ucell) override; + virtual void runner(UnitCell& ucell, const int istep) override; - virtual void after_all_runners() override; + virtual void after_all_runners(UnitCell& ucell) override; virtual double cal_energy() override; - virtual void cal_force(ModuleBase::matrix& force) override; + virtual void cal_force(UnitCell& ucell, ModuleBase::matrix& force) override; - virtual void cal_stress(ModuleBase::matrix& stress) override; + virtual void cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) override; private: // ======================= variables ========================== @@ -94,7 +94,9 @@ class ESolver_OF : public ESolver_FP void allocate_array(); // --------------------- calculate physical qualities --------------- - void cal_potential(double* ptemp_phi, double* rdLdphi); + std::function bound_cal_potential_; + void cal_potential_wrapper(double* ptemp_phi, double* rdLdphi); + void cal_potential(double* ptemp_phi, double* rdLdphi, UnitCell& ucell); void cal_dEdtheta(double** ptemp_phi, Charge* temp_rho, UnitCell& ucell, double* ptheta, double* rdEdtheta); double cal_mu(double* pphi, double* pdEdphi, double nelec); @@ -123,7 +125,7 @@ class ESolver_OF : public ESolver_FP // ---------------------- interfaces to optimization methods -------- void init_opt(); - void get_direction(); + void get_direction(UnitCell& ucell); void get_step_length(double* dEdtheta, double** ptemp_phi, UnitCell& ucell); }; } // namespace ModuleESolver diff --git a/source/module_esolver/esolver_of_interface.cpp b/source/module_esolver/esolver_of_interface.cpp index 118168a2a3..a92eb0290a 100644 --- a/source/module_esolver/esolver_of_interface.cpp +++ b/source/module_esolver/esolver_of_interface.cpp @@ -257,11 +257,16 @@ void ESolver_OF::init_opt() } } +void ESolver_OF::cal_potential_wrapper(double* ptemp_phi, double* rdLdphi) +{ + this->bound_cal_potential_(ptemp_phi, rdLdphi); +} + /** * @brief [Interface to opt] * Call optimization methods to get the optimization direction */ -void ESolver_OF::get_direction() +void ESolver_OF::get_direction(UnitCell& ucell) { for (int is = 0; is < PARAM.inp.nspin; ++is) { @@ -273,7 +278,7 @@ void ESolver_OF::get_direction() this->flag_, this->pdirect_[is], this, - &ESolver_OF::cal_potential); + &ESolver_OF::cal_potential_wrapper); } else if (this->of_method_ == "cg1") { diff --git a/source/module_esolver/esolver_of_tool.cpp b/source/module_esolver/esolver_of_tool.cpp index 6e81ad33a6..9d655d8e61 100644 --- a/source/module_esolver/esolver_of_tool.cpp +++ b/source/module_esolver/esolver_of_tool.cpp @@ -110,7 +110,7 @@ void ESolver_OF::allocate_array() * @param [in] ptemp_phi phi * @param [out] rdLdphi dL/dphi */ -void ESolver_OF::cal_potential(double* ptemp_phi, double* rdLdphi) +void ESolver_OF::cal_potential(double* ptemp_phi, double* rdLdphi, UnitCell& ucell) { double** dEdtemp_phi = new double*[PARAM.inp.nspin]; double** temp_phi = new double*[PARAM.inp.nspin]; @@ -134,9 +134,9 @@ void ESolver_OF::cal_potential(double* ptemp_phi, double* rdLdphi) if (PARAM.inp.nspin == 4) { - GlobalC::ucell.cal_ux(); + ucell.cal_ux(); } - this->pelec->pot->update_from_charge(this->ptemp_rho_, &GlobalC::ucell); + this->pelec->pot->update_from_charge(this->ptemp_rho_, &ucell); ModuleBase::matrix& vr_eff = this->pelec->pot->get_effective_v(); this->kinetic_potential(this->ptemp_rho_->rho, temp_phi, vr_eff); diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp index 27a7757ed8..9cb35acbc3 100644 --- a/source/module_esolver/esolver_sdft_pw.cpp +++ b/source/module_esolver/esolver_sdft_pw.cpp @@ -43,14 +43,14 @@ ESolver_SDFT_PW::~ESolver_SDFT_PW() } template -void ESolver_SDFT_PW::before_all_runners(const Input_para& inp, UnitCell& ucell) +void ESolver_SDFT_PW::before_all_runners(UnitCell& ucell, const Input_para& inp) { // 1) initialize parameters from int Input class this->nche_sto = inp.nche_sto; this->method_sto = inp.method_sto; // 2) run "before_all_runners" in ESolver_KS - ESolver_KS_PW::before_all_runners(inp, ucell); + ESolver_KS_PW::before_all_runners(ucell, inp); // 3) initialize the stochastic wave functions this->stowf.init(&this->kv, this->pw_wfc->npwk_max); @@ -92,9 +92,9 @@ void ESolver_SDFT_PW::before_all_runners(const Input_para& inp, UnitC } template -void ESolver_SDFT_PW::before_scf(const int istep) +void ESolver_SDFT_PW::before_scf(UnitCell& ucell, const int istep) { - ESolver_KS_PW::before_scf(istep); + ESolver_KS_PW::before_scf(ucell, istep); delete reinterpret_cast*>(this->p_hamilt); this->p_hamilt = new hamilt::HamiltSdftPW(this->pelec->pot, this->pw_wfc, @@ -111,21 +111,21 @@ void ESolver_SDFT_PW::before_scf(const int istep) } template -void ESolver_SDFT_PW::iter_finish(const int istep, int& iter) +void ESolver_SDFT_PW::iter_finish(UnitCell& ucell, const int istep, int& iter) { // call iter_finish() of ESolver_KS - ESolver_KS::iter_finish(istep, iter); + ESolver_KS::iter_finish(ucell, istep, iter); } template -void ESolver_SDFT_PW::after_scf(const int istep) +void ESolver_SDFT_PW::after_scf(UnitCell& ucell, const int istep) { // 1) call after_scf() of ESolver_KS_PW - ESolver_KS_PW::after_scf(istep); + ESolver_KS_PW::after_scf(ucell, istep); } template -void ESolver_SDFT_PW::hamilt2density_single(int istep, int iter, double ethr) +void ESolver_SDFT_PW::hamilt2density_single(UnitCell& ucell, int istep, int iter, double ethr) { ModuleBase::TITLE("ESolver_SDFT_PW", "hamilt2density"); ModuleBase::timer::tick("ESolver_SDFT_PW", "hamilt2density"); @@ -183,7 +183,7 @@ void ESolver_SDFT_PW::hamilt2density_single(int istep, int iter, doub Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) { - srho.begin(is, *(this->pelec->charge), this->pw_rho, GlobalC::ucell.symm); + srho.begin(is, *(this->pelec->charge), this->pw_rho, ucell.symm); } this->pelec->f_en.deband = this->pelec->cal_delta_eband(); } @@ -209,31 +209,31 @@ double ESolver_SDFT_PW::cal_energy() } template -void ESolver_SDFT_PW::cal_force(ModuleBase::matrix& force) +void ESolver_SDFT_PW::cal_force(UnitCell& ucell, ModuleBase::matrix& force) { - Sto_Forces ff(GlobalC::ucell.nat); + Sto_Forces ff(ucell.nat); ff.cal_stoforce(force, *this->pelec, this->pw_rho, - &GlobalC::ucell.symm, + &ucell.symm, &this->sf, &this->kv, this->pw_wfc, GlobalC::ppcell, - GlobalC::ucell, + ucell, *this->kspw_psi, this->stowf); } template -void ESolver_SDFT_PW::cal_stress(ModuleBase::matrix& stress) +void ESolver_SDFT_PW::cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) { Sto_Stress_PW ss; ss.cal_stress(stress, *this->pelec, this->pw_rho, - &GlobalC::ucell.symm, + &ucell.symm, &this->sf, &this->kv, this->pw_wfc, @@ -241,11 +241,11 @@ void ESolver_SDFT_PW::cal_stress(ModuleBase::matrix& stress) this->stowf, this->pelec->charge, &GlobalC::ppcell, - GlobalC::ucell); + ucell); } template -void ESolver_SDFT_PW::after_all_runners() +void ESolver_SDFT_PW::after_all_runners(UnitCell& ucell) { GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; GlobalV::ofs_running << std::setprecision(16); @@ -255,7 +255,7 @@ void ESolver_SDFT_PW::after_all_runners() } template <> -void ESolver_SDFT_PW, base_device::DEVICE_CPU>::after_all_runners() +void ESolver_SDFT_PW, base_device::DEVICE_CPU>::after_all_runners(UnitCell& ucell) { GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; @@ -285,7 +285,7 @@ void ESolver_SDFT_PW, base_device::DEVICE_CPU>::after_all_r // sKG cost memory, and it should be placed at the end of the program if (PARAM.inp.cal_cond) { - Sto_EleCond sto_elecond(&GlobalC::ucell, + Sto_EleCond sto_elecond(&ucell, &this->kv, this->pelec, this->pw_wfc, @@ -306,7 +306,7 @@ void ESolver_SDFT_PW, base_device::DEVICE_CPU>::after_all_r } template -void ESolver_SDFT_PW::others(const int istep) +void ESolver_SDFT_PW::others(UnitCell& ucell, const int istep) { ModuleBase::TITLE("ESolver_SDFT_PW", "others"); diff --git a/source/module_esolver/esolver_sdft_pw.h b/source/module_esolver/esolver_sdft_pw.h index 63337942fc..3287fc9ecb 100644 --- a/source/module_esolver/esolver_sdft_pw.h +++ b/source/module_esolver/esolver_sdft_pw.h @@ -19,13 +19,13 @@ class ESolver_SDFT_PW : public ESolver_KS_PW ESolver_SDFT_PW(); ~ESolver_SDFT_PW(); - void before_all_runners(const Input_para& inp, UnitCell& cell) override; + void before_all_runners(UnitCell& ucell, const Input_para& inp) override; double cal_energy() override; - void cal_force(ModuleBase::matrix& force) override; + void cal_force(UnitCell& ucell, ModuleBase::matrix& force) override; - void cal_stress(ModuleBase::matrix& stress) override; + void cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) override; public: Stochastic_WF stowf; @@ -33,17 +33,17 @@ class ESolver_SDFT_PW : public ESolver_KS_PW hamilt::HamiltSdftPW* p_hamilt_sto = nullptr; protected: - virtual void before_scf(const int istep) override; + virtual void before_scf(UnitCell& ucell, const int istep) override; - virtual void hamilt2density_single(const int istep, const int iter, const double ethr) override; + virtual void hamilt2density_single(UnitCell& ucell, const int istep, const int iter, const double ethr) override; - virtual void others(const int istep) override; + virtual void others(UnitCell& ucell, const int istep) override; - virtual void iter_finish(const int istep, int& iter) override; + virtual void iter_finish(UnitCell& ucell, const int istep, int& iter) override; - virtual void after_scf(const int istep) override; + virtual void after_scf(UnitCell& ucell, const int istep) override; - virtual void after_all_runners() override; + virtual void after_all_runners(UnitCell& ucell) override; private: int nche_sto; ///< norder of Chebyshev diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index 0c376bc421..b643ffe4c6 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -39,21 +39,21 @@ namespace ModuleESolver { template -void ESolver_KS_LCAO::before_scf(const int istep) +void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) { ModuleBase::TITLE("ESolver_KS_LCAO", "before_scf"); //! 1) call before_scf() of ESolver_FP - ESolver_FP::before_scf(istep); + ESolver_FP::before_scf(ucell, istep); - if (GlobalC::ucell.ionic_position_updated) + if (ucell.ionic_position_updated) { - this->CE.update_all_dis(GlobalC::ucell); + this->CE.update_all_dis(ucell); this->CE.extrapolate_charge( #ifdef __MPI &(GlobalC::Pgrid), #endif - GlobalC::ucell, + ucell, this->pelec->charge, &(this->sf), GlobalV::ofs_running, @@ -63,7 +63,7 @@ void ESolver_KS_LCAO::before_scf(const int istep) //---------------------------------------------------------- // about vdw, jiyy add vdwd3 and linpz add vdwd2 //---------------------------------------------------------- - auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp, &(GlobalV::ofs_running)); + auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp, &(GlobalV::ofs_running)); if (vdw_solver != nullptr) { this->pelec->f_en.evdw = vdw_solver->get_energy(); @@ -74,13 +74,13 @@ void ESolver_KS_LCAO::before_scf(const int istep) double search_radius = atom_arrange::set_sr_NL(GlobalV::ofs_running, PARAM.inp.out_level, orb_.get_rcutmax_Phi(), - GlobalC::ucell.infoNL.get_rcutmax_Beta(), + ucell.infoNL.get_rcutmax_Beta(), PARAM.globalv.gamma_only_local); atom_arrange::search(PARAM.inp.search_pbc, GlobalV::ofs_running, GlobalC::GridD, - GlobalC::ucell, + ucell, search_radius, PARAM.inp.test_atom_input); @@ -91,7 +91,7 @@ void ESolver_KS_LCAO::before_scf(const int istep) std::vector> dpsi_u; std::vector> d2psi_u; - Gint_Tools::init_orb(dr_uniform, rcuts, GlobalC::ucell, orb_, psi_u, dpsi_u, d2psi_u); + Gint_Tools::init_orb(dr_uniform, rcuts, ucell, orb_, psi_u, dpsi_u, d2psi_u); this->GridT.set_pbc_grid(this->pw_rho->nx, this->pw_rho->ny, @@ -108,7 +108,7 @@ void ESolver_KS_LCAO::before_scf(const int istep) this->pw_rho->ny, this->pw_rho->nplane, this->pw_rho->startz_current, - GlobalC::ucell, + ucell, GlobalC::GridD, dr_uniform, rcuts, @@ -211,14 +211,14 @@ void ESolver_KS_LCAO::before_scf(const int istep) const Parallel_Orbitals* pv = &this->pv; // build and save at beginning GlobalC::ld.build_psialpha(PARAM.inp.cal_force, - GlobalC::ucell, + ucell, orb_, GlobalC::GridD, *(two_center_bundle_.overlap_orb_alpha)); if (PARAM.inp.deepks_out_unittest) { - GlobalC::ld.check_psialpha(PARAM.inp.cal_force, GlobalC::ucell, orb_, GlobalC::GridD); + GlobalC::ld.check_psialpha(PARAM.inp.cal_force, ucell, orb_, GlobalC::GridD); } } #endif @@ -231,7 +231,7 @@ void ESolver_KS_LCAO::before_scf(const int istep) PARAM.inp.alpha_trial, PARAM.inp.sccut, PARAM.inp.sc_drop_thr, - GlobalC::ucell, + ucell, &(this->pv), PARAM.inp.nspin, this->kv, @@ -246,7 +246,7 @@ void ESolver_KS_LCAO::before_scf(const int istep) //========================================================= if (PARAM.inp.nspin == 4) { - GlobalC::ucell.cal_ux(); + ucell.cal_ux(); } // Peize Lin add 2016-12-03 @@ -255,16 +255,16 @@ void ESolver_KS_LCAO::before_scf(const int istep) { if (GlobalC::exx_info.info_ri.real_number) { - this->exd->exx_beforescf(istep, this->kv, *this->p_chgmix, GlobalC::ucell, orb_); + this->exd->exx_beforescf(istep, this->kv, *this->p_chgmix, ucell, orb_); } else { - this->exc->exx_beforescf(istep, this->kv, *this->p_chgmix, GlobalC::ucell, orb_); + this->exc->exx_beforescf(istep, this->kv, *this->p_chgmix, ucell, orb_); } } #endif // __EXX - this->pelec->init_scf(istep, this->sf.strucFac, GlobalC::ucell.symm); + this->pelec->init_scf(istep, this->sf.strucFac, ucell.symm); //! output the initial charge density if (PARAM.inp.out_chg[0] == 2) @@ -274,13 +274,13 @@ void ESolver_KS_LCAO::before_scf(const int istep) std::stringstream ss; ss << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_CHG_INI.cube"; ModuleIO::write_vdata_palgrid(GlobalC::Pgrid, - this->pelec->charge->rho[is], - is, - PARAM.inp.nspin, - istep, - ss.str(), - this->pelec->eferm.ef, - &(GlobalC::ucell)); + this->pelec->charge->rho[is], + is, + PARAM.inp.nspin, + istep, + ss.str(), + this->pelec->eferm.ef, + &(ucell)); } } @@ -292,15 +292,15 @@ void ESolver_KS_LCAO::before_scf(const int istep) std::stringstream ss; ss << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_POT_INI.cube"; ModuleIO::write_vdata_palgrid(GlobalC::Pgrid, - this->pelec->pot->get_effective_v(is), - is, - PARAM.inp.nspin, - istep, - ss.str(), - 0.0, // efermi - &(GlobalC::ucell), - 11, // precsion - 0); // out_fermi + this->pelec->pot->get_effective_v(is), + is, + PARAM.inp.nspin, + istep, + ss.str(), + 0.0, // efermi + &(ucell), + 11, // precsion + 0); // out_fermi } } @@ -324,11 +324,11 @@ void ESolver_KS_LCAO::before_scf(const int istep) std::string zipname = "output_DM0.npz"; elecstate::DensityMatrix* dm = dynamic_cast*>(this->pelec)->get_DM(); - ModuleIO::read_mat_npz(&(this->pv), GlobalC::ucell, zipname, *(dm->get_DMR_pointer(1))); + ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(1))); if (PARAM.inp.nspin == 2) { zipname = "output_DM1.npz"; - ModuleIO::read_mat_npz(&(this->pv), GlobalC::ucell, zipname, *(dm->get_DMR_pointer(2))); + ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(2))); } this->pelec->calculate_weights(); @@ -339,15 +339,15 @@ void ESolver_KS_LCAO::before_scf(const int istep) { std::string fn = PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_CHG.cube"; ModuleIO::write_vdata_palgrid(GlobalC::Pgrid, - this->pelec->charge->rho[is], - is, - PARAM.inp.nspin, - istep, - fn, - this->pelec->eferm.get_efval(is), - &(GlobalC::ucell), - 3, - 1); + this->pelec->charge->rho[is], + is, + PARAM.inp.nspin, + istep, + fn, + this->pelec->eferm.get_efval(is), + &(ucell), + 3, + 1); } return; @@ -358,14 +358,14 @@ void ESolver_KS_LCAO::before_scf(const int istep) Symmetry_rho srho; for (int is = 0; is < PARAM.inp.nspin; is++) { - srho.begin(is, *(this->pelec->charge), this->pw_rho, GlobalC::ucell.symm); + srho.begin(is, *(this->pelec->charge), this->pw_rho, ucell.symm); } // 1. calculate ewald energy. // mohan update 2021-02-25 if (!PARAM.inp.test_skip_ewald) { - this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(GlobalC::ucell, this->pw_rho, this->sf.strucFac); + this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(ucell, this->pw_rho, this->sf.strucFac); } this->p_hamilt->non_first_scf = istep; diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index 0f1ae92640..b4c05b2af7 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/source/module_esolver/lcao_others.cpp @@ -37,7 +37,7 @@ namespace ModuleESolver { template -void ESolver_KS_LCAO::others(const int istep) +void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) { ModuleBase::TITLE("ESolver_KS_LCAO", "others"); ModuleBase::timer::tick("ESolver_KS_LCAO", "others"); @@ -62,7 +62,7 @@ void ESolver_KS_LCAO::others(const int istep) atom_arrange::search(PARAM.inp.search_pbc, GlobalV::ofs_running, GlobalC::GridD, - GlobalC::ucell, + ucell, search_radius, PARAM.inp.test_atom_input, true); @@ -75,13 +75,13 @@ void ESolver_KS_LCAO::others(const int istep) double search_radius = atom_arrange::set_sr_NL(GlobalV::ofs_running, PARAM.inp.out_level, orb_.get_rcutmax_Phi(), - GlobalC::ucell.infoNL.get_rcutmax_Beta(), + ucell.infoNL.get_rcutmax_Beta(), PARAM.globalv.gamma_only_local); atom_arrange::search(PARAM.inp.search_pbc, GlobalV::ofs_running, GlobalC::GridD, - GlobalC::ucell, + ucell, search_radius, PARAM.inp.test_atom_input); @@ -92,7 +92,7 @@ void ESolver_KS_LCAO::others(const int istep) std::vector> dpsi_u; std::vector> d2psi_u; - Gint_Tools::init_orb(dr_uniform, rcuts, GlobalC::ucell, orb_, psi_u, dpsi_u, d2psi_u); + Gint_Tools::init_orb(dr_uniform, rcuts, ucell, orb_, psi_u, dpsi_u, d2psi_u); this->GridT.set_pbc_grid(this->pw_rho->nx, this->pw_rho->ny, @@ -109,7 +109,7 @@ void ESolver_KS_LCAO::others(const int istep) this->pw_rho->ny, this->pw_rho->nplane, this->pw_rho->startz_current, - GlobalC::ucell, + ucell, GlobalC::GridD, dr_uniform, rcuts, @@ -212,14 +212,14 @@ void ESolver_KS_LCAO::others(const int istep) const Parallel_Orbitals* pv = &this->pv; // build and save at beginning GlobalC::ld.build_psialpha(PARAM.inp.cal_force, - GlobalC::ucell, + ucell, orb_, GlobalC::GridD, *(two_center_bundle_.overlap_orb_alpha)); if (PARAM.inp.deepks_out_unittest) { - GlobalC::ld.check_psialpha(PARAM.inp.cal_force, GlobalC::ucell, orb_, GlobalC::GridD); + GlobalC::ld.check_psialpha(PARAM.inp.cal_force, ucell, orb_, GlobalC::GridD); } } #endif @@ -232,7 +232,7 @@ void ESolver_KS_LCAO::others(const int istep) PARAM.inp.alpha_trial, PARAM.inp.sccut, PARAM.inp.sc_drop_thr, - GlobalC::ucell, + ucell, &(this->pv), PARAM.inp.nspin, this->kv, @@ -247,11 +247,11 @@ void ESolver_KS_LCAO::others(const int istep) //========================================================= if (PARAM.inp.nspin == 4) { - GlobalC::ucell.cal_ux(); + ucell.cal_ux(); } // pelec should be initialized before these calculations - this->pelec->init_scf(istep, this->sf.strucFac, GlobalC::ucell.symm); + this->pelec->init_scf(istep, this->sf.strucFac, ucell.symm); // self consistent calculations for electronic ground state if (cal_type == "get_pchg") { @@ -280,7 +280,7 @@ void ESolver_KS_LCAO::others(const int istep) PARAM.globalv.nlocal, PARAM.globalv.global_out_dir, GlobalV::ofs_warning, - &GlobalC::ucell, + &ucell, &GlobalC::GridD, this->kv); } @@ -309,7 +309,7 @@ void ESolver_KS_LCAO::others(const int istep) PARAM.globalv.nlocal, PARAM.globalv.global_out_dir, GlobalV::ofs_warning, - &GlobalC::ucell, + &ucell, &GlobalC::GridD, this->kv, PARAM.inp.if_separate_k, diff --git a/source/module_esolver/pw_others.cpp b/source/module_esolver/pw_others.cpp index 5db594c225..c36775aa07 100644 --- a/source/module_esolver/pw_others.cpp +++ b/source/module_esolver/pw_others.cpp @@ -48,9 +48,9 @@ namespace ModuleESolver { - template -void ESolver_KS_PW::others(const int istep) { +void ESolver_KS_PW::others(UnitCell& ucell, const int istep) +{ ModuleBase::TITLE("ESolver_KS_PW", "others"); const std::string cal_type = PARAM.inp.calculation; diff --git a/source/module_esolver/test/esolver_dp_test.cpp b/source/module_esolver/test/esolver_dp_test.cpp index 8a2c99aade..769a04fa00 100644 --- a/source/module_esolver/test/esolver_dp_test.cpp +++ b/source/module_esolver/test/esolver_dp_test.cpp @@ -58,7 +58,7 @@ class ESolverDPTest : public ::testing::Test ucell.atom_label = new std::string[2]; ucell.atom_label[0] = "Cu"; ucell.atom_label[1] = "Al"; - esolver->before_all_runners(inp, ucell); + esolver->before_all_runners(ucell, inp); } void TearDown() override @@ -107,7 +107,7 @@ TEST_F(ESolverDPTest, RunWarningQuit) int istep = 0; testing::internal::CaptureStdout(); - EXPECT_EXIT(esolver->runner(istep, ucell), ::testing::ExitedWithCode(0), ""); + EXPECT_EXIT(esolver->runner(ucell, istep), ::testing::ExitedWithCode(0), ""); std::string output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("Please recompile with -D__DPMD")); } @@ -135,7 +135,7 @@ TEST_F(ESolverDPTest, CalForce) } } - esolver->cal_force(force); + esolver->cal_force(ucell, force); // Check the results for (int i = 0; i < ucell.nat; ++i) @@ -159,7 +159,7 @@ TEST_F(ESolverDPTest, CalStress) } } - esolver->cal_stress(stress); + esolver->cal_stress(ucell, stress); // Check the results for (int i = 0; i < 3; ++i) @@ -178,7 +178,7 @@ TEST_F(ESolverDPTest, Postprocess) // Check the results GlobalV::ofs_running.open("log"); - esolver->after_all_runners(); + esolver->after_all_runners(ucell); GlobalV::ofs_running.close(); std::string expected_output = "\n\n --------------------------------------------\n !FINAL_ETOT_IS 133.3358404 eV\n " diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp index 4e6d9dc526..863985bb1e 100644 --- a/source/module_lr/esolver_lrtd_lcao.cpp +++ b/source/module_lr/esolver_lrtd_lcao.cpp @@ -249,7 +249,7 @@ LR::ESolver_LR::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu std::transform(xc_kernel.begin(), xc_kernel.end(), xc_kernel.begin(), tolower); // necessary steps in ESolver_FP - ESolver_FP::before_all_runners(inp, ucell); + ESolver_FP::before_all_runners(ucell, inp); this->pelec = new elecstate::ElecStateLCAO(); // necessary steps in ESolver_KS::before_all_runners : symmetry and k-points @@ -407,7 +407,7 @@ LR::ESolver_LR::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu } template -void LR::ESolver_LR::runner(int istep, UnitCell& cell) +void LR::ESolver_LR::runner(UnitCell& ucell, const int istep) { ModuleBase::TITLE("ESolver_LR", "runner"); ModuleBase::timer::tick("ESolver_LR", "runner"); @@ -488,7 +488,7 @@ void LR::ESolver_LR::runner(int istep, UnitCell& cell) } template -void LR::ESolver_LR::after_all_runners() +void LR::ESolver_LR::after_all_runners(UnitCell& ucell) { ModuleBase::TITLE("ESolver_LR", "after_all_runners"); if (input.ri_hartree_benchmark != "none") { return; } //no need to calculate the spectrum in the benchmark routine diff --git a/source/module_lr/esolver_lrtd_lcao.h b/source/module_lr/esolver_lrtd_lcao.h index 8ae33ef05a..bc3db43e58 100644 --- a/source/module_lr/esolver_lrtd_lcao.h +++ b/source/module_lr/esolver_lrtd_lcao.h @@ -37,15 +37,15 @@ namespace LR ///input: input, call, basis(LCAO), psi(ground state), elecstate // initialize sth. independent of the ground state - virtual void before_all_runners(const Input_para& inp, UnitCell& cell) override {}; - virtual void runner(int istep, UnitCell& ucell) override; - virtual void after_all_runners() override; + virtual void before_all_runners(UnitCell& ucell, const Input_para& inp) override {}; + virtual void runner(UnitCell& ucell, int istep) override; + virtual void after_all_runners(UnitCell& ucell) override; virtual double cal_energy() override { return 0.0; }; - virtual void cal_force(ModuleBase::matrix& force) override {}; - virtual void cal_stress(ModuleBase::matrix& stress) override {}; + virtual void cal_force(UnitCell& ucell, ModuleBase::matrix& force) override {}; + virtual void cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) override {}; - protected: + protected: const Input_para& input; const UnitCell& ucell; std::vector orb_cutoff_; diff --git a/source/module_md/md_func.cpp b/source/module_md/md_func.cpp index 9387482c50..3d40ca3f7e 100644 --- a/source/module_md/md_func.cpp +++ b/source/module_md/md_func.cpp @@ -255,16 +255,16 @@ void force_virial(ModuleESolver::ESolver* p_esolver, ModuleBase::TITLE("MD_func", "force_virial"); ModuleBase::timer::tick("MD_func", "force_virial"); - p_esolver->runner(istep, unit_in); + p_esolver->runner(unit_in, istep); potential = p_esolver->cal_energy(); ModuleBase::matrix force_temp(unit_in.nat, 3); - p_esolver->cal_force(force_temp); + p_esolver->cal_force(unit_in, force_temp); if (cal_stress) { - p_esolver->cal_stress(virial); + p_esolver->cal_stress(unit_in, virial); } /// convert Rydberg to Hartree diff --git a/source/module_md/test/fire_test.cpp b/source/module_md/test/fire_test.cpp index 2ecdf3030b..81573498ae 100644 --- a/source/module_md/test/fire_test.cpp +++ b/source/module_md/test/fire_test.cpp @@ -49,7 +49,7 @@ class FIREtest : public testing::Test Setcell::parameters(param_in.input); p_esolver = new ModuleESolver::ESolver_LJ(); - p_esolver->before_all_runners(param_in.inp, ucell); + p_esolver->before_all_runners(ucell, param_in.inp); mdrun = new FIRE(param_in, ucell); mdrun->setup(p_esolver, PARAM.sys.global_readin_dir); diff --git a/source/module_md/test/langevin_test.cpp b/source/module_md/test/langevin_test.cpp index 58b9af787c..a06f279c4f 100644 --- a/source/module_md/test/langevin_test.cpp +++ b/source/module_md/test/langevin_test.cpp @@ -49,7 +49,7 @@ class Langevin_test : public testing::Test Setcell::parameters(param_in.input); p_esolver = new ModuleESolver::ESolver_LJ(); - p_esolver->before_all_runners(param_in.inp, ucell); + p_esolver->before_all_runners(ucell, param_in.inp); mdrun = new Langevin(param_in, ucell); mdrun->setup(p_esolver, PARAM.sys.global_readin_dir); diff --git a/source/module_md/test/lj_pot_test.cpp b/source/module_md/test/lj_pot_test.cpp index 8bc03b6762..84d91c3f87 100644 --- a/source/module_md/test/lj_pot_test.cpp +++ b/source/module_md/test/lj_pot_test.cpp @@ -47,7 +47,7 @@ class LJ_pot_test : public testing::Test TEST_F(LJ_pot_test, potential) { ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); - p_esolver->before_all_runners(input, ucell); + p_esolver->before_all_runners(ucell, input); MD_func::force_virial(p_esolver, 0, ucell, potential, force, true, stress); EXPECT_NEAR(potential, -0.011957818623534381, doublethreshold); } @@ -55,7 +55,7 @@ TEST_F(LJ_pot_test, potential) TEST_F(LJ_pot_test, force) { ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); - p_esolver->before_all_runners(input, ucell); + p_esolver->before_all_runners(ucell, input); MD_func::force_virial(p_esolver, 0, ucell, potential, force, true, stress); EXPECT_NEAR(force[0].x, 0.00049817733089377704, doublethreshold); EXPECT_NEAR(force[0].y, 0.00082237246837022328, doublethreshold); @@ -74,7 +74,7 @@ TEST_F(LJ_pot_test, force) TEST_F(LJ_pot_test, stress) { ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); - p_esolver->before_all_runners(input, ucell); + p_esolver->before_all_runners(ucell, input); MD_func::force_virial(p_esolver, 0, ucell, potential, force, true, stress); EXPECT_NEAR(stress(0, 0), 8.0360222227631859e-07, doublethreshold); EXPECT_NEAR(stress(0, 1), 1.7207745586539077e-07, doublethreshold); diff --git a/source/module_md/test/msst_test.cpp b/source/module_md/test/msst_test.cpp index 19441d52c3..27b94e5a2d 100644 --- a/source/module_md/test/msst_test.cpp +++ b/source/module_md/test/msst_test.cpp @@ -49,7 +49,7 @@ class MSST_test : public testing::Test Setcell::parameters(param_in.input); p_esolver = new ModuleESolver::ESolver_LJ(); - p_esolver->before_all_runners(param_in.inp, ucell); + p_esolver->before_all_runners(ucell, param_in.inp); mdrun = new MSST(param_in, ucell); mdrun->setup(p_esolver, PARAM.sys.global_readin_dir); diff --git a/source/module_md/test/nhchain_test.cpp b/source/module_md/test/nhchain_test.cpp index 1370cd81f5..22f0bbf9a0 100644 --- a/source/module_md/test/nhchain_test.cpp +++ b/source/module_md/test/nhchain_test.cpp @@ -47,7 +47,7 @@ class NHC_test : public testing::Test Setcell::parameters(param_in.input); p_esolver = new ModuleESolver::ESolver_LJ(); - p_esolver->before_all_runners(param_in.inp, ucell); + p_esolver->before_all_runners(ucell, param_in.inp); param_in.input.mdp.md_type = "npt"; param_in.input.mdp.md_pmode = "tri"; diff --git a/source/module_md/test/verlet_test.cpp b/source/module_md/test/verlet_test.cpp index 378c390563..578761dcc5 100644 --- a/source/module_md/test/verlet_test.cpp +++ b/source/module_md/test/verlet_test.cpp @@ -50,7 +50,7 @@ class Verlet_test : public testing::Test Setcell::parameters(param_in.input); p_esolver = new ModuleESolver::ESolver_LJ(); - p_esolver->before_all_runners(param_in.inp, ucell); + p_esolver->before_all_runners(ucell, param_in.inp); mdrun = new Verlet(param_in, ucell); mdrun->setup(p_esolver, PARAM.sys.global_readin_dir); diff --git a/source/module_relax/relax_driver.cpp b/source/module_relax/relax_driver.cpp index 9806ea975c..9a5820d92a 100644 --- a/source/module_relax/relax_driver.cpp +++ b/source/module_relax/relax_driver.cpp @@ -9,7 +9,7 @@ #include "module_io/read_exit_file.h" #include "module_io/write_wfc_r.h" #include "module_parameter/parameter.h" -void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver) +void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver, UnitCell& ucell) { ModuleBase::TITLE("Ions", "opt_ions"); ModuleBase::timer::tick("Ions", "opt_ions"); @@ -18,11 +18,11 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver) { if (!PARAM.inp.relax_new) { - rl_old.init_relax(GlobalC::ucell.nat); + rl_old.init_relax(ucell.nat); } else { - rl.init_relax(GlobalC::ucell.nat); + rl.init_relax(ucell.nat); } } @@ -48,7 +48,7 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver) #endif //__RAPIDJSON // mohan added eiter to count for the electron iteration number, 2021-01-28 - p_esolver->runner(istep - 1, GlobalC::ucell); + p_esolver->runner(ucell, istep - 1); time_t eend = time(nullptr); time_t fstart = time(nullptr); @@ -67,12 +67,12 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver) // calculate and gather all parts of total ionic forces if (PARAM.inp.cal_force) { - p_esolver->cal_force(force); + p_esolver->cal_force(ucell, force); } // calculate and gather all parts of stress if (PARAM.inp.cal_stress) { - p_esolver->cal_stress(stress); + p_esolver->cal_stress(ucell, stress); } if (PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax") @@ -85,7 +85,7 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver) { stop = rl_old.relax_step(istep, this->etot, - GlobalC::ucell, + ucell, force, stress, force_step, @@ -102,29 +102,29 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver) need_orb = need_orb || PARAM.inp.basis_type == "lcao_in_pw"; std::stringstream ss, ss1; ss << PARAM.globalv.global_out_dir << "STRU_ION_D"; - GlobalC::ucell.print_stru_file(ss.str(), - PARAM.inp.nspin, - true, - PARAM.inp.calculation == "md", - PARAM.inp.out_mul, - need_orb, - PARAM.globalv.deepks_setorb, - GlobalV::MY_RANK); + ucell.print_stru_file(ss.str(), + PARAM.inp.nspin, + true, + PARAM.inp.calculation == "md", + PARAM.inp.out_mul, + need_orb, + PARAM.globalv.deepks_setorb, + GlobalV::MY_RANK); if (Ions_Move_Basic::out_stru) { ss1 << PARAM.globalv.global_out_dir << "STRU_ION"; ss1 << istep << "_D"; - GlobalC::ucell.print_stru_file(ss1.str(), - PARAM.inp.nspin, - true, - PARAM.inp.calculation == "md", - PARAM.inp.out_mul, - need_orb, - PARAM.globalv.deepks_setorb, - GlobalV::MY_RANK); + ucell.print_stru_file(ss1.str(), + PARAM.inp.nspin, + true, + PARAM.inp.calculation == "md", + PARAM.inp.out_mul, + need_orb, + PARAM.globalv.deepks_setorb, + GlobalV::MY_RANK); ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU_NOW.cif", - GlobalC::ucell, + ucell, "# Generated by ABACUS ModuleIO::CifParser", "data_?"); } @@ -141,7 +141,7 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver) // add Json of cell coo stress force double unit_transform = ModuleBase::RYDBERG_SI / pow(ModuleBase::BOHR_RADIUS_SI, 3) * 1.0e-8; double fac = ModuleBase::Ry_to_eV / 0.529177; - Json::add_output_cell_coo_stress_force(&GlobalC::ucell, force, fac, stress, unit_transform); + Json::add_output_cell_coo_stress_force(&ucell, force, fac, stress, unit_transform); #endif //__RAPIDJSON if (stop == false) diff --git a/source/module_relax/relax_driver.h b/source/module_relax/relax_driver.h index 926f66b7af..84a5a566e0 100644 --- a/source/module_relax/relax_driver.h +++ b/source/module_relax/relax_driver.h @@ -1,6 +1,7 @@ #ifndef RELAX_DRIVER_H #define RELAX_DRIVER_H +#include "module_cell/unitcell.h" #include "module_esolver/esolver.h" #include "module_esolver/esolver_ks.h" #include "relax_new/relax.h" @@ -13,7 +14,7 @@ class Relax_Driver Relax_Driver(){}; ~Relax_Driver(){}; - void relax_driver(ModuleESolver::ESolver *p_esolver); + void relax_driver(ModuleESolver::ESolver* p_esolver, UnitCell& ucell); private: // mohan add 2021-01-28