Skip to content

feat(md): add thermodynamic output frequency #8160

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
19hello wants to merge 4 commits into
base: develop
Choose a base branch
from
+95 −108
Open

feat(md): add thermodynamic output frequency #8160

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
1021a75
feat(md): add thermodynamic output frequency
Oct 16, 2026
19f64d4
refactor(md): remove per-atom force output from md log
Oct 18, 2026
32e94cd
test(md): collect MD forces from trajectory dump
Oct 18, 2026
4a0bdfe
test(ofdft): dump forces for short MD cases
Oct 18, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
60 changes: 0 additions & 60 deletions source/source_io/module_output/output_log.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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>

Expand Down Expand Up @@ -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)
{
Expand Down
8 changes: 0 additions & 8 deletions source/source_io/module_output/output_log.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,6 @@
#include "source_base/matrix.h"
#include "source_cell/unitcell.h"

class MDCell;

namespace elecstate
{
struct Efermi;
Expand Down Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion source/source_io/module_parameter/md_parameter.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Comment thread
mohanchen marked this conversation as resolved.
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

Expand Down
23 changes: 11 additions & 12 deletions source/source_io/module_parameter/read_input_item_md.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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";
Expand All @@ -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";
Expand Down
2 changes: 1 addition & 1 deletion source/source_io/test/read_input_ptest.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -402,6 +402,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);
Expand All @@ -415,7 +416,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);
Expand Down
2 changes: 1 addition & 1 deletion source/source_io/test/support/INPUT
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
6 changes: 6 additions & 0 deletions source/source_io/test_serial/read_input_item_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2231,6 +2231,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);

@ABACUS-PR-BOT ABACUS-PR-BOT Oct 19, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Non-blocking: this verifies parsing but does not exercise the new output cadence. Please add output-level coverage for 0, 1, and a larger frequency, including startup and restart, checking the selected banner/scalar records and that forces and integration state stay unchanged. The current run_md_test.cpp only exercises cell preparation.

}
{ // md_prec_level
auto it = find_label("md_prec_level", readinput.input_lists);
TestParameters::input(param).mdp.md_prec_level = 1;
Expand Down
2 changes: 1 addition & 1 deletion source/source_md/md_base.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@ void MD_base::setup(ModuleESolver::ESolver* p_esolver, const std::string& global

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())
{
Expand Down
7 changes: 1 addition & 6 deletions source/source_md/md_func.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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");
Expand All @@ -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
Expand Down
3 changes: 1 addition & 2 deletions source/source_md/md_func.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
*
Expand Down
19 changes: 14 additions & 5 deletions source/source_md/run_md.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
{
const int md_step = mdrun->step_ + mdrun->step_rst_;
const int thermo_step = mdrun->step_ == 0 ? md_step : md_step + 1;
const 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.
Expand Down Expand Up @@ -235,8 +239,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);

@ABACUS-PR-BOT ABACUS-PR-BOT Oct 19, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[P2] Apply the disabled setting to the setup banner too. With md_thermo_freq=0, mdrun->setup() still reaches the unconditional ModuleIO::print_screen() in MD_base::setup, so both stdout and running_md.log contain STEP OF MOLECULAR DYNAMICS despite the documented disable setting. Please also gate that initial banner, including restart setup, and cover the zero-frequency case.

}
mdrun->first_half(GlobalV::ofs_running);

/// update force and virial due to the update of atom positions
Expand All @@ -246,8 +253,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())
{
Expand All @@ -268,7 +274,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)
{
Expand Down
6 changes: 3 additions & 3 deletions source/source_md/test/lj_pot_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}

Expand All @@ -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);
Expand Down Expand Up @@ -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);
Expand Down
2 changes: 1 addition & 1 deletion source/source_md/test/test_md_func.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
2 changes: 1 addition & 1 deletion tests/07_OFDFT/25_OF_MD_LDA/INPUT
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
2 changes: 1 addition & 1 deletion tests/07_OFDFT/26_OF_MD_LibxcPBE/INPUT
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion tests/07_OFDFT/30_TDOFDFT_Al/INPUT
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
2 changes: 1 addition & 1 deletion tests/07_OFDFT/31_TDOFDFT_Al_CD/INPUT
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
2 changes: 1 addition & 1 deletion tests/07_OFDFT/32_TDOFDFT_Al_mCD/INPUT
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Loading
Loading