Refactor MD setup for pre-parallel cleanup

source/source_esolver/esolver_dp.cpp

@@ -36,6 +36,8 @@ void ESolver_DP::before_all_runners(UnitCell& ucell, const Input_para& inp)
     dp_potential = 0;
     dp_force.create(ucell.nat, 3);
     dp_virial.create(3, 3);
+    dp_cell.resize(9);
+    dp_model_virial.clear();
 
     ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU.cif", 
                                ucell, 
@@ -59,16 +61,15 @@ void ESolver_DP::runner(UnitCell& ucell, const int istep)
     ModuleBase::TITLE("ESolver_DP", "runner");
     ModuleBase::timer::start("ESolver_DP", "runner");
 
-    std::vector<double> cell(9, 0.0);
-    cell[0] = ucell.latvec.e11 * ucell.lat0_angstrom;
-    cell[1] = ucell.latvec.e12 * ucell.lat0_angstrom;
-    cell[2] = ucell.latvec.e13 * ucell.lat0_angstrom;
-    cell[3] = ucell.latvec.e21 * ucell.lat0_angstrom;
-    cell[4] = ucell.latvec.e22 * ucell.lat0_angstrom;
-    cell[5] = ucell.latvec.e23 * ucell.lat0_angstrom;
-    cell[6] = ucell.latvec.e31 * ucell.lat0_angstrom;
-    cell[7] = ucell.latvec.e32 * ucell.lat0_angstrom;
-    cell[8] = ucell.latvec.e33 * ucell.lat0_angstrom;
+    dp_cell[0] = ucell.latvec.e11 * ucell.lat0_angstrom;
+    dp_cell[1] = ucell.latvec.e12 * ucell.lat0_angstrom;
+    dp_cell[2] = ucell.latvec.e13 * ucell.lat0_angstrom;
+    dp_cell[3] = ucell.latvec.e21 * ucell.lat0_angstrom;
+    dp_cell[4] = ucell.latvec.e22 * ucell.lat0_angstrom;
+    dp_cell[5] = ucell.latvec.e23 * ucell.lat0_angstrom;
+    dp_cell[6] = ucell.latvec.e31 * ucell.lat0_angstrom;
+    dp_cell[7] = ucell.latvec.e32 * ucell.lat0_angstrom;
+    dp_cell[8] = ucell.latvec.e33 * ucell.lat0_angstrom;
 
     std::vector<double> coord(3 * ucell.nat, 0.0);
     int iat = 0;
@@ -85,12 +86,13 @@ void ESolver_DP::runner(UnitCell& ucell, const int istep)
     assert(ucell.nat == iat);
 
 #ifdef __DPMD
-    std::vector<double> f, v;
     dp_potential = 0;
     dp_force.zero_out();
     dp_virial.zero_out();
+    std::vector<double> model_force;
+    dp_model_virial.clear();
 
-    dp.compute(dp_potential, f, v, coord, atype, cell, fparam, aparam);
+    dp.compute(dp_potential, model_force, dp_model_virial, coord, atype, dp_cell, fparam, aparam);
 
     // rescale the energy, force, and stress
     const double fact_e = rescaling / ModuleBase::Ry_to_eV;
@@ -103,16 +105,16 @@ void ESolver_DP::runner(UnitCell& ucell, const int istep)
 
     for (int i = 0; i < ucell.nat; ++i)
     {
-        dp_force(i, 0) = f[3 * i] * fact_f;
-        dp_force(i, 1) = f[3 * i + 1] * fact_f;
-        dp_force(i, 2) = f[3 * i + 2] * fact_f;
+        dp_force(i, 0) = model_force[3 * i] * fact_f;
+        dp_force(i, 1) = model_force[3 * i + 1] * fact_f;
+        dp_force(i, 2) = model_force[3 * i + 2] * fact_f;
     }
 
     for (int i = 0; i < 3; ++i)
     {
         for (int j = 0; j < 3; ++j)
         {
-            dp_virial(i, j) = v[3 * i + j] * fact_v;
+            dp_virial(i, j) = dp_model_virial[3 * i + j] * fact_v;
         }
     }
 #else

source/source_esolver/esolver_dp.h

@@ -115,6 +115,8 @@ class ESolver_DP : public ESolver
     double dp_potential = 0.0;       ///< computed potential energy
     ModuleBase::matrix dp_force;     ///< computed atomic forces
     ModuleBase::matrix dp_virial;    ///< computed lattice virials
