Refactor: module_md

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);
+        Run_MD::md_line(GlobalC::ucell, p_esolver, INPUT.mdp);
     }
     else // scf; cell relaxation; nscf; etc
     {

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;
 

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

source/module_hamilt_general/module_vdw/vdwd3_parameters.cpp

@@ -6,6 +6,8 @@
 
 #include "vdwd3_parameters.h"
 
+#include "module_base/constants.h"
+
 namespace vdw
 {
 

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"
@@ -189,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 +864,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 +1449,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
@@ -2791,9 +2789,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
@@ -2948,6 +2943,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);

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 <fstream>
+#include <sstream>
 #include <string>
 #include <vector>
 
+#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

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;
@@ -531,7 +532,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;
 
     //----------------------------------------------------------

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");

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);
@@ -355,6 +352,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)
@@ -417,10 +417,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);
@@ -696,6 +693,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)

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);
@@ -358,6 +355,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);

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
 {
@@ -275,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");

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;
@@ -23,58 +23,58 @@ 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)
 {
     ModuleBase::TITLE("FIRE", "setup");
     ModuleBase::timer::tick("FIRE", "setup");
 
-    MDrun::setup(p_ensolve);
+    MD_base::setup(p_esolver, my_rank, global_readin_dir);
 
     check_force();
 
     ModuleBase::timer::tick("FIRE", "setup");
 }
 
-void FIRE::first_half()
+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);
+    MD_base::update_vel(force, my_rank);
 
     check_FIRE();
 
-    MDrun::update_pos();
+    MD_base::update_pos(my_rank);
 
     ModuleBase::timer::tick("FIRE", "first_half");
 }
 
-void FIRE::second_half()
+void FIRE::second_half(const int &my_rank)
 {
     ModuleBase::TITLE("FIRE", "second_half");
     ModuleBase::timer::tick("FIRE", "second_half");
 
-    MDrun::update_vel(force);
+    MD_base::update_vel(force, my_rank);
 
     check_force();
 
     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);
+    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;
 }
 
-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 +89,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)
@@ -144,7 +144,7 @@ void FIRE::check_force()
         }
     }
 
-    if (2.0 * max < GlobalV::FORCE_THR)
+    if (2.0 * max < mdp.force_thr)
     {
         stop = true;
     }