Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
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
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
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
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");

@Critsium-xy Critsium-xy Oct 26, 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.

Removing md_out_force (shipped in v3.11.0-beta9..11) makes inputs containing it abort, and the MDCell path no longer writes forces to running_md.log. Please state this and the dump_force/md_dumpfreq replacement in the PR description, together with the new md_thermo_freq.

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
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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);
Expand Down
2 changes: 1 addition & 1 deletion source/source_io/test/support/INPUT
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
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Comment thread
19hello marked this conversation as resolved.
}
{ // md_prec_level
auto it = find_label("md_prec_level", readinput.input_lists);
TestParameters::input(param).mdp.md_prec_level = 1;
Expand Down
10 changes: 1 addition & 9 deletions source/source_md/md_base.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@
#ifdef __MPI
#include "mpi.h"
#endif
#include "source_io/module_output/print_info.h"
#include <algorithm>
#include <iomanip>

Expand Down Expand Up @@ -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())
{
Expand Down
7 changes: 1 addition & 6 deletions source/source_md/md_func.cpp
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
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
33 changes: 28 additions & 5 deletions source/source_md/run_md.cpp
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)
{
int md_step = mdrun->step_ + mdrun->step_rst_;
int thermo_step = mdrun->step_ == 0 ? md_step : md_step + 1;

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.

The gating index is inconsistent: step 0 is tested with the 0-based md_step, later steps with md_step + 1.

  • With md_thermo_freq=10, thermo is printed at steps 0, 9, 19, 29, ..., while md_dumpfreq/md_restartfreq use md_step % freq (0, 10, 20, ...), so thermo records do not line up with MD_dump frames.
  • On restart from step_rst_=10, the first resumed step is tested with 10 (printed), whereas the same step in an uninterrupted run is tested with 11 (not printed); output depends on whether the run was restarted.

Suggest using md_step % md_thermo_freq == 0 uniformly, consistent with the other MD frequencies.

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 All @@ -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)

@Critsium-xy Critsium-xy Oct 26, 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.

Ordering regression: MD_base::setup() used to print the step banner before force_virial(); now it is printed after setup() returns, so for DFT MD the first step's SCF output precedes STEP OF MOLECULAR DYNAMICS: 1. Minimal fix: keep print_screen at its original place in setup() (after restart()) and gate it with mdp.md_thermo_freq > 0 && (step_ + step_rst_) % mdp.md_thermo_freq == 0; mdp is already a member, and the recomputation block here can then be dropped.

{
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;
Expand Down Expand Up @@ -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);
Comment thread
19hello marked this conversation as resolved.
}
mdrun->first_half(GlobalV::ofs_running);

/// update force and virial due to the update of atom positions
Expand All @@ -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())
{
Expand All @@ -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)
{
Expand Down
6 changes: 3 additions & 3 deletions source/source_md/test/lj_pot_test.cpp
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
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
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
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
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
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
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/integrate/validation_tools/catch_properties.sh
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Loading
Loading