+    std::vector<double> dp_cell;     ///< DP cell buffer in Angstrom
+    std::vector<double> dp_model_virial; ///< raw virial buffer returned by DP
 };
 
 } // namespace ModuleESolver

source/source_esolver/esolver_nep.cpp

@@ -23,6 +23,7 @@
 #include "source_io/module_output/output_log.h"
 #include "source_io/module_output/cif_io.h"
 
+#include <algorithm>
 #include <numeric>
 #include <unordered_map>
 
@@ -34,6 +35,8 @@ void ESolver_NEP::before_all_runners(UnitCell& ucell, const Input_para& inp)
     nep_force.create(ucell.nat, 3);
     nep_virial.create(3, 3);
     atype.resize(ucell.nat);
+    nep_cell.resize(9);
+    nep_virial_sum.resize(9);
     _e.resize(ucell.nat);
     _f.resize(3 * ucell.nat);
     _v.resize(9 * ucell.nat);
@@ -56,16 +59,15 @@ void ESolver_NEP::runner(UnitCell& ucell, const int istep)
 
     // note that NEP are column major, thus a transpose is needed
     // cell
-    std::vector<double> cell(9, 0.0);
-    cell[0] = ucell.latvec.e11 * ucell.lat0_angstrom;
-    cell[1] = ucell.latvec.e21 * ucell.lat0_angstrom;
-    cell[2] = ucell.latvec.e31 * ucell.lat0_angstrom;
-    cell[3] = ucell.latvec.e12 * ucell.lat0_angstrom;
-    cell[4] = ucell.latvec.e22 * ucell.lat0_angstrom;
-    cell[5] = ucell.latvec.e32 * ucell.lat0_angstrom;
-    cell[6] = ucell.latvec.e13 * ucell.lat0_angstrom;
-    cell[7] = ucell.latvec.e23 * ucell.lat0_angstrom;
-    cell[8] = ucell.latvec.e33 * ucell.lat0_angstrom;
+    nep_cell[0] = ucell.latvec.e11 * ucell.lat0_angstrom;
+    nep_cell[1] = ucell.latvec.e21 * ucell.lat0_angstrom;
+    nep_cell[2] = ucell.latvec.e31 * ucell.lat0_angstrom;
+    nep_cell[3] = ucell.latvec.e12 * ucell.lat0_angstrom;
+    nep_cell[4] = ucell.latvec.e22 * ucell.lat0_angstrom;
+    nep_cell[5] = ucell.latvec.e32 * ucell.lat0_angstrom;
+    nep_cell[6] = ucell.latvec.e13 * ucell.lat0_angstrom;
+    nep_cell[7] = ucell.latvec.e23 * ucell.lat0_angstrom;
+    nep_cell[8] = ucell.latvec.e33 * ucell.lat0_angstrom;
 
     // coord
     std::vector<double> coord(3 * ucell.nat, 0.0);
@@ -88,7 +90,7 @@ void ESolver_NEP::runner(UnitCell& ucell, const int istep)
     nep_force.zero_out();
     nep_virial.zero_out();
 
-    nep.compute(atype, cell, coord, _e, _f, _v);
+    nep.compute(atype, nep_cell, coord, _e, _f, _v);
 
     // unit conversion
     const double fact_e = 1.0 / ModuleBase::Ry_to_eV;
@@ -110,13 +112,13 @@ void ESolver_NEP::runner(UnitCell& ucell, const int istep)
     }
 
     // virial
-    std::vector<double> v_sum(9, 0.0);
+    std::fill(nep_virial_sum.begin(), nep_virial_sum.end(), 0.0);
     for (int j = 0; j < 9; ++j)
     {
         for (int i = 0; i < nat; ++i)
         {
             int index = j * nat + i;
-            v_sum[j] += _v[index];
+            nep_virial_sum[j] += _v[index];
         }
     }
 
@@ -125,7 +127,7 @@ void ESolver_NEP::runner(UnitCell& ucell, const int istep)
     {
         for (int j = 0; j < 3; ++j)
         {
-            nep_virial(i, j) = v_sum[3 * i + j] * fact_v;
+            nep_virial(i, j) = nep_virial_sum[3 * i + j] * fact_v;
         }
     }
 #else

source/source_esolver/esolver_nep.h

@@ -98,11 +98,13 @@ class ESolver_NEP : public ESolver
     double nep_potential;                ///< computed potential energy
     ModuleBase::matrix nep_force;        ///< computed atomic forces
     ModuleBase::matrix nep_virial;       ///< computed lattice virials
