From 10ee8a61328feb453af66139a6e085cade8d795f Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Fri, 5 May 2023 15:50:16 +0800 Subject: [PATCH 1/5] Refactor: delete ueedless head files in module_md --- .../module_hamilt_general/module_vdw/vdwd3_parameters.cpp | 2 ++ source/module_io/input.cpp | 1 + source/module_io/write_input.cpp | 5 +++-- source/module_md/md_base.cpp | 2 -- source/module_md/md_base.h | 2 -- source/module_md/md_func.h | 3 --- source/module_md/md_para.h | 2 +- source/module_md/msst.cpp | 1 - source/module_md/run_md.cpp | 7 +++---- source/module_md/run_md.h | 1 - 10 files changed, 10 insertions(+), 16 deletions(-) diff --git a/source/module_hamilt_general/module_vdw/vdwd3_parameters.cpp b/source/module_hamilt_general/module_vdw/vdwd3_parameters.cpp index aadab9cd5c..a94991557b 100644 --- a/source/module_hamilt_general/module_vdw/vdwd3_parameters.cpp +++ b/source/module_hamilt_general/module_vdw/vdwd3_parameters.cpp @@ -6,6 +6,8 @@ #include "vdwd3_parameters.h" +#include "module_base/constants.h" + namespace vdw { diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp index a14b275ad5..cd5b689604 100644 --- a/source/module_io/input.cpp +++ b/source/module_io/input.cpp @@ -5,6 +5,7 @@ // #include "global.h" #include "module_io/input.h" +#include "module_base/constants.h" #include "module_base/global_file.h" #include "module_base/global_function.h" #include "module_base/global_variable.h" diff --git a/source/module_io/write_input.cpp b/source/module_io/write_input.cpp index 51c42ca744..a62efb7a7f 100644 --- a/source/module_io/write_input.cpp +++ b/source/module_io/write_input.cpp @@ -1,6 +1,7 @@ #include "module_io/input.h" -#include "../module_base/global_function.h" -#include "../module_base/global_variable.h" +#include "module_base/constants.h" +#include "module_base/global_function.h" +#include "module_base/global_variable.h" void Input::Print(const std::string &fn) const { diff --git a/source/module_md/md_base.cpp b/source/module_md/md_base.cpp index 0d8dc5d09a..e3b706c20c 100644 --- a/source/module_md/md_base.cpp +++ b/source/module_md/md_base.cpp @@ -4,8 +4,6 @@ #ifdef __MPI #include "mpi.h" #endif -#include "module_base/timer.h" -#include "module_esolver/esolver.h" #include "module_io/print_info.h" MDrun::MDrun(MD_parameters &MD_para_in, UnitCell &unit_in) : mdp(MD_para_in), ucell(unit_in) diff --git a/source/module_md/md_base.h b/source/module_md/md_base.h index 2f3e64dbaa..cfba171758 100644 --- a/source/module_md/md_base.h +++ b/source/module_md/md_base.h @@ -2,8 +2,6 @@ #define MDRUN_H #include "md_para.h" -#include "module_base/matrix.h" -#include "module_cell/unitcell.h" #include "module_esolver/esolver.h" class MDrun diff --git a/source/module_md/md_func.h b/source/module_md/md_func.h index 972f3e772b..ca50fe44a5 100644 --- a/source/module_md/md_func.h +++ b/source/module_md/md_func.h @@ -1,9 +1,6 @@ #ifndef MD_FUNC_H #define MD_FUNC_H -#include "md_para.h" -#include "module_base/matrix.h" -#include "module_cell/unitcell.h" #include "module_esolver/esolver.h" namespace MD_func diff --git a/source/module_md/md_para.h b/source/module_md/md_para.h index 6cb6b28dba..c869dd07ac 100644 --- a/source/module_md/md_para.h +++ b/source/module_md/md_para.h @@ -1,7 +1,7 @@ #ifndef MD_PARAMETERS_H #define MD_PARAMETERS_H -#include "module_base/constants.h" +#include class MD_parameters { diff --git a/source/module_md/msst.cpp b/source/module_md/msst.cpp index 41b6cc75fe..7ddeef108f 100644 --- a/source/module_md/msst.cpp +++ b/source/module_md/msst.cpp @@ -5,7 +5,6 @@ #include "mpi.h" #endif #include "module_base/timer.h" -#include "module_esolver/esolver.h" MSST::MSST(MD_parameters &MD_para_in, UnitCell &unit_in) : MDrun(MD_para_in, unit_in) { diff --git a/source/module_md/run_md.cpp b/source/module_md/run_md.cpp index b2143bf935..f7b1d91dd2 100644 --- a/source/module_md/run_md.cpp +++ b/source/module_md/run_md.cpp @@ -1,13 +1,12 @@ #include "run_md.h" +#include "fire.h" #include "langevin.h" #include "md_func.h" -#include "msst.h" -#include "nhchain.h" -#include "fire.h" #include "module_base/timer.h" -#include "module_io/input.h" #include "module_io/print_info.h" +#include "msst.h" +#include "nhchain.h" #include "verlet.h" Run_MD::Run_MD() diff --git a/source/module_md/run_md.h b/source/module_md/run_md.h index 06dbff5dd3..9d80efffde 100644 --- a/source/module_md/run_md.h +++ b/source/module_md/run_md.h @@ -1,7 +1,6 @@ #ifndef RUN_MD_CLASSIC_H #define RUN_MD_CLASSIC_H -#include "module_cell/unitcell.h" #include "module_esolver/esolver.h" class Run_MD From 734050f74518338120687550924f4b1187c8fe7c Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Fri, 5 May 2023 20:19:23 +0800 Subject: [PATCH 2/5] Refactor: delete GlobalV::MD_NSTEP --- source/driver_run.cpp | 2 +- source/module_base/global_variable.cpp | 1 - source/module_base/global_variable.h | 3 +- source/module_io/input.cpp | 33 ++++----- source/module_io/input.h | 12 ++-- source/module_io/input_conv.cpp | 1 - source/module_io/test/input_conv_test.cpp | 5 +- source/module_io/test/input_test.cpp | 16 ++--- source/module_io/test/input_test_para.cpp | 8 +-- source/module_io/write_input.cpp | 15 +++- source/module_md/fire.cpp | 41 ++++++----- source/module_md/fire.h | 17 +++-- source/module_md/langevin.cpp | 2 +- source/module_md/md_base.cpp | 36 +++++----- source/module_md/md_base.h | 16 ++--- source/module_md/md_func.cpp | 16 ++--- source/module_md/md_func.h | 4 +- source/module_md/md_para.h | 8 +++ source/module_md/msst.cpp | 84 +++++++++++++---------- source/module_md/msst.h | 14 ++-- source/module_md/nhchain.cpp | 44 ++++++------ source/module_md/nhchain.h | 14 ++-- source/module_md/run_md.cpp | 34 ++++----- source/module_md/run_md.h | 3 +- source/module_md/test/md_func_test.cpp | 2 +- source/module_md/test/setcell.h | 6 +- source/module_md/verlet.cpp | 36 +++++----- source/module_md/verlet.h | 14 ++-- 28 files changed, 255 insertions(+), 232 deletions(-) diff --git a/source/driver_run.cpp b/source/driver_run.cpp index 6e14250de6..224e725a96 100644 --- a/source/driver_run.cpp +++ b/source/driver_run.cpp @@ -49,7 +49,7 @@ void Driver::driver_run() if(GlobalV::CALCULATION == "md") { Run_MD run_md; - run_md.md_line(GlobalC::ucell, p_esolver); + run_md.md_line(GlobalC::ucell, p_esolver, INPUT.mdp); } else // scf; cell relaxation; nscf; etc { diff --git a/source/module_base/global_variable.cpp b/source/module_base/global_variable.cpp index 739144f023..9f811c0d9a 100644 --- a/source/module_base/global_variable.cpp +++ b/source/module_base/global_variable.cpp @@ -57,7 +57,6 @@ bool fixed_atoms = false; int OUT_FREQ_ELEC = 0; int OUT_FREQ_ION = 0; int RELAX_NMAX = 20; -int MD_NSTEP = 20; int md_prec_level = 0; int SCF_NMAX = 100; diff --git a/source/module_base/global_variable.h b/source/module_base/global_variable.h index 92cf2a20ab..cb793c2113 100644 --- a/source/module_base/global_variable.h +++ b/source/module_base/global_variable.h @@ -61,8 +61,7 @@ extern bool relax_new; extern bool fixed_atoms; extern int RELAX_NMAX; // 8.3 -extern int SCF_NMAX; // 8.4 -extern int MD_NSTEP; +extern int SCF_NMAX; // 8.4 extern int md_prec_level; // liuyu 2023-03-13 extern std::string BASIS_TYPE; // xiaohui add 2013-09-01 diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp index cd5b689604..bc20758ff3 100644 --- a/source/module_io/input.cpp +++ b/source/module_io/input.cpp @@ -190,9 +190,6 @@ void Input::Default(void) ref_cell_factor = 1.0; symmetry_prec = 1.0e-5; // LiuXh add 2021-08-12, accuracy for symmetry cal_force = 0; - dump_force = true; - dump_vel = true; - dump_virial = true; force_thr = 1.0e-3; force_thr_ev2 = 0; stress_thr = 1.0e-2; // LiuXh add 20180515 @@ -866,18 +863,6 @@ bool Input::Read(const std::string &fn) { read_bool(ifs, cal_force); } - else if (strcmp("dump_force", word) == 0) - { - read_bool(ifs, dump_force); - } - else if (strcmp("dump_vel", word) == 0) - { - read_bool(ifs, dump_vel); - } - else if (strcmp("dump_virial", word) == 0) - { - read_bool(ifs, dump_virial); - } else if (strcmp("force_thr", word) == 0) { read_value(ifs, force_thr); @@ -1463,6 +1448,18 @@ bool Input::Read(const std::string &fn) { read_value(ifs, mdp.pot_file); } + else if (strcmp("dump_force", word) == 0) + { + read_bool(ifs, mdp.dump_force); + } + else if (strcmp("dump_vel", word) == 0) + { + read_bool(ifs, mdp.dump_vel); + } + else if (strcmp("dump_virial", word) == 0) + { + read_bool(ifs, mdp.dump_virial); + } //---------------------------------------------------------- // efield and dipole correction // Yu Liu add 2022-05-18 @@ -2787,9 +2784,6 @@ void Input::Bcast() Parallel_Common::bcast_double(ref_cell_factor); Parallel_Common::bcast_double(symmetry_prec); // LiuXh add 2021-08-12, accuracy for symmetry Parallel_Common::bcast_bool(cal_force); - Parallel_Common::bcast_bool(dump_force); - Parallel_Common::bcast_bool(dump_vel); - Parallel_Common::bcast_bool(dump_virial); Parallel_Common::bcast_double(force_thr); Parallel_Common::bcast_double(force_thr_ev2); Parallel_Common::bcast_double(stress_thr); // LiuXh add 20180515 @@ -2944,6 +2938,9 @@ void Input::Bcast() Parallel_Common::bcast_double(mdp.md_pfirst); Parallel_Common::bcast_double(mdp.md_plast); Parallel_Common::bcast_double(mdp.md_pfreq); + Parallel_Common::bcast_bool(mdp.dump_force); + Parallel_Common::bcast_bool(mdp.dump_vel); + Parallel_Common::bcast_bool(mdp.dump_virial); // Yu Liu add 2022-05-18 Parallel_Common::bcast_bool(efield_flag); Parallel_Common::bcast_bool(dip_cor_flag); diff --git a/source/module_io/input.h b/source/module_io/input.h index 5d7b40d6a2..09e165da12 100644 --- a/source/module_io/input.h +++ b/source/module_io/input.h @@ -1,13 +1,14 @@ #ifndef INPUT_H #define INPUT_H -#include "module_base/vector3.h" -#include "module_md/md_para.h" - #include +#include #include #include +#include "module_base/vector3.h" +#include "module_md/md_para.h" + using namespace std; class Input @@ -45,8 +46,7 @@ class Input int nbands_istate; // number of bands around fermi level for istate calculation. int pw_seed; // random seed for initializing wave functions qianrui 2021-8-12 - bool init_vel; // read velocity from STRU or not liuyu 2021-07-14 - bool dump_vel; // output atomic velocities into the file MD_dump or not. liuyu 2023-03-01 + bool init_vel; // read velocity from STRU or not liuyu 2021-07-14 double ref_cell_factor; // construct a reference cell bigger than the initial cell liuyu 2023-03-21 /* symmetry level: @@ -110,7 +110,6 @@ class Input // Forces //========================================================== bool cal_force; - bool dump_force; // output atomic forces into the file MD_dump or not. liuyu 2023-03-01 double force_thr; // threshold of force in unit (Ry/Bohr) double force_thr_ev2; // invalid force threshold, mohan add 2011-04-17 @@ -122,7 +121,6 @@ class Input double press2; double press3; bool cal_stress; // calculate the stress - bool dump_virial; // output lattice virial into the file MD_dump or not. liuyu 2023-03-01 std::string fixed_axes; // which axes are fixed bool fixed_ibrav; //whether to keep type of lattice; must be used along with latname diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index a7ac785d28..b7bfd95745 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -530,7 +530,6 @@ void Input_Conv::Convert(void) //---------------------------------------------------------- GlobalV::SCF_NMAX = INPUT.scf_nmax; GlobalV::RELAX_NMAX = INPUT.relax_nmax; - GlobalV::MD_NSTEP = INPUT.mdp.md_nstep; GlobalV::md_prec_level = INPUT.mdp.md_prec_level; //---------------------------------------------------------- diff --git a/source/module_io/test/input_conv_test.cpp b/source/module_io/test/input_conv_test.cpp index c39d20e0a6..b4ca754b44 100644 --- a/source/module_io/test/input_conv_test.cpp +++ b/source/module_io/test/input_conv_test.cpp @@ -136,9 +136,8 @@ TEST_F(InputConvTest, Conv) EXPECT_EQ(GlobalV::out_mul,0); EXPECT_EQ(GlobalC::ppcell.cell_factor,1.2); EXPECT_EQ(GlobalV::SCF_NMAX,50); - EXPECT_EQ(GlobalV::RELAX_NMAX,1); - EXPECT_EQ(GlobalV::MD_NSTEP,10); - EXPECT_EQ(GlobalV::OUT_FREQ_ELEC,0); + EXPECT_EQ(GlobalV::RELAX_NMAX, 1); + EXPECT_EQ(GlobalV::OUT_FREQ_ELEC,0); EXPECT_EQ(GlobalV::OUT_FREQ_ION,0); EXPECT_EQ(GlobalV::init_chg,"atomic"); EXPECT_EQ(GlobalV::chg_extrap,"atomic"); diff --git a/source/module_io/test/input_test.cpp b/source/module_io/test/input_test.cpp index 63456a5f53..fb366b848f 100644 --- a/source/module_io/test/input_test.cpp +++ b/source/module_io/test/input_test.cpp @@ -83,10 +83,7 @@ TEST_F(InputTest, Default) EXPECT_FALSE(INPUT.init_vel); EXPECT_DOUBLE_EQ(INPUT.ref_cell_factor,1.0); EXPECT_DOUBLE_EQ(INPUT.symmetry_prec,1.0e-5); - EXPECT_EQ(INPUT.cal_force,0); - EXPECT_TRUE(INPUT.dump_force); - EXPECT_TRUE(INPUT.dump_vel); - EXPECT_TRUE(INPUT.dump_virial); + EXPECT_EQ(INPUT.cal_force, 0); EXPECT_DOUBLE_EQ(INPUT.force_thr,1.0e-3); EXPECT_DOUBLE_EQ(INPUT.force_thr_ev2,0); EXPECT_DOUBLE_EQ(INPUT.stress_thr,1.0e-2); @@ -354,6 +351,9 @@ TEST_F(InputTest, Default) EXPECT_DOUBLE_EQ(INPUT.mdp.msst_vel,0); EXPECT_DOUBLE_EQ(INPUT.mdp.msst_vis,0); EXPECT_EQ(INPUT.mdp.pot_file,"graph.pb"); + EXPECT_TRUE(INPUT.mdp.dump_force); + EXPECT_TRUE(INPUT.mdp.dump_vel); + EXPECT_TRUE(INPUT.mdp.dump_virial); } TEST_F(InputTest, Read) @@ -416,10 +416,7 @@ TEST_F(InputTest, Read) EXPECT_EQ(INPUT.symmetry,1); EXPECT_FALSE(INPUT.init_vel); EXPECT_DOUBLE_EQ(INPUT.symmetry_prec,1.0e-5); - EXPECT_EQ(INPUT.cal_force,0); - EXPECT_FALSE(INPUT.dump_force); - EXPECT_FALSE(INPUT.dump_vel); - EXPECT_FALSE(INPUT.dump_virial); + EXPECT_EQ(INPUT.cal_force, 0); EXPECT_NEAR(INPUT.force_thr,1.0e-3,1.0e-7); EXPECT_DOUBLE_EQ(INPUT.force_thr_ev2,0); EXPECT_DOUBLE_EQ(INPUT.stress_thr,1.0e-2); @@ -694,6 +691,9 @@ TEST_F(InputTest, Read) EXPECT_DOUBLE_EQ(INPUT.mdp.msst_vel,0); EXPECT_DOUBLE_EQ(INPUT.mdp.msst_vis,0); EXPECT_EQ(INPUT.mdp.pot_file,"graph.pb"); + EXPECT_FALSE(INPUT.mdp.dump_force); + EXPECT_FALSE(INPUT.mdp.dump_vel); + EXPECT_FALSE(INPUT.mdp.dump_virial); } TEST_F(InputTest, Default_2) diff --git a/source/module_io/test/input_test_para.cpp b/source/module_io/test/input_test_para.cpp index 0a1904b194..3109eb5cee 100644 --- a/source/module_io/test/input_test_para.cpp +++ b/source/module_io/test/input_test_para.cpp @@ -89,10 +89,7 @@ TEST_F(InputParaTest,Bcast) EXPECT_EQ(INPUT.symmetry,0); EXPECT_FALSE(INPUT.init_vel); EXPECT_DOUBLE_EQ(INPUT.symmetry_prec,1.0e-5); - EXPECT_EQ(INPUT.cal_force,0); - EXPECT_TRUE(INPUT.dump_force); - EXPECT_TRUE(INPUT.dump_vel); - EXPECT_TRUE(INPUT.dump_virial); + EXPECT_EQ(INPUT.cal_force, 0); EXPECT_DOUBLE_EQ(INPUT.force_thr,1.0e-3); EXPECT_DOUBLE_EQ(INPUT.force_thr_ev2,0); EXPECT_DOUBLE_EQ(INPUT.stress_thr,1.0e-2); @@ -357,6 +354,9 @@ TEST_F(InputParaTest,Bcast) EXPECT_DOUBLE_EQ(INPUT.mdp.msst_vel,0); EXPECT_DOUBLE_EQ(INPUT.mdp.msst_vis,0); EXPECT_EQ(INPUT.mdp.pot_file,"graph.pb"); + EXPECT_TRUE(INPUT.mdp.dump_force); + EXPECT_TRUE(INPUT.mdp.dump_vel); + EXPECT_TRUE(INPUT.mdp.dump_virial); EXPECT_FALSE(INPUT.mixing_tau); EXPECT_FALSE(INPUT.mixing_dftu); EXPECT_EQ(INPUT.out_bandgap,0); diff --git a/source/module_io/write_input.cpp b/source/module_io/write_input.cpp index a62efb7a7f..3861616baa 100644 --- a/source/module_io/write_input.cpp +++ b/source/module_io/write_input.cpp @@ -276,9 +276,18 @@ ModuleBase::GlobalFunc::OUTP(ofs, "out_bandgap", out_bandgap, "if true, print ou ModuleBase::GlobalFunc::OUTP(ofs,"md_pfirst",mdp.md_pfirst,"initial target pressure"); ModuleBase::GlobalFunc::OUTP(ofs,"md_plast",mdp.md_plast,"final target pressure"); ModuleBase::GlobalFunc::OUTP(ofs,"md_pfreq",mdp.md_pfreq,"oscillation frequency, used to determine qmass of thermostats coupled with barostat"); - ModuleBase::GlobalFunc::OUTP(ofs, "dump_force", dump_force, "output atomic forces into the file MD_dump or not"); - ModuleBase::GlobalFunc::OUTP(ofs, "dump_vel", dump_vel, "output atomic velocities into the file MD_dump or not"); - ModuleBase::GlobalFunc::OUTP(ofs, "dump_virial", dump_virial, "output lattice virial into the file MD_dump or not"); + ModuleBase::GlobalFunc::OUTP(ofs, + "dump_force", + mdp.dump_force, + "output atomic forces into the file MD_dump or not"); + ModuleBase::GlobalFunc::OUTP(ofs, + "dump_vel", + mdp.dump_vel, + "output atomic velocities into the file MD_dump or not"); + ModuleBase::GlobalFunc::OUTP(ofs, + "dump_virial", + mdp.dump_virial, + "output lattice virial into the file MD_dump or not"); ofs << "\n#Parameters (10.Electric field and dipole correction)" << std::endl; ModuleBase::GlobalFunc::OUTP(ofs,"efield_flag",efield_flag,"add electric field"); diff --git a/source/module_md/fire.cpp b/source/module_md/fire.cpp index 631ce3a23b..539d8889e6 100644 --- a/source/module_md/fire.cpp +++ b/source/module_md/fire.cpp @@ -23,58 +23,61 @@ FIRE::~FIRE() { } -void FIRE::setup(ModuleESolver::ESolver *p_ensolve) +void FIRE::setup(ModuleESolver::ESolver *p_esolver, + const int &my_rank, + const std::string &global_readin_dir, + const double &force_thr) { ModuleBase::TITLE("FIRE", "setup"); ModuleBase::timer::tick("FIRE", "setup"); - MDrun::setup(p_ensolve); + MDrun::setup(p_esolver, my_rank, global_readin_dir); - check_force(); + check_force(force_thr); ModuleBase::timer::tick("FIRE", "setup"); } -void FIRE::first_half() +void FIRE::first_half(const int &my_rank) { ModuleBase::TITLE("FIRE", "first_half"); ModuleBase::timer::tick("FIRE", "first_half"); - MDrun::update_vel(force); + MDrun::update_vel(force, my_rank); check_FIRE(); - MDrun::update_pos(); + MDrun::update_pos(my_rank); ModuleBase::timer::tick("FIRE", "first_half"); } -void FIRE::second_half() +void FIRE::second_half(const int &my_rank, const double &force_thr) { ModuleBase::TITLE("FIRE", "second_half"); ModuleBase::timer::tick("FIRE", "second_half"); - MDrun::update_vel(force); + MDrun::update_vel(force, my_rank); - check_force(); + check_force(force_thr); ModuleBase::timer::tick("FIRE", "second_half"); } -void FIRE::outputMD(std::ofstream &ofs, bool cal_stress) +void FIRE::outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank) { - MDrun::outputMD(ofs, cal_stress); + MDrun::outputMD(ofs, cal_stress, my_rank); ofs << " LARGEST GRAD (eV/A) : " << max * ModuleBase::Hartree_to_eV * ModuleBase::ANGSTROM_AU << std::endl; std::cout << " LARGEST GRAD (eV/A) : " << max * ModuleBase::Hartree_to_eV * ModuleBase::ANGSTROM_AU << std::endl; } -void FIRE::write_restart() +void FIRE::write_restart(const int &my_rank, const std::string &global_out_dir) { - if (!GlobalV::MY_RANK) + if (!my_rank) { std::stringstream ssc; - ssc << GlobalV::global_out_dir << "Restart_md.dat"; + ssc << global_out_dir << "Restart_md.dat"; std::ofstream file(ssc.str().c_str()); file << step_ + step_rst_ << std::endl; @@ -89,14 +92,14 @@ void FIRE::write_restart() #endif } -void FIRE::restart() +void FIRE::restart(const int &my_rank, const std::string &global_readin_dir) { bool ok = true; - if (!GlobalV::MY_RANK) + if (!my_rank) { std::stringstream ssc; - ssc << GlobalV::global_readin_dir << "Restart_md.dat"; + ssc << global_readin_dir << "Restart_md.dat"; std::ifstream file(ssc.str().c_str()); if (!file) @@ -129,7 +132,7 @@ void FIRE::restart() #endif } -void FIRE::check_force() +void FIRE::check_force(const double &force_thr) { max = 0; @@ -144,7 +147,7 @@ void FIRE::check_force() } } - if (2.0 * max < GlobalV::FORCE_THR) + if (2.0 * max < force_thr) { stop = true; } diff --git a/source/module_md/fire.h b/source/module_md/fire.h index 836ac79344..de51aa48fd 100644 --- a/source/module_md/fire.h +++ b/source/module_md/fire.h @@ -9,13 +9,16 @@ class FIRE : public MDrun FIRE(MD_parameters &MD_para_in, UnitCell &unit_in); ~FIRE(); - void setup(ModuleESolver::ESolver *p_ensolve); - void first_half(); - void second_half(); - void outputMD(std::ofstream &ofs, bool cal_stress); - void write_restart(); - void restart(); - void check_force(); + void setup(ModuleESolver::ESolver *p_esolver, + const int &my_rank, + const std::string &global_readin_dir, + const double &force_thr); + void first_half(const int &my_rank); + void second_half(const int &my_rank, const double &force_thr); + void outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank); + void write_restart(const int &my_rank, const std::string &global_out_dir); + void restart(const int &my_rank, const std::string &global_readin_dir); + void check_force(const double &force_thr); void check_FIRE(); double max; // max force diff --git a/source/module_md/langevin.cpp b/source/module_md/langevin.cpp index eb06295039..a951d8e3b3 100644 --- a/source/module_md/langevin.cpp +++ b/source/module_md/langevin.cpp @@ -70,7 +70,7 @@ void Langevin::post_force() { if (GlobalV::MY_RANK == 0) { - double t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_tfirst, mdp.md_tlast); + double t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, mdp.md_tfirst, mdp.md_tlast); ModuleBase::Vector3 fictitious_force; for (int i = 0; i < ucell.nat; ++i) { diff --git a/source/module_md/md_base.cpp b/source/module_md/md_base.cpp index e3b706c20c..264966d9c3 100644 --- a/source/module_md/md_base.cpp +++ b/source/module_md/md_base.cpp @@ -44,11 +44,11 @@ MDrun::~MDrun() delete[] force; } -void MDrun::setup(ModuleESolver::ESolver *p_esolver) +void MDrun::setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const std::string &global_readin_dir) { if (mdp.md_restart) { - restart(); + restart(my_rank, global_readin_dir); } Print_Info::print_screen(0, 0, step_ + step_rst_); @@ -59,20 +59,20 @@ void MDrun::setup(ModuleESolver::ESolver *p_esolver) ucell.ionic_position_updated = true; } -void MDrun::first_half() +void MDrun::first_half(const int &my_rank) { - update_vel(force); - update_pos(); + update_vel(force, my_rank); + update_pos(my_rank); } -void MDrun::second_half() +void MDrun::second_half(const int &my_rank) { - update_vel(force); + update_vel(force, my_rank); } -void MDrun::update_pos() +void MDrun::update_pos(const int &my_rank) { - if (GlobalV::MY_RANK == 0) + if (my_rank == 0) { for (int i = 0; i < ucell.nat; ++i) { @@ -98,9 +98,9 @@ void MDrun::update_pos() ucell.update_pos_taud(pos); } -void MDrun::update_vel(const ModuleBase::Vector3 *force) +void MDrun::update_vel(const ModuleBase::Vector3 *force, const int &my_rank) { - if (GlobalV::MY_RANK == 0) + if (my_rank == 0) { for (int i = 0; i < ucell.nat; ++i) { @@ -119,9 +119,9 @@ void MDrun::update_vel(const ModuleBase::Vector3 *force) #endif } -void MDrun::outputMD(std::ofstream &ofs, bool cal_stress) +void MDrun::outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank) { - if (GlobalV::MY_RANK) + if (my_rank) return; t_current = MD_func::current_temp(kinetic, ucell.nat, frozen_freedom_, allmass, vel); @@ -183,12 +183,12 @@ void MDrun::outputMD(std::ofstream &ofs, bool cal_stress) ofs << std::endl; } -void MDrun::write_restart() +void MDrun::write_restart(const int &my_rank, const std::string &global_out_dir) { - if (!GlobalV::MY_RANK) + if (!my_rank) { std::stringstream ssc; - ssc << GlobalV::global_out_dir << "Restart_md.dat"; + ssc << global_out_dir << "Restart_md.dat"; std::ofstream file(ssc.str().c_str()); file << step_ + step_rst_ << std::endl; @@ -199,7 +199,7 @@ void MDrun::write_restart() #endif } -void MDrun::restart() +void MDrun::restart(const int &my_rank, const std::string &global_readin_dir) { - step_rst_ = MD_func::current_step(GlobalV::MY_RANK, GlobalV::global_readin_dir); + step_rst_ = MD_func::current_step(my_rank, global_readin_dir); } \ No newline at end of file diff --git a/source/module_md/md_base.h b/source/module_md/md_base.h index cfba171758..30c0a9b5f2 100644 --- a/source/module_md/md_base.h +++ b/source/module_md/md_base.h @@ -14,44 +14,44 @@ class MDrun * @brief init before running md, calculate energy, force, and stress of the initial configuration. * @param p_esolver the energy solver used in md */ - virtual void setup(ModuleESolver::ESolver *p_esolver); + virtual void setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const std::string &global_readin_dir); /** * @brief the first half of equation of motion, update velocities and positions */ - virtual void first_half(); + virtual void first_half(const int &my_rank); /** * @brief the second half of equation of motion, update velocities */ - virtual void second_half(); + virtual void second_half(const int &my_rank); /** * @brief perform one step update of pos due to atomic velocity */ - virtual void update_pos(); + virtual void update_pos(const int &my_rank); /** * @brief perform half-step update of vel due to atomic force */ - virtual void update_vel(const ModuleBase::Vector3 *force); + virtual void update_vel(const ModuleBase::Vector3 *force, const int &my_rank); /** * @brief output MD information such as energy, temperature, and pressure * @param ofs determine the output files * @param cal_stress whether calculate and output stress */ - virtual void outputMD(std::ofstream &ofs, bool cal_stress); + virtual void outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank); /** * @brief write the information into files used for MD restarting */ - virtual void write_restart(); + virtual void write_restart(const int &my_rank, const std::string &global_out_dir); /** * @brief restart MD when md_restart is true */ - virtual void restart(); + virtual void restart(const int &my_rank, const std::string &global_readin_dir); MD_parameters &mdp; UnitCell &ucell; diff --git a/source/module_md/md_func.cpp b/source/module_md/md_func.cpp index 848e767e2a..cebed7aaed 100644 --- a/source/module_md/md_func.cpp +++ b/source/module_md/md_func.cpp @@ -259,7 +259,7 @@ void outStress(const ModuleBase::matrix &virial, const ModuleBase::matrix &stres void MDdump(const int &step, const UnitCell &unit_in, - const Input &inp, + const MD_parameters &mdp, const ModuleBase::matrix &virial, const ModuleBase::Vector3 *force, const ModuleBase::Vector3 *vel) @@ -294,7 +294,7 @@ void MDdump(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 (GlobalV::CAL_STRESS && inp.dump_virial) + if (GlobalV::CAL_STRESS && mdp.dump_virial) { ofs << "VIRIAL (kbar)" << std::endl; for (int i = 0; i < 3; ++i) @@ -305,11 +305,11 @@ void MDdump(const int &step, } ofs << "INDEX LABEL POSITION (Angstrom)"; - if (inp.dump_force) + if (mdp.dump_force) { ofs << " FORCE (eV/Angstrom)"; } - if (inp.dump_vel) + if (mdp.dump_vel) { ofs << " VELOCITY (Angstrom/fs)"; } @@ -323,13 +323,13 @@ void MDdump(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 (inp.dump_force) + if (mdp.dump_force) { ofs << " " << force[index].x * unit_force << " " << force[index].y * unit_force << " " << force[index].z * unit_force; } - if (inp.dump_vel) + if (mdp.dump_vel) { ofs << " " << vel[index].x * unit_vel << " " << vel[index].y * unit_vel << " " << vel[index].z * unit_vel; @@ -371,9 +371,9 @@ void getMassMbl(const UnitCell &unit_in, } } -double target_temp(const int &istep, const double &tfirst, const double &tlast) +double target_temp(const int &istep, const int &nstep, const double &tfirst, const double &tlast) { - double delta = (double)(istep) / GlobalV::MD_NSTEP; + double delta = static_cast(istep) / nstep; return tfirst + delta * (tlast - tfirst); } diff --git a/source/module_md/md_func.h b/source/module_md/md_func.h index ca50fe44a5..5ecc6749f2 100644 --- a/source/module_md/md_func.h +++ b/source/module_md/md_func.h @@ -43,7 +43,7 @@ void outStress(const ModuleBase::matrix &virial, const ModuleBase::matrix &stres void MDdump(const int &step, const UnitCell &unit_in, - const Input &inp, + const MD_parameters &mdp, const ModuleBase::matrix &virial, const ModuleBase::Vector3 *force, const ModuleBase::Vector3 *vel); @@ -53,7 +53,7 @@ void getMassMbl(const UnitCell &unit_in, ModuleBase::Vector3 &frozen, ModuleBase::Vector3 *ionmbl); -double target_temp(const int &istep, const double &tfirst, const double &tlast); +double target_temp(const int &istep, const int &nstep, const double &tfirst, const double &tlast); double current_temp(double &kinetic, const int &natom, diff --git a/source/module_md/md_para.h b/source/module_md/md_para.h index c869dd07ac..0da75b4290 100644 --- a/source/module_md/md_para.h +++ b/source/module_md/md_para.h @@ -48,6 +48,10 @@ class MD_parameters md_tolerance = 100.0; md_nraise = 1; + + dump_force = true; + dump_vel = true; + dump_virial = true; }; ~MD_parameters(){}; @@ -92,6 +96,10 @@ class MD_parameters double md_tolerance; // tolerance for velocity rescaling (K) int md_nraise; // parameters used when md_type=nvt + + bool dump_force; // output atomic forces into the file MD_dump or not. liuyu 2023-03-01 + bool dump_vel; // output atomic velocities into the file MD_dump or not. liuyu 2023-03-01 + bool dump_virial; // output lattice virial into the file MD_dump or not. liuyu 2023-03-01 }; #endif \ No newline at end of file diff --git a/source/module_md/msst.cpp b/source/module_md/msst.cpp index 7ddeef108f..a00c82d2dd 100644 --- a/source/module_md/msst.cpp +++ b/source/module_md/msst.cpp @@ -31,12 +31,12 @@ MSST::~MSST() delete[] old_v; } -void MSST::setup(ModuleESolver::ESolver *p_esolver) +void MSST::setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const std::string &global_readin_dir) { ModuleBase::TITLE("MSST", "setup"); ModuleBase::timer::tick("MSST", "setup"); - MDrun::setup(p_esolver); + MDrun::setup(p_esolver, my_rank, global_readin_dir); ucell.cell_parameter_updated = true; int sd = mdp.msst_direction; @@ -69,7 +69,7 @@ void MSST::setup(ModuleESolver::ESolver *p_esolver) ModuleBase::timer::tick("MSST", "setup"); } -void MSST::first_half() +void MSST::first_half(const int &my_rank, std::ofstream &ofs) { ModuleBase::TITLE("MSST", "first_half"); ModuleBase::timer::tick("MSST", "first_half"); @@ -91,7 +91,7 @@ void MSST::first_half() } // propagate velocity sum 1/2 step by temporarily propagating the velocities - propagate_vel(); + propagate_vel(my_rank); vsum = vel_sum(); @@ -108,21 +108,21 @@ void MSST::first_half() vol = ucell.omega + omega[sd] * dthalf; // rescale positions and change box size - rescale(vol); + rescale(ofs, vol); // propagate atom positions 1 time step - MDrun::update_pos(); + MDrun::update_pos(my_rank); // propagate volume 1/2 step vol = ucell.omega + omega[sd] * dthalf; // rescale positions and change box size - rescale(vol); + rescale(ofs, vol); ModuleBase::timer::tick("MSST", "first_half"); } -void MSST::second_half() +void MSST::second_half(const int &my_rank) { ModuleBase::TITLE("MSST", "second_half"); ModuleBase::timer::tick("MSST", "second_half"); @@ -132,7 +132,7 @@ void MSST::second_half() energy_ = potential + kinetic; // propagate velocities 1/2 step - propagate_vel(); + propagate_vel(my_rank); vsum = vel_sum(); MD_func::compute_stress(ucell, vel, allmass, virial, stress); @@ -147,17 +147,17 @@ void MSST::second_half() ModuleBase::timer::tick("MSST", "second_half"); } -void MSST::outputMD(std::ofstream &ofs, bool cal_stress) +void MSST::outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank) { - MDrun::outputMD(ofs, cal_stress); + MDrun::outputMD(ofs, cal_stress, my_rank); } -void MSST::write_restart() +void MSST::write_restart(const int &my_rank, const std::string &global_out_dir) { - if (!GlobalV::MY_RANK) + if (!my_rank) { std::stringstream ssc; - ssc << GlobalV::global_out_dir << "Restart_md.dat"; + ssc << global_out_dir << "Restart_md.dat"; std::ofstream file(ssc.str().c_str()); file << step_ + step_rst_ << std::endl; @@ -174,14 +174,14 @@ void MSST::write_restart() #endif } -void MSST::restart() +void MSST::restart(const int &my_rank, const std::string &global_readin_dir) { bool ok = true; - if (!GlobalV::MY_RANK) + if (!my_rank) { std::stringstream ssc; - ssc << GlobalV::global_readin_dir << "Restart_md.dat"; + ssc << global_readin_dir << "Restart_md.dat"; std::ifstream file(ssc.str().c_str()); if (!file) @@ -232,7 +232,7 @@ double MSST::vel_sum() return vsum; } -void MSST::rescale(double volume) +void MSST::rescale(std::ofstream &ofs, const double &volume) { int sd = mdp.msst_direction; @@ -244,7 +244,7 @@ void MSST::rescale(double volume) // ucell.latvec.e22 *= 1.01; // ucell.latvec.e33 *= 1.01; - ucell.setup_cell_after_vc(GlobalV::ofs_running); + ucell.setup_cell_after_vc(ofs); // rescale velocity for (int i = 0; i < ucell.nat; ++i) @@ -256,33 +256,41 @@ void MSST::rescale(double volume) } } -void MSST::propagate_vel() +void MSST::propagate_vel(const int &my_rank) { - 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); - - for (int i = 0; i < ucell.nat; ++i) + if (my_rank == 0) { - ModuleBase::Vector3 const_C = force[i] / allmass[i]; - ModuleBase::Vector3 const_D; - const_D.set(fac / allmass[i], fac / allmass[i], fac / allmass[i]); - const_D[sd] -= 2 * omega[sd] / ucell.omega; + 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); - for (int k = 0; k < 3; ++k) + for (int i = 0; i < ucell.nat; ++i) { - if (fabs(dthalf * const_D[k]) > 1e-6) - { - double expd = exp(dthalf * const_D[k]); - vel[i][k] = expd * (const_C[k] + const_D[k] * vel[i][k] - const_C[k] / expd) / const_D[k]; - } - else + ModuleBase::Vector3 const_C = force[i] / allmass[i]; + ModuleBase::Vector3 const_D; + const_D.set(fac / allmass[i], fac / allmass[i], fac / allmass[i]); + const_D[sd] -= 2 * omega[sd] / ucell.omega; + + for (int k = 0; k < 3; ++k) { - vel[i][k] += (const_C[k] + const_D[k] * vel[i][k]) * dthalf - + 0.5 * (const_D[k] * const_D[k] * vel[i][k] + const_C[k] * const_D[k]) * dthalf * dthalf; + if (fabs(dthalf * const_D[k]) > 1e-6) + { + double expd = exp(dthalf * const_D[k]); + vel[i][k] = expd * (const_C[k] + const_D[k] * vel[i][k] - const_C[k] / expd) / const_D[k]; + } + else + { + vel[i][k] + += (const_C[k] + const_D[k] * vel[i][k]) * dthalf + + 0.5 * (const_D[k] * const_D[k] * vel[i][k] + const_C[k] * const_D[k]) * dthalf * dthalf; + } } } } + +#ifdef __MPI + MPI_Bcast(vel, ucell.nat * 3, MPI_DOUBLE, 0, MPI_COMM_WORLD); +#endif } void MSST::propagate_voldot() diff --git a/source/module_md/msst.h b/source/module_md/msst.h index 808b34f7d6..b5fb29066e 100644 --- a/source/module_md/msst.h +++ b/source/module_md/msst.h @@ -10,15 +10,15 @@ class MSST : public MDrun ~MSST(); void setup(ModuleESolver::ESolver *p_ensolve); - void first_half(); - void second_half(); - void outputMD(std::ofstream &ofs, bool cal_stress); - void write_restart(); - void restart(); + void first_half(const int &my_rank, std::ofstream &ofs); + void second_half(const int &my_rank); + void outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank); + void write_restart(const int &my_rank, const std::string &global_out_dir); + void restart(const int &my_rank, const std::string &global_readin_dir); double extra_term(); double vel_sum(); - void rescale(double volume); - void propagate_vel(); + void rescale(std::ofstream &ofs, const double &volume); + void propagate_vel(const int &my_rank); void propagate_voldot(); ModuleBase::Vector3 *old_v; diff --git a/source/module_md/nhchain.cpp b/source/module_md/nhchain.cpp index c357888124..98ca2adf89 100644 --- a/source/module_md/nhchain.cpp +++ b/source/module_md/nhchain.cpp @@ -147,19 +147,19 @@ Nose_Hoover::~Nose_Hoover() } } -void Nose_Hoover::setup(ModuleESolver::ESolver *p_ensolve) +void Nose_Hoover::setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const std::string &global_readin_dir) { ModuleBase::TITLE("Nose_Hoover", "setup"); ModuleBase::timer::tick("Nose_Hoover", "setup"); - MDrun::setup(p_ensolve); + MDrun::setup(p_esolver, my_rank, global_readin_dir); if (mdp.md_type == "npt") { ucell.cell_parameter_updated = true; } // determine target temperature - t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_tfirst, mdp.md_tlast); + t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, mdp.md_tfirst, mdp.md_tlast); // init thermostats coupled with particles mass_eta[0] = tdof * t_target / mdp.md_tfreq / mdp.md_tfreq; @@ -204,7 +204,7 @@ void Nose_Hoover::setup(ModuleESolver::ESolver *p_ensolve) ModuleBase::timer::tick("Nose_Hoover", "setup"); } -void Nose_Hoover::first_half() +void Nose_Hoover::first_half(const int &my_rank, std::ofstream &ofs) { ModuleBase::TITLE("Nose_Hoover", "first_half"); ModuleBase::timer::tick("Nose_Hoover", "first_half"); @@ -216,7 +216,7 @@ void Nose_Hoover::first_half() } // update target T - t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_tfirst, mdp.md_tlast); + t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, mdp.md_tfirst, mdp.md_tlast); // update thermostats coupled with particles particle_thermo(); @@ -241,33 +241,33 @@ void Nose_Hoover::first_half() } // perform half-step update of vel due to atomic force - MDrun::update_vel(force); + MDrun::update_vel(force, my_rank); if (npt_flag) { // perform half-step update of volume - update_volume(); + update_volume(ofs); } // perform one step update of pos due to atomic velocity - MDrun::update_pos(); + MDrun::update_pos(my_rank); if (npt_flag) { // perform half-step update of volume - update_volume(); + update_volume(ofs); } ModuleBase::timer::tick("Nose_Hoover", "first_half"); } -void Nose_Hoover::second_half() +void Nose_Hoover::second_half(const int &my_rank) { ModuleBase::TITLE("Nose_Hoover", "second_half"); ModuleBase::timer::tick("Nose_Hoover", "second_half"); // perform half-step update of vel due to atomic force - MDrun::update_vel(force); + MDrun::update_vel(force, my_rank); if (npt_flag) { @@ -302,17 +302,17 @@ void Nose_Hoover::second_half() ModuleBase::timer::tick("Nose_Hoover", "second_half"); } -void Nose_Hoover::outputMD(std::ofstream &ofs, bool cal_stress) +void Nose_Hoover::outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank) { - MDrun::outputMD(ofs, cal_stress); + MDrun::outputMD(ofs, cal_stress, my_rank); } -void Nose_Hoover::write_restart() +void Nose_Hoover::write_restart(const int &my_rank, const std::string &global_out_dir) { - if (!GlobalV::MY_RANK) + if (!my_rank) { std::stringstream ssc; - ssc << GlobalV::global_out_dir << "Restart_md.dat"; + ssc << global_out_dir << "Restart_md.dat"; std::ofstream file(ssc.str().c_str()); file << step_ + step_rst_ << std::endl; @@ -356,16 +356,16 @@ void Nose_Hoover::write_restart() #endif } -void Nose_Hoover::restart() +void Nose_Hoover::restart(const int &my_rank, const std::string &global_readin_dir) { bool ok = true; bool ok2 = true; bool ok3 = true; - if (!GlobalV::MY_RANK) + if (!my_rank) { std::stringstream ssc; - ssc << GlobalV::global_readin_dir << "Restart_md.dat"; + ssc << global_readin_dir << "Restart_md.dat"; std::ifstream file(ssc.str().c_str()); if (!file) @@ -703,7 +703,7 @@ void Nose_Hoover::vel_baro() } } -void Nose_Hoover::update_volume() +void Nose_Hoover::update_volume(std::ofstream &ofs) { double factor; @@ -801,12 +801,12 @@ void Nose_Hoover::update_volume() } // reset ucell and pos due to change of lattice - ucell.setup_cell_after_vc(GlobalV::ofs_running); + ucell.setup_cell_after_vc(ofs); } void Nose_Hoover::target_stress() { - double delta = (double)(step_ + step_rst_) / GlobalV::MD_NSTEP; + double delta = static_cast(step_ + step_rst_) / mdp.md_nstep; p_hydro = 0; for (int i = 0; i < 3; ++i) diff --git a/source/module_md/nhchain.h b/source/module_md/nhchain.h index f0231c7737..7662ad8ba5 100644 --- a/source/module_md/nhchain.h +++ b/source/module_md/nhchain.h @@ -9,12 +9,12 @@ class Nose_Hoover : public MDrun Nose_Hoover(MD_parameters &MD_para_in, UnitCell &unit_in); ~Nose_Hoover(); - void setup(ModuleESolver::ESolver *p_ensolve); - void first_half(); - void second_half(); - void outputMD(std::ofstream &ofs, bool cal_stress); - void write_restart(); - void restart(); + void setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const std::string &global_readin_dir); + void first_half(const int &my_rank, std::ofstream &ofs); + void second_half(const int &my_rank); + void outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank); + void write_restart(const int &my_rank, const std::string &global_out_dir); + void restart(const int &my_rank, const std::string &global_readin_dir); // perform half-step update of thermostats coupled with particles void particle_thermo(); @@ -35,7 +35,7 @@ class Nose_Hoover : public MDrun void couple_stress(); // perform half-step update of volume - void update_volume(); + void update_volume(std::ofstream &ofs); const int nc_tchain = 1; const int nc_pchain = 1; diff --git a/source/module_md/run_md.cpp b/source/module_md/run_md.cpp index f7b1d91dd2..3b90c3d953 100644 --- a/source/module_md/run_md.cpp +++ b/source/module_md/run_md.cpp @@ -17,32 +17,32 @@ Run_MD::~Run_MD() { } -void Run_MD::md_line(UnitCell &unit_in, ModuleESolver::ESolver *p_esolver) +void Run_MD::md_line(UnitCell &unit_in, ModuleESolver::ESolver *p_esolver, MD_parameters &md_para) { ModuleBase::TITLE("Run_MD", "md_line"); ModuleBase::timer::tick("Run_MD", "md_line"); // determine the md_type MDrun *mdrun; - if (INPUT.mdp.md_type == "fire") + if (md_para.md_type == "fire") { - mdrun = new FIRE(INPUT.mdp, unit_in); + mdrun = new FIRE(md_para, unit_in); } - else if ((INPUT.mdp.md_type == "nvt" && INPUT.mdp.md_thermostat == "nhc") || INPUT.mdp.md_type == "npt") + else if ((md_para.md_type == "nvt" && md_para.md_thermostat == "nhc") || md_para.md_type == "npt") { - mdrun = new Nose_Hoover(INPUT.mdp, unit_in); + mdrun = new Nose_Hoover(md_para, unit_in); } - else if (INPUT.mdp.md_type == "nve" || INPUT.mdp.md_type == "nvt") + else if (md_para.md_type == "nve" || md_para.md_type == "nvt") { - mdrun = new Verlet(INPUT.mdp, unit_in); + mdrun = new Verlet(md_para, unit_in); } - else if (INPUT.mdp.md_type == "langevin") + else if (md_para.md_type == "langevin") { - mdrun = new Langevin(INPUT.mdp, unit_in); + mdrun = new Langevin(md_para, unit_in); } - else if (INPUT.mdp.md_type == "msst") + else if (md_para.md_type == "msst") { - mdrun = new MSST(INPUT.mdp, unit_in); + mdrun = new MSST(md_para, unit_in); } else { @@ -50,7 +50,7 @@ void Run_MD::md_line(UnitCell &unit_in, ModuleESolver::ESolver *p_esolver) } // md cycle - while ((mdrun->step_ + mdrun->step_rst_) <= GlobalV::MD_NSTEP && !mdrun->stop) + while ((mdrun->step_ + mdrun->step_rst_) <= md_para.md_nstep && !mdrun->stop) { if (mdrun->step_ == 0) { @@ -59,12 +59,12 @@ void Run_MD::md_line(UnitCell &unit_in, ModuleESolver::ESolver *p_esolver) else { Print_Info::print_screen(0, 0, mdrun->step_ + mdrun->step_rst_); - mdrun->first_half(); + mdrun->first_half(GlobalV::MY_RANK); // update force and virial due to the update of atom positions MD_func::force_virial(p_esolver, mdrun->step_, mdrun->ucell, mdrun->potential, mdrun->force, mdrun->virial); - mdrun->second_half(); + mdrun->second_half(GlobalV::MY_RANK); MD_func::compute_stress(mdrun->ucell, mdrun->vel, mdrun->allmass, mdrun->virial, mdrun->stress); mdrun->t_current = MD_func::current_temp(mdrun->kinetic, @@ -76,11 +76,11 @@ void Run_MD::md_line(UnitCell &unit_in, ModuleESolver::ESolver *p_esolver) if ((mdrun->step_ + mdrun->step_rst_) % mdrun->mdp.md_dumpfreq == 0) { - mdrun->outputMD(GlobalV::ofs_running, GlobalV::CAL_STRESS); + mdrun->outputMD(GlobalV::ofs_running, GlobalV::CAL_STRESS, GlobalV::MY_RANK); MD_func::MDdump(mdrun->step_ + mdrun->step_rst_, mdrun->ucell, - INPUT, + md_para, mdrun->virial, mdrun->force, mdrun->vel); @@ -92,7 +92,7 @@ void Run_MD::md_line(UnitCell &unit_in, ModuleESolver::ESolver *p_esolver) std::stringstream file; file << GlobalV::global_stru_dir << "STRU_MD_" << mdrun->step_ + mdrun->step_rst_; mdrun->ucell.print_stru_file(file.str(), 1, 1); - mdrun->write_restart(); + mdrun->write_restart(GlobalV::MY_RANK, GlobalV::global_out_dir); } mdrun->step_++; diff --git a/source/module_md/run_md.h b/source/module_md/run_md.h index 9d80efffde..248041ddc6 100644 --- a/source/module_md/run_md.h +++ b/source/module_md/run_md.h @@ -1,6 +1,7 @@ #ifndef RUN_MD_CLASSIC_H #define RUN_MD_CLASSIC_H +#include "md_para.h" #include "module_esolver/esolver.h" class Run_MD @@ -9,7 +10,7 @@ class Run_MD Run_MD(); ~Run_MD(); - void md_line(UnitCell &unit_in, ModuleESolver::ESolver *p_esolver); + void md_line(UnitCell &unit_in, ModuleESolver::ESolver *p_esolver, MD_parameters &md_para); }; #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 848a9c5ff3..84ef39ffff 100644 --- a/source/module_md/test/md_func_test.cpp +++ b/source/module_md/test/md_func_test.cpp @@ -160,7 +160,7 @@ TEST_F(MD_func_test, compute_stress) TEST_F(MD_func_test, MDdump) { - MD_func::MDdump(0, ucell, INPUT, virial, force, vel); + MD_func::MDdump(0, ucell, INPUT.mdp, virial, force, vel); std::ifstream ifs("MD_dump"); std::string output_str; getline(ifs, output_str); diff --git a/source/module_md/test/setcell.h b/source/module_md/test/setcell.h index 4aa51ce4c9..9d32f1b9b8 100644 --- a/source/module_md/test/setcell.h +++ b/source/module_md/test/setcell.h @@ -122,9 +122,9 @@ class Setcell GlobalV::SEARCH_RADIUS = 8.5 * ModuleBase::ANGSTROM_AU; GlobalV::CAL_STRESS = 1; - INPUT.dump_virial = true; - INPUT.dump_force = true; - INPUT.dump_vel = true; + INPUT.mdp.dump_virial = true; + INPUT.mdp.dump_force = true; + INPUT.mdp.dump_vel = true; INPUT.mdp.md_restart = 0; INPUT.mdp.md_dt = 1; diff --git a/source/module_md/verlet.cpp b/source/module_md/verlet.cpp index 89d840fb33..9c60db7172 100644 --- a/source/module_md/verlet.cpp +++ b/source/module_md/verlet.cpp @@ -11,7 +11,7 @@ Verlet::~Verlet() { } -void Verlet::setup(ModuleESolver::ESolver *p_ensolve) +void Verlet::setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const std::string &global_readin_dir) { ModuleBase::TITLE("Verlet", "setup"); ModuleBase::timer::tick("Verlet", "setup"); @@ -21,29 +21,29 @@ void Verlet::setup(ModuleESolver::ESolver *p_ensolve) ModuleBase::timer::tick("Verlet", "setup"); } -void Verlet::first_half() +void Verlet::first_half(const int &my_rank) { ModuleBase::TITLE("Verlet", "first_half"); ModuleBase::timer::tick("Verlet", "first_half"); - MDrun::update_vel(force); - MDrun::update_pos(); + MDrun::update_vel(force, my_rank); + MDrun::update_pos(my_rank); ModuleBase::timer::tick("Verlet", "first_half"); } -void Verlet::second_half() +void Verlet::second_half(const int &my_rank) { ModuleBase::TITLE("Verlet", "second_half"); ModuleBase::timer::tick("Verlet", "second_half"); - MDrun::update_vel(force); - apply_thermostat(); + MDrun::update_vel(force, my_rank); + apply_thermostat(my_rank); ModuleBase::timer::tick("Verlet", "second_half"); } -void Verlet::apply_thermostat() +void Verlet::apply_thermostat(const int &my_rank) { double t_target = 0; t_current = MD_func::current_temp(kinetic, ucell.nat, frozen_freedom_, allmass, vel); @@ -53,7 +53,7 @@ void Verlet::apply_thermostat() } else if (mdp.md_thermostat == "rescaling") { - t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_tfirst, mdp.md_tlast); + t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, mdp.md_tfirst, mdp.md_tlast); if (abs(t_target - t_current) * ModuleBase::Hartree_to_K > mdp.md_tolerance) { thermalize(0, t_current, t_target); @@ -63,13 +63,13 @@ void Verlet::apply_thermostat() { if ((step_ + step_rst_) % mdp.md_nraise == 0) { - t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_tfirst, mdp.md_tlast); + t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, mdp.md_tfirst, mdp.md_tlast); thermalize(0, t_current, t_target); } } else if (mdp.md_thermostat == "anderson") { - if (GlobalV::MY_RANK == 0) + if (my_rank == 0) { double deviation; for (int i = 0; i < ucell.nat; ++i) @@ -93,7 +93,7 @@ void Verlet::apply_thermostat() } else if (mdp.md_thermostat == "berendsen") { - t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_tfirst, mdp.md_tlast); + t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, mdp.md_tfirst, mdp.md_tlast); thermalize(mdp.md_nraise, t_current, t_target); } else @@ -120,17 +120,17 @@ void Verlet::thermalize(const int &nraise, const double ¤t_temp, const dou } } -void Verlet::outputMD(std::ofstream &ofs, bool cal_stress) +void Verlet::outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank) { - MDrun::outputMD(ofs, cal_stress); + MDrun::outputMD(ofs, cal_stress, my_rank); } -void Verlet::write_restart() +void Verlet::write_restart(const int &my_rank, const std::string &global_out_dir) { - MDrun::write_restart(); + MDrun::write_restart(my_rank, global_out_dir); } -void Verlet::restart() +void Verlet::restart(const int &my_rank, const std::string &global_readin_dir) { - MDrun::restart(); + MDrun::restart(my_rank, global_readin_dir); } \ No newline at end of file diff --git a/source/module_md/verlet.h b/source/module_md/verlet.h index e3db755758..5a4f1f6f61 100644 --- a/source/module_md/verlet.h +++ b/source/module_md/verlet.h @@ -9,14 +9,14 @@ class Verlet : public MDrun Verlet(MD_parameters &MD_para_in, UnitCell &unit_in); ~Verlet(); - void setup(ModuleESolver::ESolver *p_ensolve); - void first_half(); - void second_half(); - void apply_thermostat(); + void setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const std::string &global_readin_dir); + void first_half(const int &my_rank); + void second_half(const int &my_rank); + void apply_thermostat(const int &my_rank); void thermalize(const int &nraise, const double ¤t_temp, const double &target_temp); - void outputMD(std::ofstream &ofs, bool cal_stress); - void write_restart(); - void restart(); + void outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank); + void write_restart(const int &my_rank, const std::string &global_out_dir); + void restart(const int &my_rank, const std::string &global_readin_dir); }; #endif \ No newline at end of file From 5d7e057f3184767f53833669bae86b7c9fd8131c Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Sun, 7 May 2023 11:07:12 +0800 Subject: [PATCH 3/5] Refactor: remove GlobalV in module_md --- source/module_io/input_conv.cpp | 1 + source/module_md/fire.cpp | 17 ++++++-------- source/module_md/fire.h | 11 ++++----- source/module_md/langevin.cpp | 26 ++++++++++----------- source/module_md/langevin.h | 12 +++++----- source/module_md/md_base.cpp | 2 +- source/module_md/md_base.h | 15 ++++++++++++- source/module_md/md_para.h | 4 ++++ source/module_md/msst.cpp | 2 +- source/module_md/msst.h | 2 +- source/module_md/run_md.cpp | 4 ++-- source/module_md/test/fire_test.cpp | 28 +++++++++++------------ source/module_md/test/langevin_test.cpp | 14 ++++++------ source/module_md/test/md_func_test.cpp | 2 +- source/module_md/test/msst_test.cpp | 14 ++++++------ source/module_md/test/nhchain_test.cpp | 18 +++++++-------- source/module_md/test/verlet_test.cpp | 30 ++++++++++++------------- source/module_md/verlet.cpp | 4 ++-- source/module_md/verlet.h | 2 +- 19 files changed, 110 insertions(+), 98 deletions(-) diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index b7bfd95745..f0ef86d954 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -204,6 +204,7 @@ void Input_Conv::Convert(void) GlobalV::CAL_FORCE = INPUT.cal_force; GlobalV::FORCE_THR = INPUT.force_thr; + INPUT.mdp.force_thr = INPUT.force_thr; GlobalV::STRESS_THR = INPUT.stress_thr; GlobalV::PRESS1 = INPUT.press1; diff --git a/source/module_md/fire.cpp b/source/module_md/fire.cpp index 539d8889e6..62b41fb50a 100644 --- a/source/module_md/fire.cpp +++ b/source/module_md/fire.cpp @@ -23,22 +23,19 @@ FIRE::~FIRE() { } -void FIRE::setup(ModuleESolver::ESolver *p_esolver, - const int &my_rank, - const std::string &global_readin_dir, - const double &force_thr) +void FIRE::setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const std::string &global_readin_dir) { ModuleBase::TITLE("FIRE", "setup"); ModuleBase::timer::tick("FIRE", "setup"); MDrun::setup(p_esolver, my_rank, global_readin_dir); - check_force(force_thr); + check_force(); ModuleBase::timer::tick("FIRE", "setup"); } -void FIRE::first_half(const int &my_rank) +void FIRE::first_half(const int &my_rank, std::ofstream &ofs) { ModuleBase::TITLE("FIRE", "first_half"); ModuleBase::timer::tick("FIRE", "first_half"); @@ -52,14 +49,14 @@ void FIRE::first_half(const int &my_rank) ModuleBase::timer::tick("FIRE", "first_half"); } -void FIRE::second_half(const int &my_rank, const double &force_thr) +void FIRE::second_half(const int &my_rank) { ModuleBase::TITLE("FIRE", "second_half"); ModuleBase::timer::tick("FIRE", "second_half"); MDrun::update_vel(force, my_rank); - check_force(force_thr); + check_force(); ModuleBase::timer::tick("FIRE", "second_half"); } @@ -132,7 +129,7 @@ void FIRE::restart(const int &my_rank, const std::string &global_readin_dir) #endif } -void FIRE::check_force(const double &force_thr) +void FIRE::check_force() { max = 0; @@ -147,7 +144,7 @@ void FIRE::check_force(const double &force_thr) } } - if (2.0 * max < force_thr) + if (2.0 * max < mdp.force_thr) { stop = true; } diff --git a/source/module_md/fire.h b/source/module_md/fire.h index de51aa48fd..7ab02c33a4 100644 --- a/source/module_md/fire.h +++ b/source/module_md/fire.h @@ -9,16 +9,13 @@ class FIRE : public MDrun FIRE(MD_parameters &MD_para_in, UnitCell &unit_in); ~FIRE(); - void setup(ModuleESolver::ESolver *p_esolver, - const int &my_rank, - const std::string &global_readin_dir, - const double &force_thr); - void first_half(const int &my_rank); - void second_half(const int &my_rank, const double &force_thr); + void setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const std::string &global_readin_dir); + void first_half(const int &my_rank, std::ofstream &ofs); + void second_half(const int &my_rank); void outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank); void write_restart(const int &my_rank, const std::string &global_out_dir); void restart(const int &my_rank, const std::string &global_readin_dir); - void check_force(const double &force_thr); + void check_force(); void check_FIRE(); double max; // max force diff --git a/source/module_md/langevin.cpp b/source/module_md/langevin.cpp index a951d8e3b3..55fa052be9 100644 --- a/source/module_md/langevin.cpp +++ b/source/module_md/langevin.cpp @@ -17,53 +17,53 @@ Langevin::~Langevin() delete[] total_force; } -void Langevin::setup(ModuleESolver::ESolver *p_ensolve) +void Langevin::setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const std::string &global_readin_dir) { ModuleBase::TITLE("Langevin", "setup"); ModuleBase::timer::tick("Langevin", "setup"); - MDrun::setup(p_ensolve); + MDrun::setup(p_esolver, my_rank, global_readin_dir); post_force(); ModuleBase::timer::tick("Langevin", "setup"); } -void Langevin::first_half() +void Langevin::first_half(const int &my_rank, std::ofstream &ofs) { ModuleBase::TITLE("Langevin", "first_half"); ModuleBase::timer::tick("Langevin", "first_half"); - MDrun::update_vel(total_force); - MDrun::update_pos(); + MDrun::update_vel(total_force, my_rank); + MDrun::update_pos(my_rank); ModuleBase::timer::tick("Langevin", "first_half"); } -void Langevin::second_half() +void Langevin::second_half(const int &my_rank) { ModuleBase::TITLE("Langevin", "second_half"); ModuleBase::timer::tick("Langevin", "second_half"); post_force(); - MDrun::update_vel(total_force); + MDrun::update_vel(total_force, my_rank); ModuleBase::timer::tick("Langevin", "second_half"); } -void Langevin::outputMD(std::ofstream &ofs, bool cal_stress) +void Langevin::outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank) { - MDrun::outputMD(ofs, cal_stress); + MDrun::outputMD(ofs, cal_stress, my_rank); } -void Langevin::write_restart() +void Langevin::write_restart(const int &my_rank, const std::string &global_out_dir) { - MDrun::write_restart(); + MDrun::write_restart(my_rank, global_out_dir); } -void Langevin::restart() +void Langevin::restart(const int &my_rank, const std::string &global_readin_dir) { - MDrun::restart(); + MDrun::restart(my_rank, global_readin_dir); } void Langevin::post_force() diff --git a/source/module_md/langevin.h b/source/module_md/langevin.h index d24c5822f3..44ea9a6fba 100644 --- a/source/module_md/langevin.h +++ b/source/module_md/langevin.h @@ -9,12 +9,12 @@ class Langevin : public MDrun Langevin(MD_parameters &MD_para_in, UnitCell &unit_in); ~Langevin(); - void setup(ModuleESolver::ESolver *p_ensolve); - void first_half(); - void second_half(); - void outputMD(std::ofstream &ofs, bool cal_stress); - void write_restart(); - void restart(); + void setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const std::string &global_readin_dir); + void first_half(const int &my_rank, std::ofstream &ofs); + void second_half(const int &my_rank); + void outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank); + void write_restart(const int &my_rank, const std::string &global_out_dir); + void restart(const int &my_rank, const std::string &global_readin_dir); void post_force(); ModuleBase::Vector3 *total_force; // total force = true force + Langevin fictitious_force diff --git a/source/module_md/md_base.cpp b/source/module_md/md_base.cpp index 264966d9c3..0449b871f8 100644 --- a/source/module_md/md_base.cpp +++ b/source/module_md/md_base.cpp @@ -59,7 +59,7 @@ void MDrun::setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const s ucell.ionic_position_updated = true; } -void MDrun::first_half(const int &my_rank) +void MDrun::first_half(const int &my_rank, std::ofstream &ofs) { update_vel(force, my_rank); update_pos(my_rank); diff --git a/source/module_md/md_base.h b/source/module_md/md_base.h index 30c0a9b5f2..39c28d1283 100644 --- a/source/module_md/md_base.h +++ b/source/module_md/md_base.h @@ -13,26 +13,34 @@ class MDrun /** * @brief init before running md, calculate energy, force, and stress of the initial configuration. * @param p_esolver the energy solver used in md + * @param my_rank MPI rank of the processer + * @param global_readin_dir directory of files for reading */ virtual void setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const std::string &global_readin_dir); /** * @brief the first half of equation of motion, update velocities and positions + * @param my_rank MPI rank of the processer + * @param ofs determine the output files */ - virtual void first_half(const int &my_rank); + virtual void first_half(const int &my_rank, std::ofstream &ofs); /** * @brief the second half of equation of motion, update velocities + * @param my_rank MPI rank of the processer */ virtual void second_half(const int &my_rank); /** * @brief perform one step update of pos due to atomic velocity + * @param my_rank MPI rank of the processer */ virtual void update_pos(const int &my_rank); /** * @brief perform half-step update of vel due to atomic force + * @param force atomic forces + * @param my_rank MPI rank of the processer */ virtual void update_vel(const ModuleBase::Vector3 *force, const int &my_rank); @@ -40,16 +48,21 @@ class MDrun * @brief output MD information such as energy, temperature, and pressure * @param ofs determine the output files * @param cal_stress whether calculate and output stress + * @param my_rank MPI rank of the processer */ virtual void outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank); /** * @brief write the information into files used for MD restarting + * @param my_rank MPI rank of the processer + * @param global_out_dir directory of output files */ virtual void write_restart(const int &my_rank, const std::string &global_out_dir); /** * @brief restart MD when md_restart is true + * @param my_rank MPI rank of the processer + * @param global_readin_dir directory of files for reading */ virtual void restart(const int &my_rank, const std::string &global_readin_dir); diff --git a/source/module_md/md_para.h b/source/module_md/md_para.h index 0da75b4290..79d26aa8f8 100644 --- a/source/module_md/md_para.h +++ b/source/module_md/md_para.h @@ -52,6 +52,8 @@ class MD_parameters dump_force = true; dump_vel = true; dump_virial = true; + + force_thr = 1.0e-3; }; ~MD_parameters(){}; @@ -100,6 +102,8 @@ class MD_parameters bool dump_force; // output atomic forces into the file MD_dump or not. liuyu 2023-03-01 bool dump_vel; // output atomic velocities into the file MD_dump or not. liuyu 2023-03-01 bool dump_virial; // output lattice virial into the file MD_dump or not. liuyu 2023-03-01 + + double force_thr; }; #endif \ No newline at end of file diff --git a/source/module_md/msst.cpp b/source/module_md/msst.cpp index a00c82d2dd..5f5bdf00ed 100644 --- a/source/module_md/msst.cpp +++ b/source/module_md/msst.cpp @@ -102,7 +102,7 @@ void MSST::first_half(const int &my_rank, std::ofstream &ofs) } // propagate velocities 1/2 step using the new velocity sum - propagate_vel(); + propagate_vel(my_rank); // propagate volume 1/2 step vol = ucell.omega + omega[sd] * dthalf; diff --git a/source/module_md/msst.h b/source/module_md/msst.h index b5fb29066e..1d4a342b7c 100644 --- a/source/module_md/msst.h +++ b/source/module_md/msst.h @@ -9,7 +9,7 @@ class MSST : public MDrun MSST(MD_parameters &MD_para_in, UnitCell &unit_in); ~MSST(); - void setup(ModuleESolver::ESolver *p_ensolve); + void setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const std::string &global_readin_dir); void first_half(const int &my_rank, std::ofstream &ofs); void second_half(const int &my_rank); void outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank); diff --git a/source/module_md/run_md.cpp b/source/module_md/run_md.cpp index 3b90c3d953..58f383844f 100644 --- a/source/module_md/run_md.cpp +++ b/source/module_md/run_md.cpp @@ -54,12 +54,12 @@ void Run_MD::md_line(UnitCell &unit_in, ModuleESolver::ESolver *p_esolver, MD_pa { if (mdrun->step_ == 0) { - mdrun->setup(p_esolver); + mdrun->setup(p_esolver, GlobalV::MY_RANK, GlobalV::global_readin_dir); } else { Print_Info::print_screen(0, 0, mdrun->step_ + mdrun->step_rst_); - mdrun->first_half(GlobalV::MY_RANK); + mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); // update force and virial due to the update of atom positions MD_func::force_virial(p_esolver, mdrun->step_, mdrun->ucell, mdrun->potential, mdrun->force, mdrun->virial); diff --git a/source/module_md/test/fire_test.cpp b/source/module_md/test/fire_test.cpp index a8219e5ec1..c3808affb8 100644 --- a/source/module_md/test/fire_test.cpp +++ b/source/module_md/test/fire_test.cpp @@ -32,7 +32,7 @@ * - output MD information such as energy, temperature, and pressure */ -class FIRE_test : public testing::Test +class FIREtest : public testing::Test { protected: MDrun *mdrun; @@ -47,7 +47,7 @@ class FIRE_test : public testing::Test p_esolver->Init(INPUT, ucell); mdrun = new FIRE(INPUT.mdp, ucell); - mdrun->setup(p_esolver); + mdrun->setup(p_esolver, GlobalV::MY_RANK, GlobalV::global_readin_dir); } void TearDown() @@ -56,7 +56,7 @@ class FIRE_test : public testing::Test } }; -TEST_F(FIRE_test, setup) +TEST_F(FIREtest, Setup) { EXPECT_NEAR(mdrun->t_current * ModuleBase::Hartree_to_K, 299.99999999999665, doublethreshold); EXPECT_NEAR(mdrun->stress(0, 0), 6.0100555286436806e-06, doublethreshold); @@ -70,9 +70,9 @@ TEST_F(FIRE_test, setup) EXPECT_NEAR(mdrun->stress(2, 2), 1.6060561926126463e-06, doublethreshold); } -TEST_F(FIRE_test, first_half) +TEST_F(FIREtest, FirstHalf) { - mdrun->first_half(); + mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); EXPECT_NEAR(mdrun->pos[0].x, -0.00045447059554315662, doublethreshold); EXPECT_NEAR(mdrun->pos[0].y, 0.00032646833232493271, doublethreshold); @@ -101,10 +101,10 @@ TEST_F(FIRE_test, first_half) EXPECT_NEAR(mdrun->vel[3].z, -2.5498181033639999e-05, doublethreshold); } -TEST_F(FIRE_test, second_half) +TEST_F(FIREtest, SecondHalf) { - mdrun->first_half(); - mdrun->second_half(); + mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); + mdrun->second_half(GlobalV::MY_RANK); EXPECT_NEAR(mdrun->pos[0].x, -0.00045447059554315662, doublethreshold); EXPECT_NEAR(mdrun->pos[0].y, 0.00032646833232493271, doublethreshold); @@ -133,11 +133,11 @@ TEST_F(FIRE_test, second_half) EXPECT_NEAR(mdrun->vel[3].z, -2.5498181033639999e-05, doublethreshold); } -TEST_F(FIRE_test, write_restart) +TEST_F(FIREtest, WriteRestart) { mdrun->step_ = 1; mdrun->step_rst_ = 2; - mdrun->write_restart(); + mdrun->write_restart(GlobalV::MY_RANK, GlobalV::global_out_dir); std::ifstream ifs("Restart_md.dat"); std::string output_str; @@ -154,9 +154,9 @@ TEST_F(FIRE_test, write_restart) ifs.close(); } -TEST_F(FIRE_test, restart) +TEST_F(FIREtest, Restart) { - mdrun->restart(); + mdrun->restart(GlobalV::MY_RANK, GlobalV::global_readin_dir); remove("Restart_md.dat"); FIRE *fire = dynamic_cast(mdrun); @@ -167,10 +167,10 @@ TEST_F(FIRE_test, restart) EXPECT_EQ(mdrun->mdp.md_dt, 41.3414); } -TEST_F(FIRE_test, outputMD) +TEST_F(FIREtest, OutputMD) { std::ofstream ofs("running.log"); - mdrun->outputMD(ofs, true); + mdrun->outputMD(ofs, true, GlobalV::MY_RANK); ofs.close(); std::ifstream ifs("running.log"); diff --git a/source/module_md/test/langevin_test.cpp b/source/module_md/test/langevin_test.cpp index 6d08391445..b69b1434f4 100644 --- a/source/module_md/test/langevin_test.cpp +++ b/source/module_md/test/langevin_test.cpp @@ -47,7 +47,7 @@ class Langevin_test : public testing::Test p_esolver->Init(INPUT, ucell); mdrun = new Langevin(INPUT.mdp, ucell); - mdrun->setup(p_esolver); + mdrun->setup(p_esolver, GlobalV::MY_RANK, GlobalV::global_readin_dir); } void TearDown() @@ -72,7 +72,7 @@ TEST_F(Langevin_test, setup) TEST_F(Langevin_test, first_half) { - mdrun->first_half(); + mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); EXPECT_NEAR(mdrun->pos[0].x, -0.00042883345359910814, doublethreshold); EXPECT_NEAR(mdrun->pos[0].y, 0.00016393608896004904, doublethreshold); @@ -103,8 +103,8 @@ TEST_F(Langevin_test, first_half) TEST_F(Langevin_test, second_half) { - mdrun->first_half(); - mdrun->second_half(); + mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); + mdrun->second_half(GlobalV::MY_RANK);; EXPECT_NEAR(mdrun->pos[0].x, -0.00066954020090275205, doublethreshold); EXPECT_NEAR(mdrun->pos[0].y, 3.3862365219131354e-05, doublethreshold); @@ -137,7 +137,7 @@ TEST_F(Langevin_test, write_restart) { mdrun->step_ = 1; mdrun->step_rst_ = 2; - mdrun->write_restart(); + mdrun->write_restart(GlobalV::MY_RANK, GlobalV::global_out_dir); std::ifstream ifs("Restart_md.dat"); std::string output_str; @@ -148,7 +148,7 @@ TEST_F(Langevin_test, write_restart) TEST_F(Langevin_test, restart) { - mdrun->restart(); + mdrun->restart(GlobalV::MY_RANK, GlobalV::global_readin_dir); remove("Restart_md.dat"); EXPECT_EQ(mdrun->step_rst_, 3); @@ -157,7 +157,7 @@ TEST_F(Langevin_test, restart) TEST_F(Langevin_test, outputMD) { std::ofstream ofs("running.log"); - mdrun->outputMD(ofs, true); + mdrun->outputMD(ofs, true, GlobalV::MY_RANK); ofs.close(); std::ifstream ifs("running.log"); diff --git a/source/module_md/test/md_func_test.cpp b/source/module_md/test/md_func_test.cpp index 84ef39ffff..379c39f32b 100644 --- a/source/module_md/test/md_func_test.cpp +++ b/source/module_md/test/md_func_test.cpp @@ -206,7 +206,7 @@ TEST_F(MD_func_test, MDdump) ifs.close(); // append - MD_func::MDdump(1, ucell, INPUT, virial, force, vel); + MD_func::MDdump(1, 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/msst_test.cpp b/source/module_md/test/msst_test.cpp index 349fd17bbc..fdbd8c1019 100644 --- a/source/module_md/test/msst_test.cpp +++ b/source/module_md/test/msst_test.cpp @@ -47,7 +47,7 @@ class MSST_test : public testing::Test p_esolver->Init(INPUT, ucell); mdrun = new MSST(INPUT.mdp, ucell); - mdrun->setup(p_esolver); + mdrun->setup(p_esolver, GlobalV::MY_RANK, GlobalV::global_readin_dir); } void TearDown() @@ -84,7 +84,7 @@ TEST_F(MSST_test, setup) TEST_F(MSST_test, first_half) { - mdrun->first_half(); + mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); EXPECT_NEAR(ucell.lat0, 1.0, doublethreshold); EXPECT_NEAR(ucell.lat0_angstrom, 0.52917700000000001, doublethreshold); @@ -128,8 +128,8 @@ TEST_F(MSST_test, first_half) TEST_F(MSST_test, second_half) { - mdrun->first_half(); - mdrun->second_half(); + mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); + mdrun->second_half(GlobalV::MY_RANK);; EXPECT_NEAR(ucell.lat0, 1.0, doublethreshold); EXPECT_NEAR(ucell.lat0_angstrom, 0.52917700000000001, doublethreshold); @@ -175,7 +175,7 @@ TEST_F(MSST_test, write_restart) { mdrun->step_ = 1; mdrun->step_rst_ = 2; - mdrun->write_restart(); + mdrun->write_restart(GlobalV::MY_RANK, GlobalV::global_out_dir); std::ifstream ifs("Restart_md.dat"); std::string output_str; @@ -196,7 +196,7 @@ TEST_F(MSST_test, write_restart) TEST_F(MSST_test, restart) { - mdrun->restart(); + mdrun->restart(GlobalV::MY_RANK, GlobalV::global_readin_dir); remove("Restart_md.dat"); MSST *msst = dynamic_cast(mdrun); @@ -211,7 +211,7 @@ TEST_F(MSST_test, restart) TEST_F(MSST_test, outputMD) { std::ofstream ofs("running.log"); - mdrun->outputMD(ofs, true); + mdrun->outputMD(ofs, true, GlobalV::MY_RANK); ofs.close(); std::ifstream ifs("running.log"); diff --git a/source/module_md/test/nhchain_test.cpp b/source/module_md/test/nhchain_test.cpp index a9aedfa866..6299f06bc9 100644 --- a/source/module_md/test/nhchain_test.cpp +++ b/source/module_md/test/nhchain_test.cpp @@ -49,7 +49,7 @@ class NHC_test : public testing::Test INPUT.mdp.md_type = "npt"; INPUT.mdp.md_pmode = "tri"; mdrun = new Nose_Hoover(INPUT.mdp, ucell); - mdrun->setup(p_esolver); + mdrun->setup(p_esolver, GlobalV::MY_RANK, GlobalV::global_readin_dir); } void TearDown() @@ -74,7 +74,7 @@ TEST_F(NHC_test, setup) TEST_F(NHC_test, first_half) { - mdrun->first_half(); + mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); EXPECT_NEAR(mdrun->pos[0].x, -0.00023793471204889866, doublethreshold); EXPECT_NEAR(mdrun->pos[0].y, 0.00017779705725471447, doublethreshold); @@ -105,8 +105,8 @@ TEST_F(NHC_test, first_half) TEST_F(NHC_test, second_half) { - mdrun->first_half(); - mdrun->second_half(); + mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); + mdrun->second_half(GlobalV::MY_RANK);; EXPECT_NEAR(mdrun->pos[0].x, -0.00023793503786683287, doublethreshold); EXPECT_NEAR(mdrun->pos[0].y, 0.0001777972998948069, doublethreshold); @@ -137,11 +137,11 @@ TEST_F(NHC_test, second_half) TEST_F(NHC_test, write_restart) { - mdrun->first_half(); - mdrun->second_half(); + mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); + mdrun->second_half(GlobalV::MY_RANK);; mdrun->step_ = 1; mdrun->step_rst_ = 2; - mdrun->write_restart(); + mdrun->write_restart(GlobalV::MY_RANK, GlobalV::global_out_dir); std::ifstream ifs("Restart_md.dat"); std::string output_str; @@ -168,7 +168,7 @@ TEST_F(NHC_test, write_restart) TEST_F(NHC_test, restart) { - mdrun->restart(); + mdrun->restart(GlobalV::MY_RANK, GlobalV::global_readin_dir); remove("Restart_md.dat"); Nose_Hoover *nhc = dynamic_cast(mdrun); @@ -202,7 +202,7 @@ TEST_F(NHC_test, restart) TEST_F(NHC_test, outputMD) { std::ofstream ofs("running.log"); - mdrun->outputMD(ofs, true); + mdrun->outputMD(ofs, true, GlobalV::MY_RANK); ofs.close(); std::ifstream ifs("running.log"); diff --git a/source/module_md/test/verlet_test.cpp b/source/module_md/test/verlet_test.cpp index c873a730c7..5055b5fa23 100644 --- a/source/module_md/test/verlet_test.cpp +++ b/source/module_md/test/verlet_test.cpp @@ -47,7 +47,7 @@ class Verlet_test : public testing::Test p_esolver->Init(INPUT, ucell); mdrun = new Verlet(INPUT.mdp, ucell); - mdrun->setup(p_esolver); + mdrun->setup(p_esolver, GlobalV::MY_RANK, GlobalV::global_readin_dir); } void TearDown() @@ -72,7 +72,7 @@ TEST_F(Verlet_test, setup) TEST_F(Verlet_test, first_half) { - mdrun->first_half(); + mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); EXPECT_NEAR(mdrun->pos[0].x, -0.00054545529007222658, doublethreshold); EXPECT_NEAR(mdrun->pos[0].y, 0.00029590658162135359, doublethreshold); @@ -103,9 +103,9 @@ TEST_F(Verlet_test, first_half) TEST_F(Verlet_test, NVE) { - mdrun->first_half(); + mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); mdrun->mdp.md_type = "nve"; - mdrun->second_half(); + mdrun->second_half(GlobalV::MY_RANK);; EXPECT_NEAR(mdrun->pos[0].x, -0.00054545529007222658, doublethreshold); EXPECT_NEAR(mdrun->pos[0].y, 0.00029590658162135359, doublethreshold); @@ -136,10 +136,10 @@ TEST_F(Verlet_test, NVE) TEST_F(Verlet_test, Anderson) { - mdrun->first_half(); + mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); mdrun->mdp.md_type = "nvt"; mdrun->mdp.md_thermostat = "anderson"; - mdrun->second_half(); + mdrun->second_half(GlobalV::MY_RANK);; EXPECT_NEAR(mdrun->pos[0].x, -0.00054545529007222658, doublethreshold); EXPECT_NEAR(mdrun->pos[0].y, 0.00029590658162135359, doublethreshold); @@ -170,10 +170,10 @@ TEST_F(Verlet_test, Anderson) TEST_F(Verlet_test, Berendsen) { - mdrun->first_half(); + mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); mdrun->mdp.md_type = "nvt"; mdrun->mdp.md_thermostat = "berendsen"; - mdrun->second_half(); + mdrun->second_half(GlobalV::MY_RANK);; EXPECT_NEAR(mdrun->pos[0].x, -0.00054545529007222658, doublethreshold); EXPECT_NEAR(mdrun->pos[0].y, 0.00029590658162135359, doublethreshold); @@ -204,10 +204,10 @@ TEST_F(Verlet_test, Berendsen) TEST_F(Verlet_test, rescaling) { - mdrun->first_half(); + mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); mdrun->mdp.md_type = "nvt"; mdrun->mdp.md_thermostat = "rescaling"; - mdrun->second_half(); + mdrun->second_half(GlobalV::MY_RANK);; EXPECT_NEAR(mdrun->pos[0].x, -0.00054545529007222658, doublethreshold); EXPECT_NEAR(mdrun->pos[0].y, 0.00029590658162135359, doublethreshold); @@ -238,10 +238,10 @@ TEST_F(Verlet_test, rescaling) TEST_F(Verlet_test, rescale_v) { - mdrun->first_half(); + mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); mdrun->mdp.md_type = "nvt"; mdrun->mdp.md_thermostat = "rescale_v"; - mdrun->second_half(); + mdrun->second_half(GlobalV::MY_RANK);; EXPECT_NEAR(mdrun->pos[0].x, -0.00054545529007222658, doublethreshold); EXPECT_NEAR(mdrun->pos[0].y, 0.00029590658162135359, doublethreshold); @@ -274,7 +274,7 @@ TEST_F(Verlet_test, write_restart) { mdrun->step_ = 1; mdrun->step_rst_ = 2; - mdrun->write_restart(); + mdrun->write_restart(GlobalV::MY_RANK, GlobalV::global_out_dir); std::ifstream ifs("Restart_md.dat"); std::string output_str; @@ -285,7 +285,7 @@ TEST_F(Verlet_test, write_restart) TEST_F(Verlet_test, restart) { - mdrun->restart(); + mdrun->restart(GlobalV::MY_RANK, GlobalV::global_readin_dir); remove("Restart_md.dat"); EXPECT_EQ(mdrun->step_rst_, 3); @@ -294,7 +294,7 @@ TEST_F(Verlet_test, restart) TEST_F(Verlet_test, outputMD) { std::ofstream ofs("running.log"); - mdrun->outputMD(ofs, true); + mdrun->outputMD(ofs, true, GlobalV::MY_RANK); ofs.close(); std::ifstream ifs("running.log"); diff --git a/source/module_md/verlet.cpp b/source/module_md/verlet.cpp index 9c60db7172..68987aa1b0 100644 --- a/source/module_md/verlet.cpp +++ b/source/module_md/verlet.cpp @@ -16,12 +16,12 @@ void Verlet::setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const ModuleBase::TITLE("Verlet", "setup"); ModuleBase::timer::tick("Verlet", "setup"); - MDrun::setup(p_ensolve); + MDrun::setup(p_esolver, my_rank, global_readin_dir); ModuleBase::timer::tick("Verlet", "setup"); } -void Verlet::first_half(const int &my_rank) +void Verlet::first_half(const int &my_rank, std::ofstream &ofs) { ModuleBase::TITLE("Verlet", "first_half"); ModuleBase::timer::tick("Verlet", "first_half"); diff --git a/source/module_md/verlet.h b/source/module_md/verlet.h index 5a4f1f6f61..4396e2bcc3 100644 --- a/source/module_md/verlet.h +++ b/source/module_md/verlet.h @@ -10,7 +10,7 @@ class Verlet : public MDrun ~Verlet(); void setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const std::string &global_readin_dir); - void first_half(const int &my_rank); + void first_half(const int &my_rank, std::ofstream &ofs); void second_half(const int &my_rank); void apply_thermostat(const int &my_rank); void thermalize(const int &nraise, const double ¤t_temp, const double &target_temp); From b816bc2b1e06d3f0aaa398e30eced79d41846c9e Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Sun, 7 May 2023 16:12:51 +0800 Subject: [PATCH 4/5] Refactor: rename MDrun to MD_base --- source/driver_run.cpp | 3 +- source/module_md/fire.cpp | 12 +++---- source/module_md/fire.h | 7 ++-- source/module_md/langevin.cpp | 16 ++++----- source/module_md/langevin.h | 3 +- source/module_md/md_base.cpp | 20 +++++------ source/module_md/md_base.h | 46 +++++++++++++------------ source/module_md/msst.cpp | 8 ++--- source/module_md/msst.h | 3 +- source/module_md/nhchain.cpp | 12 +++---- source/module_md/nhchain.h | 3 +- source/module_md/run_md.cpp | 36 ++++++++----------- source/module_md/run_md.h | 14 +++----- source/module_md/test/fire_test.cpp | 4 ++- source/module_md/test/langevin_test.cpp | 4 ++- source/module_md/test/msst_test.cpp | 4 ++- source/module_md/test/nhchain_test.cpp | 4 ++- source/module_md/test/verlet_test.cpp | 4 ++- source/module_md/verlet.cpp | 16 ++++----- source/module_md/verlet.h | 11 +++--- 20 files changed, 117 insertions(+), 113 deletions(-) diff --git a/source/driver_run.cpp b/source/driver_run.cpp index 224e725a96..2e8d4a9295 100644 --- a/source/driver_run.cpp +++ b/source/driver_run.cpp @@ -48,8 +48,7 @@ void Driver::driver_run() //---------------------------MD/Relax------------------------- if(GlobalV::CALCULATION == "md") { - Run_MD run_md; - run_md.md_line(GlobalC::ucell, p_esolver, INPUT.mdp); + Run_MD::md_line(GlobalC::ucell, p_esolver, INPUT.mdp); } else // scf; cell relaxation; nscf; etc { diff --git a/source/module_md/fire.cpp b/source/module_md/fire.cpp index 62b41fb50a..5fc4c691db 100644 --- a/source/module_md/fire.cpp +++ b/source/module_md/fire.cpp @@ -6,7 +6,7 @@ #endif #include "module_base/timer.h" -FIRE::FIRE(MD_parameters &MD_para_in, UnitCell &unit_in) : MDrun(MD_para_in, unit_in) +FIRE::FIRE(MD_parameters &MD_para_in, UnitCell &unit_in) : MD_base(MD_para_in, unit_in) { dt_max = -1.0; alpha_start = 0.10; @@ -28,7 +28,7 @@ void FIRE::setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const st ModuleBase::TITLE("FIRE", "setup"); ModuleBase::timer::tick("FIRE", "setup"); - MDrun::setup(p_esolver, my_rank, global_readin_dir); + MD_base::setup(p_esolver, my_rank, global_readin_dir); check_force(); @@ -40,11 +40,11 @@ void FIRE::first_half(const int &my_rank, std::ofstream &ofs) ModuleBase::TITLE("FIRE", "first_half"); ModuleBase::timer::tick("FIRE", "first_half"); - MDrun::update_vel(force, my_rank); + MD_base::update_vel(force, my_rank); check_FIRE(); - MDrun::update_pos(my_rank); + MD_base::update_pos(my_rank); ModuleBase::timer::tick("FIRE", "first_half"); } @@ -54,7 +54,7 @@ void FIRE::second_half(const int &my_rank) ModuleBase::TITLE("FIRE", "second_half"); ModuleBase::timer::tick("FIRE", "second_half"); - MDrun::update_vel(force, my_rank); + MD_base::update_vel(force, my_rank); check_force(); @@ -63,7 +63,7 @@ void FIRE::second_half(const int &my_rank) void FIRE::outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank) { - MDrun::outputMD(ofs, cal_stress, my_rank); + MD_base::outputMD(ofs, cal_stress, my_rank); ofs << " LARGEST GRAD (eV/A) : " << max * ModuleBase::Hartree_to_eV * ModuleBase::ANGSTROM_AU << std::endl; std::cout << " LARGEST GRAD (eV/A) : " << max * ModuleBase::Hartree_to_eV * ModuleBase::ANGSTROM_AU << std::endl; diff --git a/source/module_md/fire.h b/source/module_md/fire.h index 7ab02c33a4..954f3cff9b 100644 --- a/source/module_md/fire.h +++ b/source/module_md/fire.h @@ -3,18 +3,19 @@ #include "md_base.h" -class FIRE : public MDrun +class FIRE : public MD_base { public: FIRE(MD_parameters &MD_para_in, UnitCell &unit_in); ~FIRE(); + private: void setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const std::string &global_readin_dir); void first_half(const int &my_rank, std::ofstream &ofs); void second_half(const int &my_rank); void outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank); - void write_restart(const int &my_rank, const std::string &global_out_dir); - void restart(const int &my_rank, const std::string &global_readin_dir); + void restart(const int& my_rank, const std::string& global_readin_dir); + void write_restart(const int& my_rank, const std::string& global_out_dir); void check_force(); void check_FIRE(); diff --git a/source/module_md/langevin.cpp b/source/module_md/langevin.cpp index 55fa052be9..17c9a1f871 100644 --- a/source/module_md/langevin.cpp +++ b/source/module_md/langevin.cpp @@ -4,7 +4,7 @@ #include "module_base/parallel_common.h" #include "module_base/timer.h" -Langevin::Langevin(MD_parameters &MD_para_in, UnitCell &unit_in) : MDrun(MD_para_in, unit_in) +Langevin::Langevin(MD_parameters &MD_para_in, UnitCell &unit_in) : MD_base(MD_para_in, unit_in) { // convert to a.u. unit mdp.md_damp /= ModuleBase::AU_to_FS; @@ -22,7 +22,7 @@ void Langevin::setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, cons ModuleBase::TITLE("Langevin", "setup"); ModuleBase::timer::tick("Langevin", "setup"); - MDrun::setup(p_esolver, my_rank, global_readin_dir); + MD_base::setup(p_esolver, my_rank, global_readin_dir); post_force(); @@ -34,8 +34,8 @@ void Langevin::first_half(const int &my_rank, std::ofstream &ofs) ModuleBase::TITLE("Langevin", "first_half"); ModuleBase::timer::tick("Langevin", "first_half"); - MDrun::update_vel(total_force, my_rank); - MDrun::update_pos(my_rank); + MD_base::update_vel(total_force, my_rank); + MD_base::update_pos(my_rank); ModuleBase::timer::tick("Langevin", "first_half"); } @@ -46,24 +46,24 @@ void Langevin::second_half(const int &my_rank) ModuleBase::timer::tick("Langevin", "second_half"); post_force(); - MDrun::update_vel(total_force, my_rank); + MD_base::update_vel(total_force, my_rank); ModuleBase::timer::tick("Langevin", "second_half"); } void Langevin::outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank) { - MDrun::outputMD(ofs, cal_stress, my_rank); + MD_base::outputMD(ofs, cal_stress, my_rank); } void Langevin::write_restart(const int &my_rank, const std::string &global_out_dir) { - MDrun::write_restart(my_rank, global_out_dir); + MD_base::write_restart(my_rank, global_out_dir); } void Langevin::restart(const int &my_rank, const std::string &global_readin_dir) { - MDrun::restart(my_rank, global_readin_dir); + MD_base::restart(my_rank, global_readin_dir); } void Langevin::post_force() diff --git a/source/module_md/langevin.h b/source/module_md/langevin.h index 44ea9a6fba..967cf49382 100644 --- a/source/module_md/langevin.h +++ b/source/module_md/langevin.h @@ -3,12 +3,13 @@ #include "md_base.h" -class Langevin : public MDrun +class Langevin : public MD_base { public: Langevin(MD_parameters &MD_para_in, UnitCell &unit_in); ~Langevin(); + private: void setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const std::string &global_readin_dir); void first_half(const int &my_rank, std::ofstream &ofs); void second_half(const int &my_rank); diff --git a/source/module_md/md_base.cpp b/source/module_md/md_base.cpp index 0449b871f8..2b25a68322 100644 --- a/source/module_md/md_base.cpp +++ b/source/module_md/md_base.cpp @@ -6,7 +6,7 @@ #endif #include "module_io/print_info.h" -MDrun::MDrun(MD_parameters &MD_para_in, UnitCell &unit_in) : mdp(MD_para_in), ucell(unit_in) +MD_base::MD_base(MD_parameters& MD_para_in, UnitCell& unit_in) : mdp(MD_para_in), ucell(unit_in) { if (mdp.md_seed >= 0) { @@ -35,7 +35,7 @@ MDrun::MDrun(MD_parameters &MD_para_in, UnitCell &unit_in) : mdp(MD_para_in), uc MD_func::InitVel(ucell, mdp.md_tfirst, allmass, frozen_freedom_, ionmbl, vel); } -MDrun::~MDrun() +MD_base::~MD_base() { delete[] allmass; delete[] pos; @@ -44,7 +44,7 @@ MDrun::~MDrun() delete[] force; } -void MDrun::setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const std::string &global_readin_dir) +void MD_base::setup(ModuleESolver::ESolver* p_esolver, const int& my_rank, const std::string& global_readin_dir) { if (mdp.md_restart) { @@ -59,18 +59,18 @@ void MDrun::setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const s ucell.ionic_position_updated = true; } -void MDrun::first_half(const int &my_rank, std::ofstream &ofs) +void MD_base::first_half(const int& my_rank, std::ofstream& ofs) { update_vel(force, my_rank); update_pos(my_rank); } -void MDrun::second_half(const int &my_rank) +void MD_base::second_half(const int& my_rank) { update_vel(force, my_rank); } -void MDrun::update_pos(const int &my_rank) +void MD_base::update_pos(const int& my_rank) { if (my_rank == 0) { @@ -98,7 +98,7 @@ void MDrun::update_pos(const int &my_rank) ucell.update_pos_taud(pos); } -void MDrun::update_vel(const ModuleBase::Vector3 *force, const int &my_rank) +void MD_base::update_vel(const ModuleBase::Vector3* force, const int& my_rank) { if (my_rank == 0) { @@ -119,7 +119,7 @@ void MDrun::update_vel(const ModuleBase::Vector3 *force, const int &my_r #endif } -void MDrun::outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank) +void MD_base::outputMD(std::ofstream& ofs, const bool& cal_stress, const int& my_rank) { if (my_rank) return; @@ -183,7 +183,7 @@ void MDrun::outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_r ofs << std::endl; } -void MDrun::write_restart(const int &my_rank, const std::string &global_out_dir) +void MD_base::write_restart(const int& my_rank, const std::string& global_out_dir) { if (!my_rank) { @@ -199,7 +199,7 @@ void MDrun::write_restart(const int &my_rank, const std::string &global_out_dir) #endif } -void MDrun::restart(const int &my_rank, const std::string &global_readin_dir) +void MD_base::restart(const int& my_rank, const std::string& global_readin_dir) { step_rst_ = MD_func::current_step(my_rank, global_readin_dir); } \ No newline at end of file diff --git a/source/module_md/md_base.h b/source/module_md/md_base.h index 39c28d1283..67079e7525 100644 --- a/source/module_md/md_base.h +++ b/source/module_md/md_base.h @@ -4,11 +4,11 @@ #include "md_para.h" #include "module_esolver/esolver.h" -class MDrun +class MD_base { public: - MDrun(MD_parameters &MD_para_in, UnitCell &unit_in); - virtual ~MDrun(); + MD_base(MD_parameters& MD_para_in, UnitCell& unit_in); + virtual ~MD_base(); /** * @brief init before running md, calculate energy, force, and stress of the initial configuration. @@ -31,19 +31,6 @@ class MDrun */ virtual void second_half(const int &my_rank); - /** - * @brief perform one step update of pos due to atomic velocity - * @param my_rank MPI rank of the processer - */ - virtual void update_pos(const int &my_rank); - - /** - * @brief perform half-step update of vel due to atomic force - * @param force atomic forces - * @param my_rank MPI rank of the processer - */ - virtual void update_vel(const ModuleBase::Vector3 *force, const int &my_rank); - /** * @brief output MD information such as energy, temperature, and pressure * @param ofs determine the output files @@ -59,24 +46,34 @@ class MDrun */ virtual void write_restart(const int &my_rank, const std::string &global_out_dir); + protected: /** * @brief restart MD when md_restart is true * @param my_rank MPI rank of the processer * @param global_readin_dir directory of files for reading */ - virtual void restart(const int &my_rank, const std::string &global_readin_dir); + virtual void restart(const int& my_rank, const std::string& global_readin_dir); + + /** + * @brief perform one step update of pos due to atomic velocity + * @param my_rank MPI rank of the processer + */ + virtual void update_pos(const int& my_rank); - MD_parameters &mdp; - UnitCell &ucell; - bool stop; // MD stop or not + /** + * @brief perform half-step update of vel due to atomic force + * @param force atomic forces + * @param my_rank MPI rank of the processer + */ + virtual void update_vel(const ModuleBase::Vector3* force, const int& my_rank); // All parameters are in a.u. unit. + public: + bool stop; // MD stop or not double t_current; // current temperature int step_; // the MD step finished in current calculation int step_rst_; // the MD step finished in previous calculations - double energy_; // total energy of the system int frozen_freedom_; // the fixed freedom of the system - double *allmass; // atom mass ModuleBase::Vector3 *pos; // atom displacements liuyu modify 2023-03-22 ModuleBase::Vector3 *vel; // atom velocity @@ -86,6 +83,11 @@ class MDrun ModuleBase::matrix stress; // stress for this lattice double potential; // potential energy double kinetic; // kinetic energy + + protected: + MD_parameters& mdp; + UnitCell& ucell; + double energy_; // total energy of the system }; #endif \ No newline at end of file diff --git a/source/module_md/msst.cpp b/source/module_md/msst.cpp index 5f5bdf00ed..e9b3bae14f 100644 --- a/source/module_md/msst.cpp +++ b/source/module_md/msst.cpp @@ -6,7 +6,7 @@ #endif #include "module_base/timer.h" -MSST::MSST(MD_parameters &MD_para_in, UnitCell &unit_in) : MDrun(MD_para_in, unit_in) +MSST::MSST(MD_parameters &MD_para_in, UnitCell &unit_in) : MD_base(MD_para_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; @@ -36,7 +36,7 @@ void MSST::setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const st ModuleBase::TITLE("MSST", "setup"); ModuleBase::timer::tick("MSST", "setup"); - MDrun::setup(p_esolver, my_rank, global_readin_dir); + MD_base::setup(p_esolver, my_rank, global_readin_dir); ucell.cell_parameter_updated = true; int sd = mdp.msst_direction; @@ -111,7 +111,7 @@ void MSST::first_half(const int &my_rank, std::ofstream &ofs) rescale(ofs, vol); // propagate atom positions 1 time step - MDrun::update_pos(my_rank); + MD_base::update_pos(my_rank); // propagate volume 1/2 step vol = ucell.omega + omega[sd] * dthalf; @@ -149,7 +149,7 @@ void MSST::second_half(const int &my_rank) void MSST::outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank) { - MDrun::outputMD(ofs, cal_stress, my_rank); + MD_base::outputMD(ofs, cal_stress, my_rank); } void MSST::write_restart(const int &my_rank, const std::string &global_out_dir) diff --git a/source/module_md/msst.h b/source/module_md/msst.h index 1d4a342b7c..dfcf0f2362 100644 --- a/source/module_md/msst.h +++ b/source/module_md/msst.h @@ -3,12 +3,13 @@ #include "md_base.h" -class MSST : public MDrun +class MSST : public MD_base { public: MSST(MD_parameters &MD_para_in, UnitCell &unit_in); ~MSST(); + private: void setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const std::string &global_readin_dir); void first_half(const int &my_rank, std::ofstream &ofs); void second_half(const int &my_rank); diff --git a/source/module_md/nhchain.cpp b/source/module_md/nhchain.cpp index 98ca2adf89..bcc0162957 100644 --- a/source/module_md/nhchain.cpp +++ b/source/module_md/nhchain.cpp @@ -6,7 +6,7 @@ #endif #include "module_base/timer.h" -Nose_Hoover::Nose_Hoover(MD_parameters &MD_para_in, UnitCell &unit_in) : MDrun(MD_para_in, unit_in) +Nose_Hoover::Nose_Hoover(MD_parameters &MD_para_in, UnitCell &unit_in) : MD_base(MD_para_in, unit_in) { const double unit_transform = ModuleBase::HARTREE_SI / pow(ModuleBase::BOHR_RADIUS_SI, 3) * 1.0e-8; mdp.md_pfirst /= unit_transform; @@ -152,7 +152,7 @@ void Nose_Hoover::setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, c ModuleBase::TITLE("Nose_Hoover", "setup"); ModuleBase::timer::tick("Nose_Hoover", "setup"); - MDrun::setup(p_esolver, my_rank, global_readin_dir); + MD_base::setup(p_esolver, my_rank, global_readin_dir); if (mdp.md_type == "npt") { ucell.cell_parameter_updated = true; @@ -241,7 +241,7 @@ void Nose_Hoover::first_half(const int &my_rank, std::ofstream &ofs) } // perform half-step update of vel due to atomic force - MDrun::update_vel(force, my_rank); + MD_base::update_vel(force, my_rank); if (npt_flag) { @@ -250,7 +250,7 @@ void Nose_Hoover::first_half(const int &my_rank, std::ofstream &ofs) } // perform one step update of pos due to atomic velocity - MDrun::update_pos(my_rank); + MD_base::update_pos(my_rank); if (npt_flag) { @@ -267,7 +267,7 @@ void Nose_Hoover::second_half(const int &my_rank) ModuleBase::timer::tick("Nose_Hoover", "second_half"); // perform half-step update of vel due to atomic force - MDrun::update_vel(force, my_rank); + MD_base::update_vel(force, my_rank); if (npt_flag) { @@ -304,7 +304,7 @@ void Nose_Hoover::second_half(const int &my_rank) void Nose_Hoover::outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank) { - MDrun::outputMD(ofs, cal_stress, my_rank); + MD_base::outputMD(ofs, cal_stress, my_rank); } void Nose_Hoover::write_restart(const int &my_rank, const std::string &global_out_dir) diff --git a/source/module_md/nhchain.h b/source/module_md/nhchain.h index 7662ad8ba5..26e18c4706 100644 --- a/source/module_md/nhchain.h +++ b/source/module_md/nhchain.h @@ -3,12 +3,13 @@ #include "md_base.h" -class Nose_Hoover : public MDrun +class Nose_Hoover : public MD_base { public: Nose_Hoover(MD_parameters &MD_para_in, UnitCell &unit_in); ~Nose_Hoover(); + private: void setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const std::string &global_readin_dir); void first_half(const int &my_rank, std::ofstream &ofs); void second_half(const int &my_rank); diff --git a/source/module_md/run_md.cpp b/source/module_md/run_md.cpp index 58f383844f..976aa869da 100644 --- a/source/module_md/run_md.cpp +++ b/source/module_md/run_md.cpp @@ -9,21 +9,16 @@ #include "nhchain.h" #include "verlet.h" -Run_MD::Run_MD() +namespace Run_MD { -} - -Run_MD::~Run_MD() -{ -} -void Run_MD::md_line(UnitCell &unit_in, ModuleESolver::ESolver *p_esolver, MD_parameters &md_para) +void md_line(UnitCell& unit_in, ModuleESolver::ESolver* p_esolver, MD_parameters& md_para) { ModuleBase::TITLE("Run_MD", "md_line"); ModuleBase::timer::tick("Run_MD", "md_line"); // determine the md_type - MDrun *mdrun; + MD_base* mdrun; if (md_para.md_type == "fire") { mdrun = new FIRE(md_para, unit_in); @@ -62,36 +57,31 @@ void Run_MD::md_line(UnitCell &unit_in, ModuleESolver::ESolver *p_esolver, MD_pa mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); // update force and virial due to the update of atom positions - MD_func::force_virial(p_esolver, mdrun->step_, mdrun->ucell, mdrun->potential, mdrun->force, mdrun->virial); + MD_func::force_virial(p_esolver, mdrun->step_, unit_in, mdrun->potential, mdrun->force, mdrun->virial); mdrun->second_half(GlobalV::MY_RANK); - MD_func::compute_stress(mdrun->ucell, mdrun->vel, mdrun->allmass, mdrun->virial, mdrun->stress); + MD_func::compute_stress(unit_in, mdrun->vel, mdrun->allmass, mdrun->virial, mdrun->stress); mdrun->t_current = MD_func::current_temp(mdrun->kinetic, - mdrun->ucell.nat, + unit_in.nat, mdrun->frozen_freedom_, mdrun->allmass, mdrun->vel); } - if ((mdrun->step_ + mdrun->step_rst_) % mdrun->mdp.md_dumpfreq == 0) + if ((mdrun->step_ + mdrun->step_rst_) % md_para.md_dumpfreq == 0) { mdrun->outputMD(GlobalV::ofs_running, GlobalV::CAL_STRESS, GlobalV::MY_RANK); - MD_func::MDdump(mdrun->step_ + mdrun->step_rst_, - mdrun->ucell, - md_para, - mdrun->virial, - mdrun->force, - mdrun->vel); + MD_func::MDdump(mdrun->step_ + mdrun->step_rst_, unit_in, md_para, mdrun->virial, mdrun->force, mdrun->vel); } - if ((mdrun->step_ + mdrun->step_rst_) % mdrun->mdp.md_restartfreq == 0) + if ((mdrun->step_ + mdrun->step_rst_) % md_para.md_restartfreq == 0) { - mdrun->ucell.update_vel(mdrun->vel); + unit_in.update_vel(mdrun->vel); std::stringstream file; file << GlobalV::global_stru_dir << "STRU_MD_" << mdrun->step_ + mdrun->step_rst_; - mdrun->ucell.print_stru_file(file.str(), 1, 1); + unit_in.print_stru_file(file.str(), 1, 1); mdrun->write_restart(GlobalV::MY_RANK, GlobalV::global_out_dir); } @@ -101,4 +91,6 @@ void Run_MD::md_line(UnitCell &unit_in, ModuleESolver::ESolver *p_esolver, MD_pa delete mdrun; ModuleBase::timer::tick("Run_MD", "md_line"); return; -} \ No newline at end of file +} + +} // namespace Run_MD \ No newline at end of file diff --git a/source/module_md/run_md.h b/source/module_md/run_md.h index 248041ddc6..3cbf676330 100644 --- a/source/module_md/run_md.h +++ b/source/module_md/run_md.h @@ -1,16 +1,12 @@ -#ifndef RUN_MD_CLASSIC_H -#define RUN_MD_CLASSIC_H +#ifndef RUN_MD_H +#define RUN_MD_H #include "md_para.h" #include "module_esolver/esolver.h" -class Run_MD +namespace Run_MD { - public: - Run_MD(); - ~Run_MD(); - - void md_line(UnitCell &unit_in, ModuleESolver::ESolver *p_esolver, MD_parameters &md_para); -}; +void md_line(UnitCell& unit_in, ModuleESolver::ESolver* p_esolver, MD_parameters& md_para); +} // namespace Run_MD #endif \ No newline at end of file diff --git a/source/module_md/test/fire_test.cpp b/source/module_md/test/fire_test.cpp index c3808affb8..9f8807767e 100644 --- a/source/module_md/test/fire_test.cpp +++ b/source/module_md/test/fire_test.cpp @@ -1,3 +1,5 @@ +#define private public +#define protected public #include "module_md/fire.h" #include "gmock/gmock.h" @@ -35,7 +37,7 @@ class FIREtest : public testing::Test { protected: - MDrun *mdrun; + MD_base *mdrun; UnitCell ucell; void SetUp() diff --git a/source/module_md/test/langevin_test.cpp b/source/module_md/test/langevin_test.cpp index b69b1434f4..bc2d408d63 100644 --- a/source/module_md/test/langevin_test.cpp +++ b/source/module_md/test/langevin_test.cpp @@ -1,3 +1,5 @@ +#define private public +#define protected public #include "module_md/langevin.h" #include "gmock/gmock.h" @@ -35,7 +37,7 @@ class Langevin_test : public testing::Test { protected: - MDrun *mdrun; + MD_base *mdrun; UnitCell ucell; void SetUp() diff --git a/source/module_md/test/msst_test.cpp b/source/module_md/test/msst_test.cpp index fdbd8c1019..47631d2a96 100644 --- a/source/module_md/test/msst_test.cpp +++ b/source/module_md/test/msst_test.cpp @@ -1,3 +1,5 @@ +#define private public +#define protected public #include "module_md/msst.h" #include "gmock/gmock.h" @@ -35,7 +37,7 @@ class MSST_test : public testing::Test { protected: - MDrun *mdrun; + MD_base *mdrun; UnitCell ucell; void SetUp() diff --git a/source/module_md/test/nhchain_test.cpp b/source/module_md/test/nhchain_test.cpp index 6299f06bc9..23acd65cbc 100644 --- a/source/module_md/test/nhchain_test.cpp +++ b/source/module_md/test/nhchain_test.cpp @@ -1,3 +1,5 @@ +#define private public +#define protected public #include "module_md/nhchain.h" #include "gmock/gmock.h" @@ -35,7 +37,7 @@ class NHC_test : public testing::Test { protected: - MDrun *mdrun; + MD_base *mdrun; UnitCell ucell; void SetUp() diff --git a/source/module_md/test/verlet_test.cpp b/source/module_md/test/verlet_test.cpp index 5055b5fa23..9c6b8f0324 100644 --- a/source/module_md/test/verlet_test.cpp +++ b/source/module_md/test/verlet_test.cpp @@ -1,3 +1,5 @@ +#define private public +#define protected public #include "module_md/verlet.h" #include "gmock/gmock.h" @@ -35,7 +37,7 @@ class Verlet_test : public testing::Test { protected: - MDrun *mdrun; + MD_base *mdrun; UnitCell ucell; void SetUp() diff --git a/source/module_md/verlet.cpp b/source/module_md/verlet.cpp index 68987aa1b0..3f08b69c09 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_parameters &MD_para_in, UnitCell &unit_in) : MDrun(MD_para_in, unit_in) +Verlet::Verlet(MD_parameters &MD_para_in, UnitCell &unit_in) : MD_base(MD_para_in, unit_in) { } @@ -16,7 +16,7 @@ void Verlet::setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const ModuleBase::TITLE("Verlet", "setup"); ModuleBase::timer::tick("Verlet", "setup"); - MDrun::setup(p_esolver, my_rank, global_readin_dir); + MD_base::setup(p_esolver, my_rank, global_readin_dir); ModuleBase::timer::tick("Verlet", "setup"); } @@ -26,8 +26,8 @@ void Verlet::first_half(const int &my_rank, std::ofstream &ofs) ModuleBase::TITLE("Verlet", "first_half"); ModuleBase::timer::tick("Verlet", "first_half"); - MDrun::update_vel(force, my_rank); - MDrun::update_pos(my_rank); + MD_base::update_vel(force, my_rank); + MD_base::update_pos(my_rank); ModuleBase::timer::tick("Verlet", "first_half"); } @@ -37,7 +37,7 @@ void Verlet::second_half(const int &my_rank) ModuleBase::TITLE("Verlet", "second_half"); ModuleBase::timer::tick("Verlet", "second_half"); - MDrun::update_vel(force, my_rank); + MD_base::update_vel(force, my_rank); apply_thermostat(my_rank); ModuleBase::timer::tick("Verlet", "second_half"); @@ -122,15 +122,15 @@ void Verlet::thermalize(const int &nraise, const double ¤t_temp, const dou void Verlet::outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank) { - MDrun::outputMD(ofs, cal_stress, my_rank); + MD_base::outputMD(ofs, cal_stress, my_rank); } void Verlet::write_restart(const int &my_rank, const std::string &global_out_dir) { - MDrun::write_restart(my_rank, global_out_dir); + MD_base::write_restart(my_rank, global_out_dir); } void Verlet::restart(const int &my_rank, const std::string &global_readin_dir) { - MDrun::restart(my_rank, global_readin_dir); + MD_base::restart(my_rank, global_readin_dir); } \ No newline at end of file diff --git a/source/module_md/verlet.h b/source/module_md/verlet.h index 4396e2bcc3..17b3c5400d 100644 --- a/source/module_md/verlet.h +++ b/source/module_md/verlet.h @@ -3,20 +3,21 @@ #include "md_base.h" -class Verlet : public MDrun +class Verlet : public MD_base { public: Verlet(MD_parameters &MD_para_in, UnitCell &unit_in); ~Verlet(); + private: void setup(ModuleESolver::ESolver *p_esolver, const int &my_rank, const std::string &global_readin_dir); void first_half(const int &my_rank, std::ofstream &ofs); void second_half(const int &my_rank); - void apply_thermostat(const int &my_rank); + void restart(const int& my_rank, const std::string& global_readin_dir); + void outputMD(std::ofstream& ofs, const bool& cal_stress, const int& my_rank); + void apply_thermostat(const int& my_rank); void thermalize(const int &nraise, const double ¤t_temp, const double &target_temp); - void outputMD(std::ofstream &ofs, const bool &cal_stress, const int &my_rank); - void write_restart(const int &my_rank, const std::string &global_out_dir); - void restart(const int &my_rank, const std::string &global_readin_dir); + void write_restart(const int& my_rank, const std::string& global_out_dir); }; #endif \ No newline at end of file From d8978622bbe19748b0a52d59acfe515f568a77f3 Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Sun, 7 May 2023 18:04:32 +0800 Subject: [PATCH 5/5] Test: update head files in module_md tests --- source/module_md/test/fire_test.cpp | 14 ++++++------- source/module_md/test/langevin_test.cpp | 15 +++++++------- source/module_md/test/msst_test.cpp | 17 ++++++++-------- source/module_md/test/nhchain_test.cpp | 20 +++++++++--------- source/module_md/test/verlet_test.cpp | 27 +++++++++++++++---------- 5 files changed, 51 insertions(+), 42 deletions(-) diff --git a/source/module_md/test/fire_test.cpp b/source/module_md/test/fire_test.cpp index 9f8807767e..9b3c38bd11 100644 --- a/source/module_md/test/fire_test.cpp +++ b/source/module_md/test/fire_test.cpp @@ -1,12 +1,12 @@ -#define private public -#define protected public -#include "module_md/fire.h" - #include "gmock/gmock.h" #include "gtest/gtest.h" #include "module_esolver/esolver_lj.h" #include "setcell.h" +#define private public +#define protected public +#include "module_md/fire.h" + #define doublethreshold 1e-12 /************************************************ @@ -37,7 +37,7 @@ class FIREtest : public testing::Test { protected: - MD_base *mdrun; + MD_base* mdrun; UnitCell ucell; void SetUp() @@ -45,7 +45,7 @@ class FIREtest : public testing::Test Setcell::setupcell(ucell); Setcell::parameters(); - ModuleESolver::ESolver *p_esolver = new ModuleESolver::ESolver_LJ(); + ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); p_esolver->Init(INPUT, ucell); mdrun = new FIRE(INPUT.mdp, ucell); @@ -161,7 +161,7 @@ TEST_F(FIREtest, Restart) mdrun->restart(GlobalV::MY_RANK, GlobalV::global_readin_dir); remove("Restart_md.dat"); - FIRE *fire = dynamic_cast(mdrun); + FIRE* fire = dynamic_cast(mdrun); EXPECT_EQ(mdrun->step_rst_, 3); EXPECT_EQ(fire->alpha, 0.1); EXPECT_EQ(fire->negative_count, 0); diff --git a/source/module_md/test/langevin_test.cpp b/source/module_md/test/langevin_test.cpp index bc2d408d63..79c29caee0 100644 --- a/source/module_md/test/langevin_test.cpp +++ b/source/module_md/test/langevin_test.cpp @@ -1,12 +1,12 @@ -#define private public -#define protected public -#include "module_md/langevin.h" - #include "gmock/gmock.h" #include "gtest/gtest.h" #include "module_esolver/esolver_lj.h" #include "setcell.h" +#define private public +#define protected public +#include "module_md/langevin.h" + #define doublethreshold 1e-12 /************************************************ @@ -37,7 +37,7 @@ class Langevin_test : public testing::Test { protected: - MD_base *mdrun; + MD_base* mdrun; UnitCell ucell; void SetUp() @@ -45,7 +45,7 @@ class Langevin_test : public testing::Test Setcell::setupcell(ucell); Setcell::parameters(); - ModuleESolver::ESolver *p_esolver = new ModuleESolver::ESolver_LJ(); + ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); p_esolver->Init(INPUT, ucell); mdrun = new Langevin(INPUT.mdp, ucell); @@ -106,7 +106,8 @@ TEST_F(Langevin_test, first_half) TEST_F(Langevin_test, second_half) { mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); - mdrun->second_half(GlobalV::MY_RANK);; + mdrun->second_half(GlobalV::MY_RANK); + ; EXPECT_NEAR(mdrun->pos[0].x, -0.00066954020090275205, doublethreshold); EXPECT_NEAR(mdrun->pos[0].y, 3.3862365219131354e-05, doublethreshold); diff --git a/source/module_md/test/msst_test.cpp b/source/module_md/test/msst_test.cpp index 47631d2a96..3ad126f114 100644 --- a/source/module_md/test/msst_test.cpp +++ b/source/module_md/test/msst_test.cpp @@ -1,12 +1,12 @@ -#define private public -#define protected public -#include "module_md/msst.h" - #include "gmock/gmock.h" #include "gtest/gtest.h" #include "module_esolver/esolver_lj.h" #include "setcell.h" +#define private public +#define protected public +#include "module_md/msst.h" + #define doublethreshold 1e-12 /************************************************ @@ -37,7 +37,7 @@ class MSST_test : public testing::Test { protected: - MD_base *mdrun; + MD_base* mdrun; UnitCell ucell; void SetUp() @@ -45,7 +45,7 @@ class MSST_test : public testing::Test Setcell::setupcell(ucell); Setcell::parameters(); - ModuleESolver::ESolver *p_esolver = new ModuleESolver::ESolver_LJ(); + ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); p_esolver->Init(INPUT, ucell); mdrun = new MSST(INPUT.mdp, ucell); @@ -131,7 +131,8 @@ TEST_F(MSST_test, first_half) TEST_F(MSST_test, second_half) { mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); - mdrun->second_half(GlobalV::MY_RANK);; + mdrun->second_half(GlobalV::MY_RANK); + ; EXPECT_NEAR(ucell.lat0, 1.0, doublethreshold); EXPECT_NEAR(ucell.lat0_angstrom, 0.52917700000000001, doublethreshold); @@ -201,7 +202,7 @@ TEST_F(MSST_test, restart) mdrun->restart(GlobalV::MY_RANK, GlobalV::global_readin_dir); remove("Restart_md.dat"); - MSST *msst = dynamic_cast(mdrun); + MSST* msst = dynamic_cast(mdrun); EXPECT_EQ(mdrun->step_rst_, 3); EXPECT_EQ(msst->omega[mdrun->mdp.msst_direction], -0.00977662); EXPECT_EQ(msst->e0, -0.00768262); diff --git a/source/module_md/test/nhchain_test.cpp b/source/module_md/test/nhchain_test.cpp index 23acd65cbc..a607109dfe 100644 --- a/source/module_md/test/nhchain_test.cpp +++ b/source/module_md/test/nhchain_test.cpp @@ -1,12 +1,12 @@ -#define private public -#define protected public -#include "module_md/nhchain.h" - #include "gmock/gmock.h" #include "gtest/gtest.h" #include "module_esolver/esolver_lj.h" #include "setcell.h" +#define private public +#define protected public +#include "module_md/nhchain.h" + #define doublethreshold 1e-12 /************************************************ @@ -37,7 +37,7 @@ class NHC_test : public testing::Test { protected: - MD_base *mdrun; + MD_base* mdrun; UnitCell ucell; void SetUp() @@ -45,7 +45,7 @@ class NHC_test : public testing::Test Setcell::setupcell(ucell); Setcell::parameters(); - ModuleESolver::ESolver *p_esolver = new ModuleESolver::ESolver_LJ(); + ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); p_esolver->Init(INPUT, ucell); INPUT.mdp.md_type = "npt"; @@ -108,7 +108,8 @@ TEST_F(NHC_test, first_half) TEST_F(NHC_test, second_half) { mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); - mdrun->second_half(GlobalV::MY_RANK);; + mdrun->second_half(GlobalV::MY_RANK); + ; EXPECT_NEAR(mdrun->pos[0].x, -0.00023793503786683287, doublethreshold); EXPECT_NEAR(mdrun->pos[0].y, 0.0001777972998948069, doublethreshold); @@ -140,7 +141,8 @@ TEST_F(NHC_test, second_half) TEST_F(NHC_test, write_restart) { mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); - mdrun->second_half(GlobalV::MY_RANK);; + mdrun->second_half(GlobalV::MY_RANK); + ; mdrun->step_ = 1; mdrun->step_rst_ = 2; mdrun->write_restart(GlobalV::MY_RANK, GlobalV::global_out_dir); @@ -173,7 +175,7 @@ TEST_F(NHC_test, restart) mdrun->restart(GlobalV::MY_RANK, GlobalV::global_readin_dir); remove("Restart_md.dat"); - Nose_Hoover *nhc = dynamic_cast(mdrun); + Nose_Hoover* nhc = dynamic_cast(mdrun); EXPECT_EQ(mdrun->step_rst_, 3); EXPECT_EQ(mdrun->mdp.md_tchain, 4); EXPECT_EQ(mdrun->mdp.md_pchain, 4); diff --git a/source/module_md/test/verlet_test.cpp b/source/module_md/test/verlet_test.cpp index 9c6b8f0324..fcea61103e 100644 --- a/source/module_md/test/verlet_test.cpp +++ b/source/module_md/test/verlet_test.cpp @@ -1,12 +1,12 @@ -#define private public -#define protected public -#include "module_md/verlet.h" - #include "gmock/gmock.h" #include "gtest/gtest.h" #include "module_esolver/esolver_lj.h" #include "setcell.h" +#define private public +#define protected public +#include "module_md/verlet.h" + #define doublethreshold 1e-12 /************************************************ @@ -37,7 +37,7 @@ class Verlet_test : public testing::Test { protected: - MD_base *mdrun; + MD_base* mdrun; UnitCell ucell; void SetUp() @@ -45,7 +45,7 @@ class Verlet_test : public testing::Test Setcell::setupcell(ucell); Setcell::parameters(); - ModuleESolver::ESolver *p_esolver = new ModuleESolver::ESolver_LJ(); + ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); p_esolver->Init(INPUT, ucell); mdrun = new Verlet(INPUT.mdp, ucell); @@ -107,7 +107,8 @@ TEST_F(Verlet_test, NVE) { mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); mdrun->mdp.md_type = "nve"; - mdrun->second_half(GlobalV::MY_RANK);; + mdrun->second_half(GlobalV::MY_RANK); + ; EXPECT_NEAR(mdrun->pos[0].x, -0.00054545529007222658, doublethreshold); EXPECT_NEAR(mdrun->pos[0].y, 0.00029590658162135359, doublethreshold); @@ -141,7 +142,8 @@ TEST_F(Verlet_test, Anderson) mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); mdrun->mdp.md_type = "nvt"; mdrun->mdp.md_thermostat = "anderson"; - mdrun->second_half(GlobalV::MY_RANK);; + mdrun->second_half(GlobalV::MY_RANK); + ; EXPECT_NEAR(mdrun->pos[0].x, -0.00054545529007222658, doublethreshold); EXPECT_NEAR(mdrun->pos[0].y, 0.00029590658162135359, doublethreshold); @@ -175,7 +177,8 @@ TEST_F(Verlet_test, Berendsen) mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); mdrun->mdp.md_type = "nvt"; mdrun->mdp.md_thermostat = "berendsen"; - mdrun->second_half(GlobalV::MY_RANK);; + mdrun->second_half(GlobalV::MY_RANK); + ; EXPECT_NEAR(mdrun->pos[0].x, -0.00054545529007222658, doublethreshold); EXPECT_NEAR(mdrun->pos[0].y, 0.00029590658162135359, doublethreshold); @@ -209,7 +212,8 @@ TEST_F(Verlet_test, rescaling) mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); mdrun->mdp.md_type = "nvt"; mdrun->mdp.md_thermostat = "rescaling"; - mdrun->second_half(GlobalV::MY_RANK);; + mdrun->second_half(GlobalV::MY_RANK); + ; EXPECT_NEAR(mdrun->pos[0].x, -0.00054545529007222658, doublethreshold); EXPECT_NEAR(mdrun->pos[0].y, 0.00029590658162135359, doublethreshold); @@ -243,7 +247,8 @@ TEST_F(Verlet_test, rescale_v) mdrun->first_half(GlobalV::MY_RANK, GlobalV::ofs_running); mdrun->mdp.md_type = "nvt"; mdrun->mdp.md_thermostat = "rescale_v"; - mdrun->second_half(GlobalV::MY_RANK);; + mdrun->second_half(GlobalV::MY_RANK); + ; EXPECT_NEAR(mdrun->pos[0].x, -0.00054545529007222658, doublethreshold); EXPECT_NEAR(mdrun->pos[0].y, 0.00029590658162135359, doublethreshold);