From f5ca1248ce8eecd110d1d0dfd72c8ca0d47290c3 Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Fri, 26 Jul 2024 20:35:21 +0800 Subject: [PATCH] Refactor: update md parameters to apply const param_in --- source/driver_run.cpp | 14 ++- .../hamilt_lcaodft/LCAO_set_st.cpp | 4 +- .../module_tddft/propagator.cpp | 18 ++-- source/module_io/input.h | 4 +- source/module_io/input_conv.cpp | 4 +- source/module_io/input_conv_tmp.cpp | 1 - source/module_io/print_info.cpp | 26 ++--- source/module_io/read_input_item_md.cpp | 12 +++ source/module_io/read_input_item_relax.cpp | 4 - source/module_md/fire.cpp | 25 ++--- source/module_md/fire.h | 4 +- source/module_md/langevin.cpp | 12 +-- source/module_md/langevin.h | 3 +- source/module_md/md_base.cpp | 38 ++++--- source/module_md/md_base.h | 12 ++- source/module_md/md_func.cpp | 21 ++-- source/module_md/md_func.h | 4 +- source/module_md/msst.cpp | 45 +++++---- source/module_md/msst.h | 5 +- source/module_md/nhchain.cpp | 99 ++++++++++--------- source/module_md/nhchain.h | 6 +- source/module_md/run_md.cpp | 34 +++---- source/module_md/run_md.h | 4 +- source/module_md/test/md_func_test.cpp | 4 +- source/module_md/test/setcell.h | 2 +- source/module_md/verlet.cpp | 12 +-- source/module_md/verlet.h | 2 +- source/module_parameter/md_parameter.h | 6 -- source/module_parameter/parameter.cpp | 1 - 29 files changed, 220 insertions(+), 206 deletions(-) diff --git a/source/driver_run.cpp b/source/driver_run.cpp index fd6362f3c2..68dfab2539 100644 --- a/source/driver_run.cpp +++ b/source/driver_run.cpp @@ -59,13 +59,17 @@ void Driver::driver_run() { const std::string cal_type = GlobalV::CALCULATION; //! 4: different types of calculations - if (cal_type == "md") { - Run_MD::md_line(GlobalC::ucell, p_esolver, INPUT.mdp); - } else if (cal_type == "scf" || cal_type == "relax" - || cal_type == "cell-relax") { + if (cal_type == "md") + { + Run_MD::md_line(GlobalC::ucell, p_esolver, PARAM); + } + else if (cal_type == "scf" || cal_type == "relax" || cal_type == "cell-relax") + { Relax_Driver rl_driver; rl_driver.relax_driver(p_esolver); - } else { + } + else + { //! supported "other" functions: //! nscf(PW,LCAO), //! get_pchg(LCAO), diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_set_st.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_set_st.cpp index b877f23098..0983550127 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_set_st.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_set_st.cpp @@ -1,6 +1,6 @@ #include "module_base/timer.h" #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" // only for INPUT +#include "module_parameter/parameter.h" namespace LCAO_domain { @@ -362,7 +362,7 @@ void build_ST_new(ForceStressArrays& fsr, { for (int k = 0; k < 3; k++) { - tau1[k] = tau1[k] - atom1->vel[I1][k] * INPUT.mdp.md_dt / ucell.lat0; + tau1[k] = tau1[k] - atom1->vel[I1][k] * PARAM.mdp.md_dt / ucell.lat0; } } diff --git a/source/module_hamilt_lcao/module_tddft/propagator.cpp b/source/module_hamilt_lcao/module_tddft/propagator.cpp index 8ac07aa2a6..1371780d11 100644 --- a/source/module_hamilt_lcao/module_tddft/propagator.cpp +++ b/source/module_hamilt_lcao/module_tddft/propagator.cpp @@ -1,11 +1,11 @@ #include "propagator.h" -#include -#include - #include "module_base/lapack_connector.h" #include "module_base/scalapack_connector.h" -#include "module_io/input.h" +#include "module_parameter/parameter.h" + +#include +#include namespace module_tddft { @@ -97,11 +97,11 @@ void Propagator::compute_propagator_cn2(const int nlocal, // ->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> // (2) compute Numerator & Denominator by GEADD - // Numerator = Stmp - i*para * Htmp; beta1 = - para = -0.25 * INPUT.mdp.md_dt - // Denominator = Stmp + i*para * Htmp; beta2 = para = 0.25 * INPUT.mdp.md_dt + // Numerator = Stmp - i*para * Htmp; beta1 = - para = -0.25 * PARAM.mdp.md_dt + // Denominator = Stmp + i*para * Htmp; beta2 = para = 0.25 * PARAM.mdp.md_dt std::complex alpha = {1.0, 0.0}; - std::complex beta1 = {0.0, -0.25 * INPUT.mdp.md_dt}; - std::complex beta2 = {0.0, 0.25 * INPUT.mdp.md_dt}; + std::complex beta1 = {0.0, -0.25 * PARAM.mdp.md_dt}; + std::complex beta2 = {0.0, 0.25 * PARAM.mdp.md_dt}; ScalapackConnector::geadd('N', nlocal, @@ -350,7 +350,7 @@ void Propagator::compute_propagator_taylor(const int nlocal, } // loop ipcol } // loop iprow - std::complex beta = {0.0, -0.5 * INPUT.mdp.md_dt / tag}; // for ETRS tag=2 , for taylor tag=1 + std::complex beta = {0.0, -0.5 * PARAM.mdp.md_dt / tag}; // for ETRS tag=2 , for taylor tag=1 //->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> // invert Stmp diff --git a/source/module_io/input.h b/source/module_io/input.h index 99d4393775..c4b7eb6669 100644 --- a/source/module_io/input.h +++ b/source/module_io/input.h @@ -19,10 +19,8 @@ class Input } // They will be removed. - int cond_dtbatch; + int cond_dtbatch; int nche_sto; - double md_tfirst; - MD_para mdp; int* orbital_corr = nullptr; ///< which correlated orbitals need corrected ; double* hubbard_u = nullptr; ///< Hubbard Coulomb interaction parameter U(ev) std::string stru_file; // file contains atomic positions -- xiaohui modify diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index 4d2d2cdb02..471cae2007 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -183,8 +183,8 @@ void Input_Conv::Convert() if (PARAM.inp.calculation == "md" && PARAM.mdp.md_restart) // md restart liuyu add 2023-04-12 { int istep = 0; - MD_func::current_md_info(GlobalV::MY_RANK, GlobalV::global_readin_dir, istep, INPUT.md_tfirst); - INPUT.md_tfirst *= ModuleBase::Hartree_to_K; + double temperature = 0.0; + MD_func::current_md_info(GlobalV::MY_RANK, GlobalV::global_readin_dir, istep, temperature); if (PARAM.inp.read_file_dir == "auto") { GlobalV::stru_file = INPUT.stru_file = GlobalV::global_stru_dir + "STRU_MD_" + std::to_string(istep); diff --git a/source/module_io/input_conv_tmp.cpp b/source/module_io/input_conv_tmp.cpp index af3cd1cb0e..7049c91d71 100644 --- a/source/module_io/input_conv_tmp.cpp +++ b/source/module_io/input_conv_tmp.cpp @@ -10,7 +10,6 @@ void Input_Conv::tmp_convert() INPUT.stru_file = PARAM.inp.stru_file; INPUT.cond_dtbatch = PARAM.inp.cond_dtbatch; INPUT.nche_sto = PARAM.inp.nche_sto; - INPUT.mdp = PARAM.mdp; const int ntype = PARAM.inp.ntype; delete[] INPUT.orbital_corr; diff --git a/source/module_io/print_info.cpp b/source/module_io/print_info.cpp index e98324b013..bc7552f8b2 100644 --- a/source/module_io/print_info.cpp +++ b/source/module_io/print_info.cpp @@ -1,7 +1,7 @@ #include "print_info.h" -#include "module_io/input.h" -#include "../module_base/global_variable.h" -//#include "../module_cell/klist.h" + +#include "module_base/global_variable.h" +#include "module_parameter/parameter.h" Print_Info::Print_Info(){} @@ -38,32 +38,34 @@ void Print_Info::setup_parameters(UnitCell &ucell, K_Vectors &kv) std::cout << " ---------------------------------------------------------" << std::endl; - if(INPUT.mdp.md_type == "fire") + if (PARAM.mdp.md_type == "fire") { std::cout << " ENSEMBLE : " << "FIRE" << std::endl; } - else if(INPUT.mdp.md_type == "nve") + else if (PARAM.mdp.md_type == "nve") { std::cout << " ENSEMBLE : " << "NVE" << std::endl; } - else if(INPUT.mdp.md_type == "nvt") + else if (PARAM.mdp.md_type == "nvt") { - std::cout << " ENSEMBLE : " << "NVT mode: " << INPUT.mdp.md_thermostat << std::endl; + std::cout << " ENSEMBLE : " + << "NVT mode: " << PARAM.mdp.md_thermostat << std::endl; } - else if(INPUT.mdp.md_type == "npt") + else if (PARAM.mdp.md_type == "npt") { - std::cout << " ENSEMBLE : " << "NPT mode: " << INPUT.mdp.md_pmode << std::endl; + std::cout << " ENSEMBLE : " + << "NPT mode: " << PARAM.mdp.md_pmode << std::endl; } - else if(INPUT.mdp.md_type == "langevin") + else if (PARAM.mdp.md_type == "langevin") { std::cout << " ENSEMBLE : " << "Langevin" << std::endl; } - else if(INPUT.mdp.md_type == "msst") + else if (PARAM.mdp.md_type == "msst") { std::cout << " ENSEMBLE : " << "MSST" << std::endl; } - std::cout << " Time interval(fs) : " << INPUT.mdp.md_dt << std::endl; + std::cout << " Time interval(fs) : " << PARAM.mdp.md_dt << std::endl; } std::cout << " ---------------------------------------------------------" << std::endl; diff --git a/source/module_io/read_input_item_md.cpp b/source/module_io/read_input_item_md.cpp index 48a311c952..98e284ecc5 100644 --- a/source/module_io/read_input_item_md.cpp +++ b/source/module_io/read_input_item_md.cpp @@ -58,6 +58,12 @@ void ReadInput::item_md() { Input_Item item("md_tlast"); item.annotation = "temperature last"; + item.reset_value = [](const Input_Item& item, Parameter& para) { + if (para.mdp.md_tlast < 0) + { + para.input.mdp.md_tlast = para.mdp.md_tfirst; + } + }; read_sync_double(input.mdp.md_tlast); this->add_item(item); } @@ -296,6 +302,12 @@ void ReadInput::item_md() { Input_Item item("md_pcouple"); item.annotation = "whether couple different components: xyz, xy, yz, xz, none"; + item.reset_value = [](const Input_Item& item, Parameter& para) { + if (para.mdp.md_pmode == "iso") + { + para.input.mdp.md_pcouple = "xyz"; + } + }; read_sync_string(input.mdp.md_pcouple); this->add_item(item); } diff --git a/source/module_io/read_input_item_relax.cpp b/source/module_io/read_input_item_relax.cpp index ede0587189..8bb9834eba 100644 --- a/source/module_io/read_input_item_relax.cpp +++ b/source/module_io/read_input_item_relax.cpp @@ -200,10 +200,8 @@ void ReadInput::item_relax() ModuleBase::WARNING("ReadInput", "both force_thr and force_thr_ev are set, use force_thr"); para.input.force_thr_ev = para.input.force_thr * 13.6058 / 0.529177; } - para.input.mdp.force_thr = para.input.force_thr; // temperaory }; sync_double(input.force_thr); - add_double_bcast(input.mdp.force_thr); this->add_item(item); } { @@ -295,10 +293,8 @@ void ReadInput::item_relax() { para.input.cal_stress = true; } - para.input.mdp.cal_stress = para.input.cal_stress; // temperaory }; read_sync_bool(input.cal_stress); - add_bool_bcast(input.mdp.cal_stress); this->add_item(item); } { diff --git a/source/module_md/fire.cpp b/source/module_md/fire.cpp index 3518f71c6a..3299e05795 100644 --- a/source/module_md/fire.cpp +++ b/source/module_md/fire.cpp @@ -6,8 +6,9 @@ #endif #include "module_base/timer.h" -FIRE::FIRE(MD_para& MD_para_in, UnitCell& unit_in) : MD_base(MD_para_in, unit_in) +FIRE::FIRE(const Parameter& param_in, UnitCell& unit_in) : MD_base(param_in, unit_in) { + force_thr = param_in.inp.force_thr; dt_max = -1.0; alpha_start = 0.10; alpha = alpha_start; @@ -82,18 +83,18 @@ void FIRE::print_md(std::ofstream& ofs, const bool& cal_stress) void FIRE::write_restart(const std::string& global_out_dir) { - if (!mdp.my_rank) + if (!my_rank) { std::stringstream ssc; ssc << global_out_dir << "Restart_md.dat"; std::ofstream file(ssc.str().c_str()); file << step_ + step_rst_ << std::endl; - file << mdp.md_tfirst << std::endl; + file << md_tfirst << std::endl; file << alpha << std::endl; file << negative_count << std::endl; file << dt_max << std::endl; - file << mdp.md_dt << std::endl; + file << md_dt << std::endl; file.close(); } #ifdef __MPI @@ -108,7 +109,7 @@ void FIRE::restart(const std::string& global_readin_dir) { bool ok = true; - if (!mdp.my_rank) + if (!my_rank) { std::stringstream ssc; ssc << global_readin_dir << "Restart_md.dat"; @@ -121,7 +122,7 @@ void FIRE::restart(const std::string& global_readin_dir) if (ok) { - file >> step_rst_ >> mdp.md_tfirst >> alpha >> negative_count >> dt_max >> mdp.md_dt; + file >> step_rst_ >> md_tfirst >> alpha >> negative_count >> dt_max >> md_dt; file.close(); } } @@ -137,11 +138,11 @@ void FIRE::restart(const std::string& global_readin_dir) #ifdef __MPI MPI_Bcast(&step_rst_, 1, MPI_INT, 0, MPI_COMM_WORLD); - MPI_Bcast(&mdp.md_tfirst, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&md_tfirst, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(&alpha, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(&negative_count, 1, MPI_INT, 0, MPI_COMM_WORLD); MPI_Bcast(&dt_max, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast(&mdp.md_dt, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&md_dt, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); #endif return; @@ -163,7 +164,7 @@ void FIRE::check_force(void) } } - if (2.0 * max < mdp.force_thr) + if (2.0 * max < force_thr) { stop = true; } @@ -181,7 +182,7 @@ void FIRE::check_fire(void) /// initial dt_max if (dt_max < 0) { - dt_max = 2.5 * mdp.md_dt; + dt_max = 2.5 * md_dt; } for (int i = 0; i < ucell.nat; ++i) @@ -207,13 +208,13 @@ void FIRE::check_fire(void) negative_count++; if (negative_count >= n_min) { - mdp.md_dt = std::min(mdp.md_dt * finc, dt_max); + md_dt = std::min(md_dt * finc, dt_max); alpha *= f_alpha; } } else { - mdp.md_dt *= fdec; + md_dt *= fdec; negative_count = 0; for (int i = 0; i < ucell.nat; ++i) diff --git a/source/module_md/fire.h b/source/module_md/fire.h index 7f506fbd02..2528e44b20 100644 --- a/source/module_md/fire.h +++ b/source/module_md/fire.h @@ -13,8 +13,7 @@ class FIRE : public MD_base { public: - - FIRE(MD_para& MD_para_in, UnitCell& unit_in); + FIRE(const Parameter& param_in, UnitCell& unit_in); ~FIRE(); @@ -51,6 +50,7 @@ class FIRE : public MD_base int n_min; ///< n_min double dt_max; ///< dt_max int negative_count; ///< Negative count + double force_thr = 1.0e-3; ///< force convergence threshold in FIRE method }; #endif diff --git a/source/module_md/langevin.cpp b/source/module_md/langevin.cpp index 38e2ad16fe..16ba53041a 100644 --- a/source/module_md/langevin.cpp +++ b/source/module_md/langevin.cpp @@ -4,12 +4,12 @@ #include "module_base/parallel_common.h" #include "module_base/timer.h" -Langevin::Langevin(MD_para& MD_para_in, UnitCell& unit_in) : MD_base(MD_para_in, unit_in) +Langevin::Langevin(const Parameter& param_in, UnitCell& unit_in) : MD_base(param_in, unit_in) { /// convert to a.u. unit assert(ModuleBase::AU_to_FS!=0.0); - mdp.md_damp /= ModuleBase::AU_to_FS; + md_damp = mdp.md_damp / ModuleBase::AU_to_FS; assert(ucell.nat>0); @@ -85,16 +85,16 @@ void Langevin::restart(const std::string& global_readin_dir) void Langevin::post_force(void) { - if (mdp.my_rank == 0) + if (my_rank == 0) { - double t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, mdp.md_tfirst, mdp.md_tlast); + double t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, md_tfirst, md_tlast); ModuleBase::Vector3 fictitious_force; for (int i = 0; i < ucell.nat; ++i) { - fictitious_force = -allmass[i] * vel[i] / mdp.md_damp; + fictitious_force = -allmass[i] * vel[i] / md_damp; for (int j = 0; j < 3; ++j) { - fictitious_force[j] += sqrt(24.0 * t_target * allmass[i] / mdp.md_damp / mdp.md_dt) + fictitious_force[j] += sqrt(24.0 * t_target * allmass[i] / md_damp / md_dt) * (static_cast(std::rand()) / RAND_MAX - 0.5); } total_force[i] = force[i] + fictitious_force; diff --git a/source/module_md/langevin.h b/source/module_md/langevin.h index 7f3645ba93..e89122d9f1 100644 --- a/source/module_md/langevin.h +++ b/source/module_md/langevin.h @@ -14,7 +14,7 @@ class Langevin : public MD_base { public: - Langevin(MD_para& MD_para_in, UnitCell& unit_in); + Langevin(const Parameter& param_in, UnitCell& unit_in); ~Langevin(); @@ -38,6 +38,7 @@ class Langevin : public MD_base void post_force(); ModuleBase::Vector3* total_force; ///< total force = true force + Langevin fictitious_force + double md_damp; ///< damping factor }; #endif diff --git a/source/module_md/md_base.cpp b/source/module_md/md_base.cpp index 665dfec20f..17a88e64ad 100644 --- a/source/module_md/md_base.cpp +++ b/source/module_md/md_base.cpp @@ -6,9 +6,10 @@ #endif #include "module_io/print_info.h" -MD_base::MD_base(MD_para& MD_para_in, UnitCell& unit_in) -: mdp(MD_para_in), ucell(unit_in) +MD_base::MD_base(const Parameter& param_in, UnitCell& unit_in) : mdp(param_in.mdp), ucell(unit_in) { + my_rank = param_in.globalv.myrank; + cal_stress = param_in.inp.cal_stress; if (mdp.md_seed >= 0) { srand(mdp.md_seed); @@ -30,20 +31,15 @@ MD_base::MD_base(MD_para& MD_para_in, UnitCell& unit_in) assert(ModuleBase::Hartree_to_K!=0.0); /// convert to a.u. unit - mdp.md_dt /= ModuleBase::AU_to_FS; - mdp.md_tfirst /= ModuleBase::Hartree_to_K; - mdp.md_tlast /= ModuleBase::Hartree_to_K; - mdp.md_tfreq *= ModuleBase::AU_to_FS; + md_dt = mdp.md_dt / ModuleBase::AU_to_FS; + md_tfirst = mdp.md_tfirst / ModuleBase::Hartree_to_K; + md_tlast = mdp.md_tlast / ModuleBase::Hartree_to_K; step_ = 0; step_rst_ = 0; - MD_func::init_vel(ucell, mdp.my_rank, mdp.md_restart, mdp.md_tfirst, allmass, frozen_freedom_, ionmbl, vel); + MD_func::init_vel(ucell, my_rank, mdp.md_restart, md_tfirst, allmass, frozen_freedom_, ionmbl, vel); t_current = MD_func::current_temp(kinetic, ucell.nat, frozen_freedom_, allmass, vel); - if (mdp.md_tlast < 0) - { - mdp.md_tlast = mdp.md_tfirst; - } } @@ -66,8 +62,8 @@ void MD_base::setup(ModuleESolver::ESolver* p_esolver, const std::string& global Print_Info::print_screen(0, 0, step_ + step_rst_); - MD_func::force_virial(p_esolver, step_, ucell, potential, force, mdp.cal_stress, virial); - MD_func::compute_stress(ucell, vel, allmass, mdp.cal_stress, virial, stress); + MD_func::force_virial(p_esolver, step_, ucell, potential, force, cal_stress, virial); + MD_func::compute_stress(ucell, vel, allmass, cal_stress, virial, stress); ucell.ionic_position_updated = true; return; @@ -93,7 +89,7 @@ void MD_base::second_half(void) void MD_base::update_pos() { - if (mdp.my_rank == 0) + if (my_rank == 0) { for (int i = 0; i < ucell.nat; ++i) { @@ -101,7 +97,7 @@ void MD_base::update_pos() { if (ionmbl[i][k]) { - pos[i][k] = vel[i][k] * mdp.md_dt / ucell.lat0; + pos[i][k] = vel[i][k] * md_dt / ucell.lat0; } else { @@ -124,7 +120,7 @@ void MD_base::update_pos() void MD_base::update_vel(const ModuleBase::Vector3* force) { - if (mdp.my_rank == 0) + if (my_rank == 0) { for (int i = 0; i < ucell.nat; ++i) { @@ -132,7 +128,7 @@ void MD_base::update_vel(const ModuleBase::Vector3* force) { if (ionmbl[i][k]) { - vel[i][k] += 0.5 * force[i][k] * mdp.md_dt / allmass[i]; + vel[i][k] += 0.5 * force[i][k] * md_dt / allmass[i]; } } } @@ -147,7 +143,7 @@ void MD_base::update_vel(const ModuleBase::Vector3* force) void MD_base::print_md(std::ofstream& ofs, const bool& cal_stress) { - if (mdp.my_rank) + if (my_rank) { return; } @@ -228,14 +224,14 @@ void MD_base::print_md(std::ofstream& ofs, const bool& cal_stress) void MD_base::write_restart(const std::string& global_out_dir) { - if (!mdp.my_rank) + if (!my_rank) { std::stringstream ssc; ssc << global_out_dir << "Restart_md.dat"; std::ofstream file(ssc.str().c_str()); file << step_ + step_rst_ << std::endl; - file << mdp.md_tfirst << std::endl; + file << md_tfirst << std::endl; file.close(); } #ifdef __MPI @@ -248,7 +244,7 @@ void MD_base::write_restart(const std::string& global_out_dir) void MD_base::restart(const std::string& global_readin_dir) { - MD_func::current_md_info(mdp.my_rank, global_readin_dir, step_rst_, mdp.md_tfirst); + MD_func::current_md_info(my_rank, global_readin_dir, step_rst_, md_tfirst); return; } diff --git a/source/module_md/md_base.h b/source/module_md/md_base.h index ed8c6467cc..ca6eb0dcf7 100644 --- a/source/module_md/md_base.h +++ b/source/module_md/md_base.h @@ -2,7 +2,7 @@ #define MD_BASE_H #include "module_esolver/esolver.h" -#include "module_parameter/md_parameter.h" +#include "module_parameter/parameter.h" /** * @brief base class of md @@ -15,7 +15,7 @@ class MD_base { public: - MD_base(MD_para& MD_para_in, UnitCell& unit_in); + MD_base(const Parameter& param_in, UnitCell& unit_in); virtual ~MD_base(); /** @@ -86,9 +86,15 @@ class MD_base double kinetic; ///< kinetic energy protected: - MD_para& mdp; ///< input parameters used in md + const MD_para& mdp; ///< input parameters used in md UnitCell& ucell; ///< unitcell information double energy_; ///< total energy of the system + + bool cal_stress; ///< whether calculate stress + int my_rank; ///< MPI rank of the processor + double md_dt; ///< Time increment (hbar/E_hartree) + double md_tfirst; ///< Temperature (in Hartree, 1 Hartree ~ 3E5 K) + double md_tlast; ///< Target temperature }; #endif // MD_BASE_H \ No newline at end of file diff --git a/source/module_md/md_func.cpp b/source/module_md/md_func.cpp index e362bb361a..da6131d198 100644 --- a/source/module_md/md_func.cpp +++ b/source/module_md/md_func.cpp @@ -302,19 +302,18 @@ void print_stress(std::ofstream& ofs, const ModuleBase::matrix& virial, const Mo return; } - void dump_info(const int& step, const std::string& global_out_dir, const UnitCell& unit_in, - const MD_para& mdp, + const Parameter& param_in, const ModuleBase::matrix& virial, const ModuleBase::Vector3* force, const ModuleBase::Vector3* vel) { - if (mdp.my_rank) - { - return; - } + if (param_in.globalv.myrank) + { + return; + } std::stringstream file; file << global_out_dir << "MD_dump"; @@ -343,7 +342,7 @@ void dump_info(const int& step, ofs << " " << unit_in.latvec.e21 << " " << unit_in.latvec.e22 << " " << unit_in.latvec.e23 << std::endl; ofs << " " << unit_in.latvec.e31 << " " << unit_in.latvec.e32 << " " << unit_in.latvec.e33 << std::endl; - if (mdp.cal_stress && mdp.dump_virial) + if (param_in.inp.cal_stress && param_in.mdp.dump_virial) { ofs << "VIRIAL (kbar)" << std::endl; for (int i = 0; i < 3; ++i) @@ -354,11 +353,11 @@ void dump_info(const int& step, } ofs << "INDEX LABEL POSITION (Angstrom)"; - if (mdp.dump_force) + if (param_in.mdp.dump_force) { ofs << " FORCE (eV/Angstrom)"; } - if (mdp.dump_vel) + if (param_in.mdp.dump_vel) { ofs << " VELOCITY (Angstrom/fs)"; } @@ -372,13 +371,13 @@ void dump_info(const int& step, ofs << " " << index << " " << unit_in.atom_label[it] << " " << unit_in.atoms[it].tau[ia].x * unit_pos << " " << unit_in.atoms[it].tau[ia].y * unit_pos << " " << unit_in.atoms[it].tau[ia].z * unit_pos; - if (mdp.dump_force) + if (param_in.mdp.dump_force) { ofs << " " << force[index].x * unit_force << " " << force[index].y * unit_force << " " << force[index].z * unit_force; } - if (mdp.dump_vel) + if (param_in.mdp.dump_vel) { ofs << " " << vel[index].x * unit_vel << " " << vel[index].y * unit_vel << " " << vel[index].z * unit_vel; diff --git a/source/module_md/md_func.h b/source/module_md/md_func.h index f2b93a9d56..7d91e05236 100644 --- a/source/module_md/md_func.h +++ b/source/module_md/md_func.h @@ -131,7 +131,7 @@ void print_stress(std::ofstream& ofs, const ModuleBase::matrix& virial, const Mo * @param step current md step * @param global_out_dir directory of output files * @param unit_in unitcell information - * @param mdp input parameters used in md + * @param param_in input parameters used in md * @param virial lattice virial tensor * @param force atomic forces * @param vel atomic velocities @@ -139,7 +139,7 @@ void print_stress(std::ofstream& ofs, const ModuleBase::matrix& virial, const Mo void dump_info(const int& step, const std::string& global_out_dir, const UnitCell& unit_in, - const MD_para& mdp, + const Parameter& param_in, const ModuleBase::matrix& virial, const ModuleBase::Vector3* force, const ModuleBase::Vector3* vel); diff --git a/source/module_md/msst.cpp b/source/module_md/msst.cpp index 1a7e3c7aef..67484f0f88 100644 --- a/source/module_md/msst.cpp +++ b/source/module_md/msst.cpp @@ -6,12 +6,11 @@ #endif #include "module_base/timer.h" -MSST::MSST(MD_para& MD_para_in, UnitCell& unit_in) -: MD_base(MD_para_in, unit_in) +MSST::MSST(const Parameter& param_in, UnitCell& unit_in) : MD_base(param_in, unit_in) { - mdp.msst_qmass = mdp.msst_qmass / pow(ModuleBase::ANGSTROM_AU, 4) / pow(ModuleBase::AU_to_MASS, 2); - mdp.msst_vel = mdp.msst_vel * ModuleBase::ANGSTROM_AU * ModuleBase::AU_to_FS; - mdp.msst_vis = mdp.msst_vis / ModuleBase::AU_to_MASS / ModuleBase::ANGSTROM_AU * ModuleBase::AU_to_FS; + msst_qmass = mdp.msst_qmass / pow(ModuleBase::ANGSTROM_AU, 4) / pow(ModuleBase::AU_to_MASS, 2); + msst_vel = mdp.msst_vel * ModuleBase::ANGSTROM_AU * ModuleBase::AU_to_FS; + msst_vis = mdp.msst_vis / ModuleBase::AU_to_MASS / ModuleBase::ANGSTROM_AU * ModuleBase::AU_to_FS; assert(ucell.nat>0); @@ -53,7 +52,7 @@ void MSST::setup(ModuleESolver::ESolver* p_esolver, const std::string& global_re if (kinetic > 0 && mdp.msst_tscale > 0) { - double fac1 = mdp.msst_tscale * totmass * 2.0 * kinetic / mdp.msst_qmass; + double fac1 = mdp.msst_tscale * totmass * 2.0 * kinetic / msst_qmass; omega[sd] = -1.0 * sqrt(fac1); double fac2 = omega[sd] / v0; @@ -65,7 +64,7 @@ void MSST::setup(ModuleESolver::ESolver* p_esolver, const std::string& global_re } } - MD_func::compute_stress(ucell, vel, allmass, mdp.cal_stress, virial, stress); + MD_func::compute_stress(ucell, vel, allmass, cal_stress, virial, stress); t_current = MD_func::current_temp(kinetic, ucell.nat, frozen_freedom_, allmass, vel); } @@ -80,7 +79,7 @@ void MSST::first_half(std::ofstream& ofs) ModuleBase::timer::tick("MSST", "first_half"); const int sd = mdp.msst_direction; - const double dthalf = 0.5 * mdp.md_dt; + const double dthalf = 0.5 * md_dt; double vol; energy_ = potential + kinetic; @@ -135,21 +134,21 @@ void MSST::second_half(void) ModuleBase::timer::tick("MSST", "second_half"); const int sd = mdp.msst_direction; - const double dthalf = 0.5 * mdp.md_dt; + const double dthalf = 0.5 * md_dt; energy_ = potential + kinetic; /// propagate velocities 1/2 step propagate_vel(); vsum = vel_sum(); - MD_func::compute_stress(ucell, vel, allmass, mdp.cal_stress, virial, stress); + MD_func::compute_stress(ucell, vel, allmass, cal_stress, virial, stress); t_current = MD_func::current_temp(kinetic, ucell.nat, frozen_freedom_, allmass, vel); /// propagate the time derivative of volume 1/2 step propagate_voldot(); /// calculate Lagrangian position - lag_pos -= mdp.msst_vel * ucell.omega / v0 * mdp.md_dt; + lag_pos -= msst_vel * ucell.omega / v0 * md_dt; ModuleBase::timer::tick("MSST", "second_half"); @@ -166,14 +165,14 @@ void MSST::print_md(std::ofstream& ofs, const bool& cal_stress) void MSST::write_restart(const std::string& global_out_dir) { - if (!mdp.my_rank) + if (!my_rank) { std::stringstream ssc; ssc << global_out_dir << "Restart_md.dat"; std::ofstream file(ssc.str().c_str()); file << step_ + step_rst_ << std::endl; - file << mdp.md_tfirst << std::endl; + file << md_tfirst << std::endl; file << omega[mdp.msst_direction] << std::endl; file << e0 << std::endl; file << v0 << std::endl; @@ -194,7 +193,7 @@ void MSST::restart(const std::string& global_readin_dir) { bool ok = true; - if (!mdp.my_rank) + if (!my_rank) { std::stringstream ssc; ssc << global_readin_dir << "Restart_md.dat"; @@ -207,7 +206,7 @@ void MSST::restart(const std::string& global_readin_dir) if (ok) { - file >> step_rst_ >> mdp.md_tfirst >> omega[mdp.msst_direction] >> e0 >> v0 >> p0 >> lag_pos; + file >> step_rst_ >> md_tfirst >> omega[mdp.msst_direction] >> e0 >> v0 >> p0 >> lag_pos; file.close(); } } @@ -223,7 +222,7 @@ void MSST::restart(const std::string& global_readin_dir) #ifdef __MPI MPI_Bcast(&step_rst_, 1, MPI_INT, 0, MPI_COMM_WORLD); - MPI_Bcast(&mdp.md_tfirst, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&md_tfirst, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(&omega[mdp.msst_direction], 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(&e0, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(&v0, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); @@ -269,11 +268,11 @@ void MSST::rescale(std::ofstream& ofs, const double& volume) void MSST::propagate_vel(void) { - if (mdp.my_rank == 0) + if (my_rank == 0) { const int sd = mdp.msst_direction; - const double dthalf = 0.5 * mdp.md_dt; - const double fac = mdp.msst_vis * pow(omega[sd], 2) / (vsum * ucell.omega); + const double dthalf = 0.5 * md_dt; + const double fac = msst_vis * pow(omega[sd], 2) / (vsum * ucell.omega); for (int i = 0; i < ucell.nat; ++i) { @@ -310,11 +309,11 @@ void MSST::propagate_vel(void) void MSST::propagate_voldot(void) { const int sd = mdp.msst_direction; - const double dthalf = 0.5 * mdp.md_dt; + const double dthalf = 0.5 * md_dt; double p_current = stress(sd, sd); - double p_msst = mdp.msst_vel * mdp.msst_vel * totmass * (v0 - ucell.omega) / (v0 * v0); - double const_A = totmass * (p_current - p0 - p_msst) / mdp.msst_qmass; - double const_B = totmass * mdp.msst_vis / (mdp.msst_qmass * ucell.omega); + double p_msst = msst_vel * msst_vel * totmass * (v0 - ucell.omega) / (v0 * v0); + double const_A = totmass * (p_current - p0 - p_msst) / msst_qmass; + double const_B = totmass * msst_vis / (msst_qmass * ucell.omega); /// prevent the increase of volume if (ucell.omega > v0 && const_A > 0) diff --git a/source/module_md/msst.h b/source/module_md/msst.h index d3143c0b0e..6b42016955 100644 --- a/source/module_md/msst.h +++ b/source/module_md/msst.h @@ -14,7 +14,7 @@ class MSST : public MD_base { public: - MSST(MD_para& MD_para_in, UnitCell& unit_in); + MSST(const Parameter& param_in, UnitCell& unit_in); ~MSST(); private: @@ -61,6 +61,9 @@ class MSST : public MD_base double totmass; ///< total mass of the cell double lag_pos; ///< Lagrangian location of cell double vsum; ///< sum over v^2 + double msst_vel; ///< shock msst_vel (\AA/fs) + double msst_qmass; ///< cell mass-like parameter (mass^2/length^4) + double msst_vis; ///< artificial msst_vis (mass/length/time) }; #endif diff --git a/source/module_md/nhchain.cpp b/source/module_md/nhchain.cpp index db4dce9d3a..215e94327a 100644 --- a/source/module_md/nhchain.cpp +++ b/source/module_md/nhchain.cpp @@ -6,17 +6,19 @@ #endif #include "module_base/timer.h" -Nose_Hoover::Nose_Hoover(MD_para& MD_para_in, UnitCell& unit_in) : MD_base(MD_para_in, unit_in) +Nose_Hoover::Nose_Hoover(const Parameter& param_in, UnitCell& unit_in) : MD_base(param_in, unit_in) { const double unit_transform = ModuleBase::HARTREE_SI / pow(ModuleBase::BOHR_RADIUS_SI, 3) * 1.0e-8; + md_tfreq = mdp.md_tfreq * ModuleBase::AU_to_FS; + assert(unit_transform>0.0); - mdp.md_pfirst /= unit_transform; - mdp.md_plast /= unit_transform; - mdp.md_pfreq *= ModuleBase::AU_to_FS; + md_pfirst = mdp.md_pfirst / unit_transform; + md_plast = mdp.md_plast / unit_transform; + md_pfreq = mdp.md_pfreq * ModuleBase::AU_to_FS; - if (mdp.md_tfirst == 0) + if (md_tfirst == 0) { ModuleBase::WARNING_QUIT("Nose_Hoover", " md_tfirst must be larger than 0 in NHC"); } @@ -32,10 +34,9 @@ Nose_Hoover::Nose_Hoover(MD_para& MD_para_in, UnitCell& unit_in) : MD_base(MD_pa /// determine the NPT methods if (mdp.md_pmode == "iso") { - mdp.md_pcouple = "xyz"; - pstart[0] = pstart[1] = pstart[2] = mdp.md_pfirst; - pstop[0] = pstop[1] = pstop[2] = mdp.md_plast; - pfreq[0] = pfreq[1] = pfreq[2] = mdp.md_pfreq; + pstart[0] = pstart[1] = pstart[2] = md_pfirst; + pstop[0] = pstop[1] = pstop[2] = md_plast; + pfreq[0] = pfreq[1] = pfreq[2] = md_pfreq; pflag[0] = pflag[1] = pflag[2] = 1; } else if (mdp.md_pmode == "aniso") @@ -44,9 +45,9 @@ Nose_Hoover::Nose_Hoover(MD_para& MD_para_in, UnitCell& unit_in) : MD_base(MD_pa { ModuleBase::WARNING_QUIT("Nose_Hoover", "md_pcouple==xyz will convert aniso to iso!"); } - pstart[0] = pstart[1] = pstart[2] = mdp.md_pfirst; - pstop[0] = pstop[1] = pstop[2] = mdp.md_plast; - pfreq[0] = pfreq[1] = pfreq[2] = mdp.md_pfreq; + pstart[0] = pstart[1] = pstart[2] = md_pfirst; + pstop[0] = pstop[1] = pstop[2] = md_plast; + pfreq[0] = pfreq[1] = pfreq[2] = md_pfreq; pflag[0] = pflag[1] = pflag[2] = 1; } /** @@ -62,14 +63,14 @@ Nose_Hoover::Nose_Hoover(MD_para& MD_para_in, UnitCell& unit_in) : MD_base(MD_pa { ModuleBase::WARNING_QUIT("Nose_Hoover", "the lattice must be lower-triangular when md_pmode == tri!"); } - pstart[0] = pstart[1] = pstart[2] = mdp.md_pfirst; - pstop[0] = pstop[1] = pstop[2] = mdp.md_plast; - pfreq[0] = pfreq[1] = pfreq[2] = mdp.md_pfreq; + pstart[0] = pstart[1] = pstart[2] = md_pfirst; + pstop[0] = pstop[1] = pstop[2] = md_plast; + pfreq[0] = pfreq[1] = pfreq[2] = md_pfreq; pflag[0] = pflag[1] = pflag[2] = 1; pstart[3] = pstart[4] = pstart[5] = 0; pstop[3] = pstop[4] = pstop[5] = 0; - pfreq[3] = pfreq[4] = pfreq[5] = mdp.md_pfreq; + pfreq[3] = pfreq[4] = pfreq[5] = md_pfreq; pflag[3] = pflag[4] = pflag[5] = 1; } else @@ -164,13 +165,13 @@ void Nose_Hoover::setup(ModuleESolver::ESolver* p_esolver, const std::string& gl } /// determine target temperature - t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, mdp.md_tfirst, mdp.md_tlast); + t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, md_tfirst, md_tlast); /// init thermostats coupled with particles - mass_eta[0] = tdof * t_target / mdp.md_tfreq / mdp.md_tfreq; + mass_eta[0] = tdof * t_target / md_tfreq / md_tfreq; for (int m = 1; m < mdp.md_tchain; ++m) { - mass_eta[m] = t_target / mdp.md_tfreq / mdp.md_tfreq; + mass_eta[m] = t_target / md_tfreq / md_tfreq; g_eta[m] = (mass_eta[m - 1] * v_eta[m - 1] * v_eta[m - 1] - t_target) / mass_eta[m]; } @@ -197,10 +198,10 @@ void Nose_Hoover::setup(ModuleESolver::ESolver* p_esolver, const std::string& gl /// init thermostats coupled with barostat if (mdp.md_pchain) { - mass_peta[0] = t_target / mdp.md_pfreq / mdp.md_pfreq; + mass_peta[0] = t_target / md_pfreq / md_pfreq; for (int m = 1; m < mdp.md_pchain; ++m) { - mass_peta[m] = t_target / mdp.md_pfreq / mdp.md_pfreq; + mass_peta[m] = t_target / md_pfreq / md_pfreq; g_peta[m] = (mass_peta[m - 1] * v_peta[m - 1] * v_peta[m - 1] - t_target) / mass_peta[m]; } } @@ -223,7 +224,7 @@ void Nose_Hoover::first_half(std::ofstream& ofs) } /// update target T - t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, mdp.md_tfirst, mdp.md_tlast); + t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, md_tfirst, md_tlast); /// update thermostats coupled with particles particle_thermo(); @@ -232,7 +233,7 @@ void Nose_Hoover::first_half(std::ofstream& ofs) { /// update temperature and stress due to velocity rescaling t_current = MD_func::current_temp(kinetic, ucell.nat, frozen_freedom_, allmass, vel); - MD_func::compute_stress(ucell, vel, allmass, mdp.cal_stress,virial, stress); + MD_func::compute_stress(ucell, vel, allmass, cal_stress, virial, stress); /// couple stress component due to md_pcouple couple_stress(); @@ -291,7 +292,7 @@ void Nose_Hoover::second_half(void) if (npt_flag) { /// update stress due to velocity rescaling - MD_func::compute_stress(ucell, vel, allmass, mdp.cal_stress,virial, stress); + MD_func::compute_stress(ucell, vel, allmass, cal_stress, virial, stress); /// couple stress component due to md_pcouple couple_stress(); @@ -322,14 +323,14 @@ void Nose_Hoover::print_md(std::ofstream& ofs, const bool& cal_stress) void Nose_Hoover::write_restart(const std::string& global_out_dir) { - if (!mdp.my_rank) + if (!my_rank) { std::stringstream ssc; ssc << global_out_dir << "Restart_md.dat"; std::ofstream file(ssc.str().c_str()); file << step_ + step_rst_ << std::endl; - file << mdp.md_tfirst << std::endl; + file << md_tfirst << std::endl; file << mdp.md_tchain << std::endl; for (int i = 0; i < mdp.md_tchain; ++i) { @@ -376,7 +377,7 @@ void Nose_Hoover::restart(const std::string& global_readin_dir) bool ok2 = true; bool ok3 = true; - if (!mdp.my_rank) + if (!my_rank) { std::stringstream ssc; ssc << global_readin_dir << "Restart_md.dat"; @@ -390,7 +391,7 @@ void Nose_Hoover::restart(const std::string& global_readin_dir) if (ok) { double Mnum; - file >> step_rst_ >> mdp.md_tfirst >> Mnum; + file >> step_rst_ >> md_tfirst >> Mnum; if (Mnum != mdp.md_tchain) { @@ -461,7 +462,7 @@ void Nose_Hoover::restart(const std::string& global_readin_dir) #ifdef __MPI MPI_Bcast(&step_rst_, 1, MPI_INT, 0, MPI_COMM_WORLD); - MPI_Bcast(&mdp.md_tfirst, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&md_tfirst, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(eta, mdp.md_tchain, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(v_eta, mdp.md_tchain, MPI_DOUBLE, 0, MPI_COMM_WORLD); if (npt_flag) @@ -476,10 +477,10 @@ void Nose_Hoover::restart(const std::string& global_readin_dir) void Nose_Hoover::particle_thermo() { /// update mass_eta - mass_eta[0] = tdof * t_target / mdp.md_tfreq / mdp.md_tfreq; + mass_eta[0] = tdof * t_target / md_tfreq / md_tfreq; for (int m = 1; m < mdp.md_tchain; ++m) { - mass_eta[m] = t_target / mdp.md_tfreq / mdp.md_tfreq; + mass_eta[m] = t_target / md_tfreq / md_tfreq; } /// propogate g_eta @@ -500,7 +501,7 @@ void Nose_Hoover::particle_thermo() { for (int j = 0; j < nys; ++j) { - double delta = w[j] * mdp.md_dt / nc_tchain; + double delta = w[j] * md_dt / nc_tchain; /// propogate v_eta for (int m = mdp.md_tchain - 1; m >= 0; --m) @@ -589,7 +590,7 @@ void Nose_Hoover::baro_thermo() { for (int j = 0; j < nys; ++j) { - double delta = w[j] * mdp.md_dt / nc_pchain; + double delta = w[j] * md_dt / nc_pchain; /// propogate v_peta for (int m = mdp.md_pchain - 1; m >= 0; --m) @@ -668,7 +669,7 @@ void Nose_Hoover::update_baro() if (pflag[i]) { g_omega = (p_current[i] - p_hydro) * ucell.omega / mass_omega[i] + term_one / mass_omega[i]; - v_omega[i] += g_omega * mdp.md_dt / 2.0; + v_omega[i] += g_omega * md_dt / 2.0; term_two += v_omega[i]; } } @@ -679,7 +680,7 @@ void Nose_Hoover::update_baro() if (pflag[i]) { g_omega = p_current[i] * ucell.omega / mass_omega[i]; - v_omega[i] += g_omega * mdp.md_dt / 2.0; + v_omega[i] += g_omega * md_dt / 2.0; } } @@ -691,7 +692,7 @@ void Nose_Hoover::vel_baro() double factor[3]; for (int i = 0; i < 3; ++i) { - factor[i] = exp(-(v_omega[i] + mtk_term) * mdp.md_dt / 4); + factor[i] = exp(-(v_omega[i] + mtk_term) * md_dt / 4); } for (int i = 0; i < ucell.nat; ++i) @@ -704,11 +705,11 @@ void Nose_Hoover::vel_baro() /// Note: I am not sure whether fixed atoms should update here if (ionmbl[i][0]) { - vel[i][0] -= (vel[i][1] * v_omega[5] + vel[i][2] * v_omega[4]) * mdp.md_dt / 2; + vel[i][0] -= (vel[i][1] * v_omega[5] + vel[i][2] * v_omega[4]) * md_dt / 2; } if (ionmbl[i][1]) { - vel[i][1] -= vel[i][2] * v_omega[3] * mdp.md_dt / 2; + vel[i][1] -= vel[i][2] * v_omega[3] * md_dt / 2; } for (int j = 0; j < 3; ++j) @@ -725,7 +726,7 @@ void Nose_Hoover::update_volume(std::ofstream& ofs) /// tri mode, off-diagonal components, first half if (pflag[4]) { - factor = exp(v_omega[0] * mdp.md_dt / 16); + factor = exp(v_omega[0] * md_dt / 16); ucell.latvec.e31 *= factor; ucell.latvec.e31 += (v_omega[5] * ucell.latvec.e32 + v_omega[4] * ucell.latvec.e33); ucell.latvec.e31 *= factor; @@ -733,7 +734,7 @@ void Nose_Hoover::update_volume(std::ofstream& ofs) if (pflag[3]) { - factor = exp(v_omega[1] * mdp.md_dt / 8); + factor = exp(v_omega[1] * md_dt / 8); ucell.latvec.e32 *= factor; ucell.latvec.e32 += (v_omega[3] * ucell.latvec.e33); ucell.latvec.e32 *= factor; @@ -741,7 +742,7 @@ void Nose_Hoover::update_volume(std::ofstream& ofs) if (pflag[5]) { - factor = exp(v_omega[0] * mdp.md_dt / 8); + factor = exp(v_omega[0] * md_dt / 8); ucell.latvec.e21 *= factor; ucell.latvec.e21 += (v_omega[5] * ucell.latvec.e22); ucell.latvec.e21 *= factor; @@ -749,7 +750,7 @@ void Nose_Hoover::update_volume(std::ofstream& ofs) if (pflag[4]) { - factor = exp(v_omega[0] * mdp.md_dt / 16); + factor = exp(v_omega[0] * md_dt / 16); ucell.latvec.e31 *= factor; ucell.latvec.e31 += (v_omega[5] * ucell.latvec.e32 + v_omega[4] * ucell.latvec.e33); ucell.latvec.e31 *= factor; @@ -758,26 +759,26 @@ void Nose_Hoover::update_volume(std::ofstream& ofs) /// Diagonal components if (pflag[0]) { - factor = exp(v_omega[0] * mdp.md_dt / 2); + factor = exp(v_omega[0] * md_dt / 2); ucell.latvec.e11 *= factor; } if (pflag[1]) { - factor = exp(v_omega[1] * mdp.md_dt / 2); + factor = exp(v_omega[1] * md_dt / 2); ucell.latvec.e22 *= factor; } if (pflag[2]) { - factor = exp(v_omega[2] * mdp.md_dt / 2); + factor = exp(v_omega[2] * md_dt / 2); ucell.latvec.e33 *= factor; } /// tri mode, off-diagonal components, second half if (pflag[4]) { - factor = exp(v_omega[0] * mdp.md_dt / 16); + factor = exp(v_omega[0] * md_dt / 16); ucell.latvec.e31 *= factor; ucell.latvec.e31 += (v_omega[5] * ucell.latvec.e32 + v_omega[4] * ucell.latvec.e33); ucell.latvec.e31 *= factor; @@ -785,7 +786,7 @@ void Nose_Hoover::update_volume(std::ofstream& ofs) if (pflag[3]) { - factor = exp(v_omega[1] * mdp.md_dt / 8); + factor = exp(v_omega[1] * md_dt / 8); ucell.latvec.e32 *= factor; ucell.latvec.e32 += (v_omega[3] * ucell.latvec.e33); ucell.latvec.e32 *= factor; @@ -793,7 +794,7 @@ void Nose_Hoover::update_volume(std::ofstream& ofs) if (pflag[5]) { - factor = exp(v_omega[0] * mdp.md_dt / 8); + factor = exp(v_omega[0] * md_dt / 8); ucell.latvec.e21 *= factor; ucell.latvec.e21 += (v_omega[5] * ucell.latvec.e22); ucell.latvec.e21 *= factor; @@ -801,7 +802,7 @@ void Nose_Hoover::update_volume(std::ofstream& ofs) if (pflag[4]) { - factor = exp(v_omega[0] * mdp.md_dt / 16); + factor = exp(v_omega[0] * md_dt / 16); ucell.latvec.e31 *= factor; ucell.latvec.e31 += (v_omega[5] * ucell.latvec.e32 + v_omega[4] * ucell.latvec.e33); ucell.latvec.e31 *= factor; diff --git a/source/module_md/nhchain.h b/source/module_md/nhchain.h index 2de76293d5..b804d8784b 100644 --- a/source/module_md/nhchain.h +++ b/source/module_md/nhchain.h @@ -12,7 +12,7 @@ class Nose_Hoover : public MD_base { public: - Nose_Hoover(MD_para& MD_para_in, UnitCell& unit_in); + Nose_Hoover(const Parameter& param_in, UnitCell& unit_in); ~Nose_Hoover(); private: @@ -94,6 +94,10 @@ class Nose_Hoover : public MD_base double* v_peta; ///< velocity of thermostats coupled with barostat double* g_peta; ///< acceleration of thermostats coupled with barostat double mtk_term; ///< mtk correction + double md_tfreq; ///< Oscillation frequency, used to determine qmass of thermostats coupled with particles + double md_pfirst; ///< Initial pressure + double md_plast; ///< Final pressure + double md_pfreq; ///< Oscillation frequency, used to determine qmass of thermostats coupled with barostat }; #endif \ No newline at end of file diff --git a/source/module_md/run_md.cpp b/source/module_md/run_md.cpp index 84b04e690a..45a6c98a4a 100644 --- a/source/module_md/run_md.cpp +++ b/source/module_md/run_md.cpp @@ -13,32 +13,32 @@ namespace Run_MD { -void md_line(UnitCell& unit_in, ModuleESolver::ESolver* p_esolver, MD_para& md_para) +void md_line(UnitCell& unit_in, ModuleESolver::ESolver* p_esolver, const Parameter& param_in) { ModuleBase::TITLE("Run_MD", "md_line"); ModuleBase::timer::tick("Run_MD", "md_line"); /// determine the md_type MD_base* mdrun; - if (md_para.md_type == "fire") + if (param_in.mdp.md_type == "fire") { - mdrun = new FIRE(md_para, unit_in); + mdrun = new FIRE(param_in, unit_in); } - else if ((md_para.md_type == "nvt" && md_para.md_thermostat == "nhc") || md_para.md_type == "npt") + else if ((param_in.mdp.md_type == "nvt" && param_in.mdp.md_thermostat == "nhc") || param_in.mdp.md_type == "npt") { - mdrun = new Nose_Hoover(md_para, unit_in); + mdrun = new Nose_Hoover(param_in, unit_in); } - else if (md_para.md_type == "nve" || md_para.md_type == "nvt") + else if (param_in.mdp.md_type == "nve" || param_in.mdp.md_type == "nvt") { - mdrun = new Verlet(md_para, unit_in); + mdrun = new Verlet(param_in, unit_in); } - else if (md_para.md_type == "langevin") + else if (param_in.mdp.md_type == "langevin") { - mdrun = new Langevin(md_para, unit_in); + mdrun = new Langevin(param_in, unit_in); } - else if (md_para.md_type == "msst") + else if (param_in.mdp.md_type == "msst") { - mdrun = new MSST(md_para, unit_in); + mdrun = new MSST(param_in, unit_in); } else { @@ -46,7 +46,7 @@ void md_line(UnitCell& unit_in, ModuleESolver::ESolver* p_esolver, MD_para& md_p } /// md cycle - while ((mdrun->step_ + mdrun->step_rst_) <= md_para.md_nstep && !mdrun->stop) + while ((mdrun->step_ + mdrun->step_rst_) <= param_in.mdp.md_nstep && !mdrun->stop) { if (mdrun->step_ == 0) { @@ -63,7 +63,7 @@ void md_line(UnitCell& unit_in, ModuleESolver::ESolver* p_esolver, MD_para& md_p unit_in, mdrun->potential, mdrun->force, - md_para.cal_stress, + param_in.inp.cal_stress, mdrun->virial); mdrun->second_half(); @@ -71,7 +71,7 @@ void md_line(UnitCell& unit_in, ModuleESolver::ESolver* p_esolver, MD_para& md_p MD_func::compute_stress(unit_in, mdrun->vel, mdrun->allmass, - md_para.cal_stress, + param_in.inp.cal_stress, mdrun->virial, mdrun->stress); mdrun->t_current = MD_func::current_temp(mdrun->kinetic, @@ -81,20 +81,20 @@ void md_line(UnitCell& unit_in, ModuleESolver::ESolver* p_esolver, MD_para& md_p mdrun->vel); } - if ((mdrun->step_ + mdrun->step_rst_) % md_para.md_dumpfreq == 0) + if ((mdrun->step_ + mdrun->step_rst_) % param_in.mdp.md_dumpfreq == 0) { mdrun->print_md(GlobalV::ofs_running, GlobalV::CAL_STRESS); MD_func::dump_info(mdrun->step_ + mdrun->step_rst_, GlobalV::global_out_dir, unit_in, - md_para, + param_in, mdrun->virial, mdrun->force, mdrun->vel); } - if ((mdrun->step_ + mdrun->step_rst_) % md_para.md_restartfreq == 0) + if ((mdrun->step_ + mdrun->step_rst_) % param_in.mdp.md_restartfreq == 0) { unit_in.update_vel(mdrun->vel); std::stringstream file; diff --git a/source/module_md/run_md.h b/source/module_md/run_md.h index 4e877275e0..fdc27f3c0d 100644 --- a/source/module_md/run_md.h +++ b/source/module_md/run_md.h @@ -2,7 +2,7 @@ #define RUN_MD_H #include "module_esolver/esolver.h" -#include "module_parameter/md_parameter.h" +#include "module_parameter/parameter.h" /** * @brief the md loop line @@ -17,7 +17,7 @@ namespace Run_MD * @param p_esolver energy solver * @param md_para input parameters used in md */ -void md_line(UnitCell& unit_in, ModuleESolver::ESolver* p_esolver, MD_para& md_para); +void md_line(UnitCell& unit_in, ModuleESolver::ESolver* p_esolver, const Parameter& param_in); } // namespace Run_MD #endif \ No newline at end of file diff --git a/source/module_md/test/md_func_test.cpp b/source/module_md/test/md_func_test.cpp index c3368510e5..44ede5fa53 100644 --- a/source/module_md/test/md_func_test.cpp +++ b/source/module_md/test/md_func_test.cpp @@ -210,7 +210,7 @@ TEST_F(MD_func_test, compute_stress) TEST_F(MD_func_test, dump_info) { - MD_func::dump_info(0, GlobalV::global_out_dir, ucell, input.mdp, virial, force, vel); + MD_func::dump_info(0, GlobalV::MY_RANK, GlobalV::global_out_dir, ucell, input.mdp, virial, force, vel); std::ifstream ifs("MD_dump"); std::string output_str; getline(ifs, output_str); @@ -256,7 +256,7 @@ TEST_F(MD_func_test, dump_info) ifs.close(); // append - MD_func::dump_info(1, GlobalV::global_out_dir, ucell, input.mdp, virial, force, vel); + MD_func::dump_info(1, GlobalV::MY_RANK, GlobalV::global_out_dir, ucell, input.mdp, virial, force, vel); std::ifstream ifs2("MD_dump"); getline(ifs2, output_str); EXPECT_THAT(output_str, testing::HasSubstr("MDSTEP: 0")); diff --git a/source/module_md/test/setcell.h b/source/module_md/test/setcell.h index dbd0e43b57..0b3d60f2f7 100644 --- a/source/module_md/test/setcell.h +++ b/source/module_md/test/setcell.h @@ -134,7 +134,7 @@ class Setcell input.mdp.dump_virial = true; input.mdp.dump_force = true; input.mdp.dump_vel = true; - input.mdp.cal_stress = true; + input.cal_stress = true; input.mdp.md_restart = false; input.mdp.md_dt = 1; diff --git a/source/module_md/verlet.cpp b/source/module_md/verlet.cpp index c8db859593..2eec7a0e20 100644 --- a/source/module_md/verlet.cpp +++ b/source/module_md/verlet.cpp @@ -3,7 +3,7 @@ #include "md_func.h" #include "module_base/timer.h" -Verlet::Verlet(MD_para& MD_para_in, UnitCell& unit_in) : MD_base(MD_para_in, unit_in) +Verlet::Verlet(const Parameter& param_in, UnitCell& unit_in) : MD_base(param_in, unit_in) { } @@ -57,7 +57,7 @@ void Verlet::apply_thermostat(void) } else if (mdp.md_thermostat == "rescaling") { - t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, mdp.md_tfirst, mdp.md_tlast); + t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, md_tfirst, md_tlast); if (std::abs(t_target - t_current) * ModuleBase::Hartree_to_K > mdp.md_tolerance) { thermalize(0, t_current, t_target); @@ -67,20 +67,20 @@ void Verlet::apply_thermostat(void) { if ((step_ + step_rst_) % mdp.md_nraise == 0) { - t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, mdp.md_tfirst, mdp.md_tlast); + t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, md_tfirst, md_tlast); thermalize(0, t_current, t_target); } } else if (mdp.md_thermostat == "anderson") { - if (mdp.my_rank == 0) + if (my_rank == 0) { double deviation; for (int i = 0; i < ucell.nat; ++i) { if (static_cast(std::rand()) / RAND_MAX <= 1.0 / mdp.md_nraise) { - deviation = sqrt(mdp.md_tlast / allmass[i]); + deviation = sqrt(md_tlast / allmass[i]); for (int k = 0; k < 3; ++k) { if (ionmbl[i][k]) @@ -97,7 +97,7 @@ void Verlet::apply_thermostat(void) } else if (mdp.md_thermostat == "berendsen") { - t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, mdp.md_tfirst, mdp.md_tlast); + t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, md_tfirst, md_tlast); thermalize(mdp.md_nraise, t_current, t_target); } else diff --git a/source/module_md/verlet.h b/source/module_md/verlet.h index d9c6d39030..72e0edcd6e 100644 --- a/source/module_md/verlet.h +++ b/source/module_md/verlet.h @@ -10,7 +10,7 @@ class Verlet : public MD_base { public: - Verlet(MD_para& MD_para_in, UnitCell& unit_in); + Verlet(const Parameter& param_in, UnitCell& unit_in); ~Verlet(); private: diff --git a/source/module_parameter/md_parameter.h b/source/module_parameter/md_parameter.h index 8c72c565bc..7ca6c2d8ca 100644 --- a/source/module_parameter/md_parameter.h +++ b/source/module_parameter/md_parameter.h @@ -58,12 +58,6 @@ struct MD_para ///< not. liuyu 2023-03-01 bool dump_virial = true; ///< output lattice virial into the file MD_dump or ///< not. liuyu 2023-03-01 - - // !!! - // They are repeated in the MD_para and Parameter, which is not good. - double force_thr = 1.0e-3; ///< force convergence threshold in FIRE method - bool cal_stress = false; ///< whether calculate stress - int my_rank = 0; ///< MPI rank of the processor }; #endif // MD_PARA_H \ No newline at end of file diff --git a/source/module_parameter/parameter.cpp b/source/module_parameter/parameter.cpp index 86da7b6580..45c9a04477 100644 --- a/source/module_parameter/parameter.cpp +++ b/source/module_parameter/parameter.cpp @@ -6,7 +6,6 @@ void Parameter::set_rank_nproc(const int& myrank, const int& nproc) { sys.myrank = myrank; sys.nproc = nproc; - input.mdp.my_rank = myrank; } void Parameter::set_start_time(const std::time_t& start_time)