+    std::vector<double> nep_cell;        ///< NEP cell buffer in Angstrom, column-major
+    std::vector<double> nep_virial_sum;  ///< summed per-atom virial components
     std::vector<double> _e;              ///< temporary storage for energy computation
     std::vector<double> _f;              ///< temporary storage for force computation
     std::vector<double> _v;              ///< temporary storage for virial computation
 };
 
 } // namespace ModuleESolver
 
-#endif
+#endif

source/source_md/md_func.cpp

@@ -52,6 +52,43 @@ double kinetic_energy(const int& natom, const ModuleBase::Vector3<double>* vel,
     return ke;
 }
 
+MDKineticState calc_kinetic_state(const int& natom,
+                                  const int& frozen_freedom,
+                                  const double* allmass,
+                                  const ModuleBase::Vector3<double>* vel)
+{
+    MDKineticState state;
+    if (3 * natom == frozen_freedom)
+    {
+        return state;
+    }
+
+    state.kinetic = kinetic_energy(natom, vel, allmass);
+    state.temperature = 2 * state.kinetic / (3 * natom - frozen_freedom);
+    return state;
+}
+
+MDStressState calc_stress_state(const int& natom,
+                                const double& omega,
+                                const ModuleBase::Vector3<double>* vel,
+                                const double* allmass,
+                                const ModuleBase::matrix& virial)
+{
+    MDStressState state;
+    temp_vector(natom, vel, allmass, state.temperature_tensor);
+    state.stress.create(3, 3);
+
+    for (int i = 0; i < 3; ++i)
+    {
+        for (int j = 0; j < 3; ++j)
+        {
+            state.stress(i, j) = virial(i, j) + state.temperature_tensor(i, j) / omega;
+        }
+    }
+
+    return state;
+}
+
 void compute_stress(const UnitCell& unit_in,
                     const ModuleBase::Vector3<double>* vel,
                     const double* allmass,
@@ -61,17 +98,7 @@ void compute_stress(const UnitCell& unit_in,
 {
     if (cal_stress)
     {
-        ModuleBase::matrix t_vector;
-
-        temp_vector(unit_in.nat, vel, allmass, t_vector);
-
-        for (int i = 0; i < 3; ++i)
-        {
-            for (int j = 0; j < 3; ++j)
-            {
-                stress(i, j) = virial(i, j) + t_vector(i, j) / unit_in.omega;
-            }
-        }
+        stress = calc_stress_state(unit_in.nat, unit_in.omega, vel, allmass, virial).stress;
     }
 
     return;
@@ -463,8 +490,9 @@ double current_temp(double& kinetic,
     }
     else
     {
-        kinetic = kinetic_energy(natom, vel, allmass);
-        return 2 * kinetic / (3 * natom - frozen_freedom);
+        const MDKineticState state = calc_kinetic_state(natom, frozen_freedom, allmass, vel);
+        kinetic = state.kinetic;
+        return state.temperature;
     }
 }
 

source/source_md/md_func.h

@@ -1,6 +1,7 @@
 #ifndef MD_FUNC_H
 #define MD_FUNC_H
 
+#include "md_statistics.h"
 #include "source_esolver/esolver.h"
 
 class Parameter;
@@ -117,6 +118,14 @@ void force_virial(ModuleESolver::ESolver* p_esolver,
  */
 double kinetic_energy(const int& natom, const ModuleBase::Vector3<double>* vel, const double* allmass);
 
+/**
+ * @brief calculate kinetic energy and temperature without writing caller-owned state
+ */
+MDKineticState calc_kinetic_state(const int& natom,
+                                  const int& frozen_freedom,
+                                  const double* allmass,
+                                  const ModuleBase::Vector3<double>* vel);
+
 /**
  * @brief calculate the total stress tensor
  *
@@ -134,6 +143,15 @@ void compute_stress(const UnitCell& unit_in,
                     const ModuleBase::matrix& virial,
                     ModuleBase::matrix& stress);
 
+/**
+ * @brief calculate stress and ionic temperature tensor without writing caller-owned state
+ */
+MDStressState calc_stress_state(const int& natom,
+                                const double& omega,
+                                const ModuleBase::Vector3<double>* vel,
+                                const double* allmass,
+                                const ModuleBase::matrix& virial);
+
 /**
  * @brief output the stress information
  *

source/source_md/md_statistics.h

@@ -0,0 +1,23 @@
+#ifndef MD_STATISTICS_H
+#define MD_STATISTICS_H
+
+#include "source_base/matrix.h"
+
+namespace MD_func
+{
+
+struct MDKineticState
+{
+    double kinetic = 0.0;
+    double temperature = 0.0;
+};
+
+struct MDStressState
+{
+    ModuleBase::matrix stress;
+    ModuleBase::matrix temperature_tensor;
+};
+
+} // namespace MD_func
+
+#endif // MD_STATISTICS_H

source/source_md/run_md.cpp

@@ -12,41 +12,48 @@
 #include "verlet.h"
 #include "source_cell/update_cell.h"
 #include "source_cell/print_cell.h"
-namespace Run_MD
-{
+#include <memory>
 
-void md_line(UnitCell& unit_in, ModuleESolver::ESolver* p_esolver, const Parameter& param_in)
+namespace
+{
+std::unique_ptr<MD_base> create_md_runner(const Parameter& param_in, UnitCell& unit_in)
 {
-    ModuleBase::TITLE("Run_MD", "md_line");
-    ModuleBase::timer::start("Run_MD", "md_line");
-
-    /// determine the md_type
-    MD_base* mdrun = nullptr;
     if (param_in.mdp.md_type == "fire")
     {
-        mdrun = new FIRE(param_in, unit_in);
+        return std::unique_ptr<MD_base>(new FIRE(param_in, unit_in));
     }
-    else if ((param_in.mdp.md_type == "nvt" && param_in.mdp.md_thermostat == "nhc") || param_in.mdp.md_type == "npt")
+    if ((param_in.mdp.md_type == "nvt" && param_in.mdp.md_thermostat == "nhc") || param_in.mdp.md_type == "npt")
     {
-        mdrun = new Nose_Hoover(param_in, unit_in);
+        return std::unique_ptr<MD_base>(new Nose_Hoover(param_in, unit_in));
     }
-    else if (param_in.mdp.md_type == "nve" || param_in.mdp.md_type == "nvt")
+    if (param_in.mdp.md_type == "nve" || param_in.mdp.md_type == "nvt")
     {
-        mdrun = new Verlet(param_in, unit_in);
+        return std::unique_ptr<MD_base>(new Verlet(param_in, unit_in));
     }
-    else if (param_in.mdp.md_type == "langevin")
+    if (param_in.mdp.md_type == "langevin")
     {
-        mdrun = new Langevin(param_in, unit_in);
+        return std::unique_ptr<MD_base>(new Langevin(param_in, unit_in));
     }
-    else if (param_in.mdp.md_type == "msst")
+    if (param_in.mdp.md_type == "msst")
     {
-        mdrun = new MSST(param_in, unit_in);
-    }
-    else
-    {
-        ModuleBase::WARNING_QUIT("md_line", "no such md_type!");
+        return std::unique_ptr<MD_base>(new MSST(param_in, unit_in));
     }
 
+    ModuleBase::WARNING_QUIT("md_line", "no such md_type!");
+    return nullptr;
+}
+} // namespace
+
+namespace Run_MD
+{
+
+void md_line(UnitCell& unit_in, ModuleESolver::ESolver* p_esolver, const Parameter& param_in)
+{
+    ModuleBase::TITLE("Run_MD", "md_line");
+    ModuleBase::timer::start("Run_MD", "md_line");
+
+    std::unique_ptr<MD_base> mdrun = create_md_runner(param_in, unit_in);
+
     /// md cycle, mohan update 2026-01-04, change '<=' to '<'
     while ((mdrun->step_ + mdrun->step_rst_) < param_in.mdp.md_nstep && !mdrun->stop)
     {
@@ -129,7 +136,6 @@ void md_line(UnitCell& unit_in, ModuleESolver::ESolver* p_esolver, const Paramet
         mdrun->step_++;
     }
 
-    delete mdrun;
     ModuleBase::timer::end("Run_MD", "md_line");
     return;
 }

source/source_md/test/CMakeLists.txt

@@ -37,7 +37,6 @@ list(APPEND depend_files
   ../../source_base/realarray.cpp
   ../../source_base/complexarray.cpp
   ../../source_base/complexmatrix.cpp
-  ../../source_base/global_variable.cpp
   ../../source_base/libm/branred.cpp
   ../../source_base/libm/sincos.cpp
   ../../source_base/math_integral.cpp