feat(md): add thermodynamic output frequency

source/source_io/module_output/output_log.cpp

@@ -7,8 +7,6 @@
 #include "source_base/global_variable.h"
 #include "source_base/parallel_reduce.h"
 #include "source_base/parallel_comm.h"
-#include "source_cell/mdcell.h"
-
 #include <cmath>
 #include <iomanip>
 
@@ -286,64 +284,6 @@ void print_force(std::ofstream& ofs_running,
     }
 }
 
-void print_force(std::ofstream& ofs, const MDCell& cell, const std::string& name)
-{
-    const double output_acc = 1.0e-8;
-    const double force_unit = ModuleBase::Hartree_to_eV / ModuleBase::BOHR_TO_A;
-    const std::vector<LocalAtom>& owned_atoms = cell.owned_atoms();
-    const std::vector<std::string>& type_labels = cell.type_labels();
-
-    const auto print_atom = [&ofs, &type_labels, output_acc, force_unit](const LocalAtom& atom) {
-        const std::string& label = type_labels[static_cast<std::size_t>(atom.type)];
-        const double fx = std::abs(atom.force.x) > output_acc ? atom.force.x * force_unit : 0.0;
-        const double fy = std::abs(atom.force.y) > output_acc ? atom.force.y * force_unit : 0.0;
-        const double fz = std::abs(atom.force.z) > output_acc ? atom.force.z * force_unit : 0.0;
-        ofs << std::setw(9) << label + std::to_string(atom.type_index + 1)
-            << std::setw(20) << std::fixed << std::setprecision(10) << fx
-            << std::setw(20) << fy << std::setw(20) << fz << std::endl;
-    };
-
-#ifdef __MPI
-    int rank = 0;
-    int size = 1;
-    MPI_Comm_rank(cell.communicator(), &rank);
-    MPI_Comm_size(cell.communicator(), &size);
-    if (rank != 0)
-    {
-        const int nowned_atoms = cell.owned_atoms().size();
-        MPI_Send(&nowned_atoms, 1, MPI_INT, 0, 0, cell.communicator());
-        for (const LocalAtom& atom : owned_atoms)
-        {
-            MPI_Send(&atom.type, 1, MPI_INT, 0, 1, cell.communicator());
-            MPI_Send(&atom.type_index, 1, MPI_INT64_T, 0, 2, cell.communicator());
-            MPI_Send(&atom.force.x, 3, MPI_DOUBLE, 0, 3, cell.communicator());
-        }
-        return;
-    }
-#endif
-
-    ofs << "\n #" << name << "#" << std::endl;
-    ofs << std::setw(9) << "Atoms" << std::setw(20) << "Force_x" << std::setw(20)
-        << "Force_y" << std::setw(20) << "Force_z" << std::endl;
-    for (const LocalAtom& atom : owned_atoms) print_atom(atom);
-
-#ifdef __MPI
-    for (int source = 1; source < size; ++source)
-    {
-        int nlocal = 0;
-        MPI_Recv(&nlocal, 1, MPI_INT, source, 0, cell.communicator(), MPI_STATUS_IGNORE);
-        for (int iat = 0; iat < nlocal; ++iat)
-        {
-            LocalAtom atom;
-            MPI_Recv(&atom.type, 1, MPI_INT, source, 1, cell.communicator(), MPI_STATUS_IGNORE);
-            MPI_Recv(&atom.type_index, 1, MPI_INT64_T, source, 2, cell.communicator(), MPI_STATUS_IGNORE);
-            MPI_Recv(&atom.force.x, 3, MPI_DOUBLE, source, 3, cell.communicator(), MPI_STATUS_IGNORE);
-            print_atom(atom);
-        }
-    }
-#endif
-}
-
 void print_stress(const std::string& name, const ModuleBase::matrix& scs, 
   const bool screen, const bool ry, std::ofstream &ofs)
 {

source/source_io/module_output/output_log.h

@@ -7,8 +7,6 @@
 #include "source_base/matrix.h"
 #include "source_cell/unitcell.h"
 
-class MDCell;
-
 namespace elecstate
 {
 struct Efermi;
@@ -67,12 +65,6 @@ void print_force(std::ofstream& ofs,
                  const ModuleBase::matrix& force,
                  bool ry = true);
 
-/// @brief output forces stored on locally owned atoms of an MDCell
-/// @param ofs the output stream on rank zero
-/// @param cell the MD cell
-/// @param name force term name
-void print_force(std::ofstream& ofs, const MDCell& cell, const std::string& name);
-
 /// @brief output stress components
 /// @param name stress term name
 /// @param f stress components

source/source_io/module_parameter/md_parameter.h

@@ -19,9 +19,9 @@ struct MD_para
     double md_tfirst = -1.0;           ///< Temperature (in Hartree, 1 Hartree ~ 3E5 K)
     double md_tlast = -1.0;            ///< Target temperature
     int md_dumpfreq = 1;               ///< The period to dump MD information
+    int md_thermo_freq = 1;            ///< The period to output MD thermodynamic information
     double md_neighbor_skin = 0.0;     ///< Extra neighbor-list radius in Angstrom for MDCell solvers
     int md_restartfreq = 5;            ///< The period to output MD restart information
-    bool md_out_force = true;          ///< output all atomic forces into running_md.log
     int md_seed = -1;                  ///< random seed for MD
     int md_prec_level = 0;             ///< precision level for vc-md
 

source/source_io/module_parameter/read_input_item_md.cpp

@@ -186,6 +186,17 @@ Note that md_tlast is only used in NVT/NPT simulations. If md_tlast is unset or
         read_sync_int(input.mdp.md_dumpfreq);
         this->add_item(item);
     }
+    {
+        Input_Item item("md_thermo_freq");
+        item.annotation = "The period to output MD thermodynamic information";
+        item.category = "Molecular dynamics";
+        item.type = "Integer";
+        item.description = "The output frequency of MD step, energy, temperature, and pressure information in OUT.${suffix}/running_md.log and on screen. Set to 0 to disable this output.";
+        item.default_value = "1";
+        item.unit = "";
+        read_sync_int(input.mdp.md_thermo_freq);
+        this->add_item(item);
+    }
     {
         Input_Item item("md_neighbor_skin");
         item.annotation = "extra MDCell neighbor-list radius";
@@ -200,18 +211,6 @@ Note that md_tlast is only used in NVT/NPT simulations. If md_tlast is unset or
         read_sync_double(input.mdp.md_neighbor_skin);
         this->add_item(item);
     }
-    {
-        Input_Item item("md_out_force");
-        item.annotation = "output all atomic forces into running_md.log or not";
-        item.category = "Molecular dynamics";
-        item.type = "Boolean";
-        item.description = "Whether to output the TOTAL-FORCE table in OUT.${suffix}/running_md.log for MDCell molecular dynamics. This does not affect force calculation or molecular dynamics integration.";
-        item.default_value = "True";
-        item.unit = "";
-        item.set_availability("");
-        read_sync_bool(input.mdp.md_out_force);
-        this->add_item(item);
-    }
     {
         Input_Item item("dump_force");
         item.annotation = "output atomic forces into the file MD_dump or not";

source/source_io/test/read_input_ptest.cpp

@@ -403,6 +403,7 @@ TEST_F(InputParaTest, ParaRead)
     EXPECT_EQ(param.inp.mdp.md_damp, 1);
     EXPECT_EQ(param.inp.mdp.md_dt, 1);
     EXPECT_EQ(param.inp.mdp.md_dumpfreq, 1);
+    EXPECT_EQ(param.inp.mdp.md_thermo_freq, 1);
     EXPECT_EQ(param.inp.mdp.md_nraise, 1);
     EXPECT_EQ(param.inp.cal_syns[0], 0);
     EXPECT_EQ(param.inp.cal_syns[1], 8);
@@ -416,7 +417,6 @@ TEST_F(InputParaTest, ParaRead)
     EXPECT_EQ(param.inp.mdp.md_pmode, "iso");
     EXPECT_EQ(param.inp.mdp.md_restart, 0);
     EXPECT_EQ(param.inp.mdp.md_restartfreq, 5);
-    EXPECT_FALSE(param.inp.mdp.md_out_force);
     EXPECT_FALSE(param.inp.mdp.plumed);
     EXPECT_EQ(param.inp.mdp.plumed_file, "plumed.dat");
     EXPECT_EQ(param.inp.mdp.md_seed, -1);

source/source_io/test/support/INPUT

@@ -186,8 +186,8 @@ md_tchain                      1 #number of Nose-Hoover chains
 md_tfirst                      -1 #temperature first
 md_tlast                       -1 #temperature last
 md_dumpfreq                    1 #The period to dump MD information
+md_thermo_freq                 1 #The period to output MD thermodynamic information
 md_restartfreq                 5 #The period to output MD restart information
-md_out_force                   0 #output all atomic forces into running_md.log or not
 md_seed                        -1 #random seed for MD
 md_prec_level                  2 #precision level for vc-md
 ref_cell_factor                1.0 #construct a reference cell bigger than the initial cell (currently disabled, must be 1.0)

source/source_io/test_serial/read_input_item_test.cpp

@@ -2364,6 +2364,12 @@ TEST_F(InputTest, Item_test2)
         it->second.reset_value(it->second, param);
         EXPECT_EQ(TestParameters::input(param).mdp.md_nstep, 50);
     }
+    { // md_thermo_freq
+        auto it = find_label("md_thermo_freq", readinput.input_lists);
+        it->second.str_values = {"100"};
+        it->second.read_value(it->second, param);
+        EXPECT_EQ(TestParameters::input(param).mdp.md_thermo_freq, 100);
+    }
     { // md_prec_level
         auto it = find_label("md_prec_level", readinput.input_lists);
         TestParameters::input(param).mdp.md_prec_level = 1;

source/source_md/md_base.cpp

@@ -4,7 +4,6 @@
 #ifdef __MPI
 #include "mpi.h"
 #endif
-#include "source_io/module_output/print_info.h"
 #include <algorithm>
 #include <iomanip>
 
@@ -54,14 +53,7 @@ void MD_base::setup(ModuleESolver::ESolver* p_esolver, const std::string& global
         restart(global_readin_dir);
     }
 
-    // mohan add 2026-01-04
-    const int stress_step = 0;
-    const int force_step = 0;
-    const int istep_print = step_ + step_rst_ + 1;
-
-	ModuleIO::print_screen(stress_step, force_step, istep_print);
-
-    MD_func::force_virial(p_esolver, step_, mdcell, decomp, potential, cal_stress, virial, mdp.md_out_force);
+    MD_func::force_virial(p_esolver, step_, mdcell, decomp, potential, cal_stress, virial);
     MD_func::compute_stress(mdcell, cal_stress, virial, stress);
     if (mdcell.has_backing_unitcell())
     {

source/source_md/md_func.cpp

@@ -329,8 +329,7 @@ void force_virial(ModuleESolver::ESolver* p_esolver,
                   DomainDecomposition& decomp,
                   double& potential,
                   const bool& cal_stress,
-                  ModuleBase::matrix& virial,
-                  const bool& md_out_force)
+                  ModuleBase::matrix& virial)
 {
     ModuleBase::TITLE("MD_func", "force_virial");
     ModuleBase::timer::start("MD_func", "force_virial");
@@ -341,10 +340,6 @@ void force_virial(ModuleESolver::ESolver* p_esolver,
         decomp.accumulate_ghost_forces(mdcell);
         potential = 0.5 * p_esolver->cal_energy();
         for (LocalAtom& atom : mdcell.owned_atoms()) atom.force *= 0.5;
-        if (md_out_force)
-        {
-            ModuleIO::print_force(GlobalV::ofs_running, mdcell, "TOTAL-FORCE (eV/Angstrom)");
-        }
         if (cal_stress) { p_esolver->cal_stress(static_cast<BaseCell&>(mdcell), virial); virial *= 0.5; }
     }
     else

source/source_md/md_func.h

@@ -95,8 +95,7 @@ void force_virial(ModuleESolver::ESolver* p_esolver,
                   DomainDecomposition& decomp,
                   double& potential,
                   const bool& cal_stress,
-                  ModuleBase::matrix& virial,
-                  const bool& md_out_force);
+                  ModuleBase::matrix& virial);
 /**
  * @brief calculate the ionic kinetic energy
  *

source/source_md/run_md.cpp

@@ -191,6 +191,10 @@ void md_line(MDCell& mdcell,
     /// md cycle, mohan update 2026-01-04, change '<=' to '<'
     while ((mdrun->step_ + mdrun->step_rst_) < param_in.mdp.md_nstep && !mdrun->stop)
     {
+        int md_step = mdrun->step_ + mdrun->step_rst_;
+        int thermo_step = mdrun->step_ == 0 ? md_step : md_step + 1;
+        bool print_thermo = param_in.mdp.md_thermo_freq > 0
+                            && thermo_step % param_in.mdp.md_thermo_freq == 0;
 #ifdef __JSON
         // JSON output currently follows the UnitCell-backed electronic-structure path.
         // Start one output record before the solver appends SCF information for this MD step.
@@ -204,6 +208,20 @@ void md_line(MDCell& mdcell,
         {
             mdrun->setup(p_esolver, param_in.globalv.global_readin_dir, decomp);
 
+            // setup() may load step_rst_ from Restart_md.txt.
+            md_step = mdrun->step_ + mdrun->step_rst_;
+            thermo_step = md_step;
+            print_thermo = param_in.mdp.md_thermo_freq > 0
+                           && thermo_step % param_in.mdp.md_thermo_freq == 0;
+
+            if (print_thermo)
+            {
+                const int stress_step = 0;
+                const int force_step = 0;
+                const int istep_print = md_step + 1;
+                ModuleIO::print_screen(stress_step, force_step, istep_print);
+            }
+
             if (param_in.mdp.plumed)
             {
                 int nproc = 1;
@@ -235,8 +253,11 @@ void md_line(MDCell& mdcell,
             // mohan add 2026-01-04
             const int stress_step = 0;
             const int force_step = 0;
-            const int istep_print = mdrun->step_ + mdrun->step_rst_ + 1;
-            ModuleIO::print_screen(stress_step, force_step, istep_print);
+            const int istep_print = thermo_step;
+            if (print_thermo)
+            {
+                ModuleIO::print_screen(stress_step, force_step, istep_print);
+            }
             mdrun->first_half(GlobalV::ofs_running);
 
             /// update force and virial due to the update of atom positions
@@ -246,8 +267,7 @@ void md_line(MDCell& mdcell,
                                   decomp,
                                   mdrun->potential,
                                   param_in.inp.cal_stress,
-                                  mdrun->virial,
-                                  param_in.mdp.md_out_force);
+                                  mdrun->virial);
 
             if (plumed.active())
             {
@@ -268,7 +288,10 @@ void md_line(MDCell& mdcell,
                                                      mdrun->frozen_freedom_);
         }
 
-        mdrun->print_md(GlobalV::ofs_running, PARAM.inp.cal_stress);
+        if (print_thermo)
+        {
+            mdrun->print_md(GlobalV::ofs_running, param_in.inp.cal_stress);
+        }
         if (param_in.mdp.md_dumpfreq > 0
             && (mdrun->step_ + mdrun->step_rst_) % param_in.mdp.md_dumpfreq == 0)
         {

source/source_md/test/lj_pot_test.cpp

@@ -48,7 +48,7 @@ TEST_F(LJ_pot_test, potential)
     DomainDecomposition decomp;
     decomp.init(ModuleBase::world_comm_domain(), mdcell.latvec(), mdcell.lat0(), 0.0, 0.0);
     p_esolver->before_all_runners(mdcell, inp);
-    MD_func::force_virial(p_esolver, 0, mdcell, decomp, potential, true, stress, false);
+    MD_func::force_virial(p_esolver, 0, mdcell, decomp, potential, true, stress);
     EXPECT_NEAR(potential, -0.011957818623534381, doublethreshold);
 }
 
@@ -73,7 +73,7 @@ TEST_F(LJ_pot_test, force)
     DomainDecomposition decomp;
     decomp.init(ModuleBase::world_comm_domain(), mdcell.latvec(), mdcell.lat0(), 0.0, 0.0);
     p_esolver->before_all_runners(mdcell, inp);
-    MD_func::force_virial(p_esolver, 0, mdcell, decomp, potential, true, stress, false);
+    MD_func::force_virial(p_esolver, 0, mdcell, decomp, potential, true, stress);
     const std::vector<LocalAtom>& atoms = mdcell.owned_atoms();
     EXPECT_NEAR(atoms[0].force.x, 0.00049817733089377704, doublethreshold);
     EXPECT_NEAR(atoms[0].force.y, 0.00082237246837022328, doublethreshold);
@@ -119,7 +119,7 @@ TEST_F(LJ_pot_test, stress)
     DomainDecomposition decomp;
     decomp.init(ModuleBase::world_comm_domain(), mdcell.latvec(), mdcell.lat0(), 0.0, 0.0);
     p_esolver->before_all_runners(mdcell, inp);
-    MD_func::force_virial(p_esolver, 0, mdcell, decomp, potential, true, stress, false);
+    MD_func::force_virial(p_esolver, 0, mdcell, decomp, potential, true, stress);
     EXPECT_NEAR(stress(0, 0), 8.0360222227631859e-07, doublethreshold);
     EXPECT_NEAR(stress(0, 1), 1.7207745586539077e-07, doublethreshold);
     EXPECT_NEAR(stress(0, 2), 0.0, doublethreshold);

source/source_md/test/test_md_func.cpp

@@ -64,7 +64,7 @@ TEST(MdForceVirialTest, ReturnsGhostForcesBeforeConvertingUnitsExactlyOnce)
     ModuleBase::matrix virial(3, 3);
     for (int step = 0; step < 2; ++step)
     {
-        MD_func::force_virial(&solver, step, cell, decomp, potential, true, virial, false);
+        MD_func::force_virial(&solver, step, cell, decomp, potential, true, virial);
         EXPECT_DOUBLE_EQ(potential, 3.0);
         EXPECT_DOUBLE_EQ(virial(0, 0), 4.0);
         ASSERT_EQ(cell.owned_atoms().size(), 1);

tests/07_OFDFT/25_OF_MD_LDA/INPUT

@@ -27,7 +27,7 @@ md_type           nvt
 md_nstep          3
 md_dt             0.25
 md_tfirst         58022.52706
-md_dumpfreq       10
+md_dumpfreq       1
 md_tfreq          1.08
 md_tchain         1
 nbspline          10 # be sure fft dimension is odd

tests/07_OFDFT/26_OF_MD_LibxcPBE/INPUT

@@ -27,7 +27,7 @@ md_type           nvt
 md_nstep          3
 md_dt             0.25
 md_tfirst         58022.52706
-md_dumpfreq       10
+md_dumpfreq       1
 md_tfreq          1.08
 md_tchain         1
 nbspline          10 # be sure fft dimension is odd

tests/07_OFDFT/30_TDOFDFT_Al/INPUT

@@ -30,7 +30,7 @@ md_type           nve
 md_nstep          3
 md_dt             0.25
 md_tfirst         58022.52706
-md_dumpfreq       10
+md_dumpfreq       1
 md_tfreq          1.08
 md_tchain         1
 nbspline          10 # be sure fft dimension is odd

tests/07_OFDFT/31_TDOFDFT_Al_CD/INPUT

@@ -31,7 +31,7 @@ md_type           nve
 md_nstep          3
 md_dt             0.25
 md_tfirst         58022.52706
-md_dumpfreq       10
+md_dumpfreq       1
 md_tfreq          1.08
 md_tchain         1
 nbspline          10 # be sure fft dimension is odd

tests/07_OFDFT/32_TDOFDFT_Al_mCD/INPUT

@@ -31,7 +31,7 @@ md_type           nve
 md_nstep          3
 md_dt             0.25
 md_tfirst         58022.52706
-md_dumpfreq       10
+md_dumpfreq       1
 md_tfreq          1.08
 md_tchain         1
 nbspline          10 # be sure fft dimension is odd

tests/integrate/validation_tools/catch_properties.sh

@@ -26,7 +26,7 @@ source "$PROPS_SCRIPT_DIR/props_deepks.sh"
 # trailing totaltimeref entry. Add a new category by dropping a
 # props_<cat>.sh next to the others and calling run_<cat>_props here.
 run_lr_force_props
-run_basic_props
+run_basic_props || exit 1
 
 #-------------------------------
 # matrix collectors: out_dm1 (props_mat.sh)