diff --git a/source/source_io/module_output/output_log.cpp b/source/source_io/module_output/output_log.cpp index 8220e16f710..9abc4aae88d 100644 --- a/source/source_io/module_output/output_log.cpp +++ b/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 #include @@ -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& owned_atoms = cell.owned_atoms(); - const std::vector& 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(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) { diff --git a/source/source_io/module_output/output_log.h b/source/source_io/module_output/output_log.h index 100e6cbf222..dfafae795f9 100644 --- a/source/source_io/module_output/output_log.h +++ b/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 diff --git a/source/source_io/module_parameter/md_parameter.h b/source/source_io/module_parameter/md_parameter.h index 706218e6960..2511379d39d 100644 --- a/source/source_io/module_parameter/md_parameter.h +++ b/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 diff --git a/source/source_io/module_parameter/read_input_item_md.cpp b/source/source_io/module_parameter/read_input_item_md.cpp index d5f6dc6bd67..76da32e9c93 100644 --- a/source/source_io/module_parameter/read_input_item_md.cpp +++ b/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"; diff --git a/source/source_io/test/read_input_ptest.cpp b/source/source_io/test/read_input_ptest.cpp index 30398cd7796..104c3de7a6f 100644 --- a/source/source_io/test/read_input_ptest.cpp +++ b/source/source_io/test/read_input_ptest.cpp @@ -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); @@ -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); diff --git a/source/source_io/test/support/INPUT b/source/source_io/test/support/INPUT index 8e8d3902c9c..05fb6eb3c6c 100644 --- a/source/source_io/test/support/INPUT +++ b/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) diff --git a/source/source_io/test_serial/read_input_item_test.cpp b/source/source_io/test_serial/read_input_item_test.cpp index 6107585e733..e58a2037ed2 100644 --- a/source/source_io/test_serial/read_input_item_test.cpp +++ b/source/source_io/test_serial/read_input_item_test.cpp @@ -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); + } { // md_prec_level auto it = find_label("md_prec_level", readinput.input_lists); TestParameters::input(param).mdp.md_prec_level = 1; diff --git a/source/source_md/md_base.cpp b/source/source_md/md_base.cpp index cb00647a865..2de3e906841 100644 --- a/source/source_md/md_base.cpp +++ b/source/source_md/md_base.cpp @@ -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()) { diff --git a/source/source_md/md_func.cpp b/source/source_md/md_func.cpp index a549e87732c..10b3b2e94ba 100644 --- a/source/source_md/md_func.cpp +++ b/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(mdcell), virial); virial *= 0.5; } } else diff --git a/source/source_md/md_func.h b/source/source_md/md_func.h index c8095d069e8..0eb09f622d4 100644 --- a/source/source_md/md_func.h +++ b/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 * diff --git a/source/source_md/run_md.cpp b/source/source_md/run_md.cpp index 114f9cbb3a5..8a1ddce9edc 100644 --- a/source/source_md/run_md.cpp +++ b/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) { + 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. @@ -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); + } mdrun->first_half(GlobalV::ofs_running); /// update force and virial due to the update of atom positions @@ -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()) { @@ -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) { diff --git a/source/source_md/test/lj_pot_test.cpp b/source/source_md/test/lj_pot_test.cpp index 0ae49fb2b41..555bc8dfe6c 100644 --- a/source/source_md/test/lj_pot_test.cpp +++ b/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& 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); diff --git a/source/source_md/test/test_md_func.cpp b/source/source_md/test/test_md_func.cpp index 9e0aef2f6ee..0939a2c204a 100644 --- a/source/source_md/test/test_md_func.cpp +++ b/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); diff --git a/tests/07_OFDFT/25_OF_MD_LDA/INPUT b/tests/07_OFDFT/25_OF_MD_LDA/INPUT index f859f72dffa..89c049d33ba 100644 --- a/tests/07_OFDFT/25_OF_MD_LDA/INPUT +++ b/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 diff --git a/tests/07_OFDFT/26_OF_MD_LibxcPBE/INPUT b/tests/07_OFDFT/26_OF_MD_LibxcPBE/INPUT index 99034663a66..aa24f0a20db 100644 --- a/tests/07_OFDFT/26_OF_MD_LibxcPBE/INPUT +++ b/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 diff --git a/tests/07_OFDFT/30_TDOFDFT_Al/INPUT b/tests/07_OFDFT/30_TDOFDFT_Al/INPUT index 8c5bc4f8b91..569246add7e 100644 --- a/tests/07_OFDFT/30_TDOFDFT_Al/INPUT +++ b/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 diff --git a/tests/07_OFDFT/31_TDOFDFT_Al_CD/INPUT b/tests/07_OFDFT/31_TDOFDFT_Al_CD/INPUT index 0e8a48f6803..885e0f80072 100644 --- a/tests/07_OFDFT/31_TDOFDFT_Al_CD/INPUT +++ b/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 diff --git a/tests/07_OFDFT/32_TDOFDFT_Al_mCD/INPUT b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/INPUT index cfccd2470f3..c6fccd263fd 100644 --- a/tests/07_OFDFT/32_TDOFDFT_Al_mCD/INPUT +++ b/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 diff --git a/tests/integrate/validation_tools/props_basic.sh b/tests/integrate/validation_tools/props_basic.sh index a73902d0978..6980b6d69a4 100644 --- a/tests/integrate/validation_tools/props_basic.sh +++ b/tests/integrate/validation_tools/props_basic.sh @@ -51,8 +51,55 @@ fi #---------------------------- if ! test -z "$has_force" && [ $has_force == 1 ]; then nn3=`echo "$natom + 3" |bc` - # check the last step result - grep -A$nn3 "TOTAL-FORCE" $running_path |awk 'NF==4{print $2,$3,$4}' | tail -$natom > force.txt + # MD atomic forces are stored in the last MD_dump frame. They are no + # longer printed in running_md.log, which is reserved for scalar + # thermodynamic information. Keep the running-log path for non-MD + # calculations and for legacy outputs that have no trajectory dump. + if [ "$calculation" == "md" ] && [ -f OUT.autotest/MD_dump ]; then + awk -v natom="$natom" ' + function reset_frame() { + in_atoms = 0 + atom_count = 0 + last_count = 0 + delete force_x + delete force_y + delete force_z + } + /^MDSTEP:/ { + reset_frame() + next + } + /^INDEX[[:space:]]+LABEL[[:space:]]+POSITION.*FORCE/ { + in_atoms = 1 + atom_count = 0 + delete force_x + delete force_y + delete force_z + next + } + in_atoms && $1 ~ /^[0-9]+$/ && NF >= 8 { + atom_count++ + force_x[atom_count] = $6 + force_y[atom_count] = $7 + force_z[atom_count] = $8 + if (atom_count == natom) { + last_count = atom_count + in_atoms = 0 + } + } + END { + if (last_count != natom) { + exit 1 + } + for (iat = 1; iat <= last_count; ++iat) { + print force_x[iat], force_y[iat], force_z[iat] + } + } + ' OUT.autotest/MD_dump > force.txt + else + # Check the last step result for non-MD calculations. + grep -A$nn3 "TOTAL-FORCE" $running_path | awk 'NF==4{print $2,$3,$4}' | tail -$natom > force.txt + fi total_force=`sum_file force.txt` rm force.txt echo "totalforceref $total_force" >>$props_result_file