From 58e888e6ecbafa7cfff143906d62b6d788071904 Mon Sep 17 00:00:00 2001 From: jiebin chen Date: Mon, 8 Jun 2026 17:28:02 +0000 Subject: [PATCH 1/3] feat(md): add CSVR thermostat for NVT molecular dynamics Implement the Canonical Sampling through Velocity Rescaling (CSVR) thermostat as described in: G. Bussi, D. Donadio, M. Parrinello, J. Chem. Phys. 126, 014101 (2007) Features: - New thermostat option: md_thermostat = csvr - New parameter: md_csvr_tau (characteristic time scale) - Properly samples the canonical (NVT) ensemble - Simple implementation with only one parameter Implements #6941 --- .../source_io/module_parameter/md_parameter.h | 2 + .../module_parameter/read_input_item_md.cpp | 14 ++++ source/source_md/verlet.cpp | 71 +++++++++++++++++++ source/source_md/verlet.h | 8 +++ 4 files changed, 95 insertions(+) diff --git a/source/source_io/module_parameter/md_parameter.h b/source/source_io/module_parameter/md_parameter.h index 80d04b0e94d..a84b23a3aaf 100644 --- a/source/source_io/module_parameter/md_parameter.h +++ b/source/source_io/module_parameter/md_parameter.h @@ -58,6 +58,8 @@ struct MD_para double md_damp = 1.0; ///< Langevin damping parameter (time units) + double md_csvr_tau = 100.0; ///< CSVR thermostat characteristic time scale (in MD time units) + double md_tolerance = 100.0; ///< tolerance for velocity rescaling (K) int md_nraise = 1; ///< parameters used when md_type=nvt 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 30b6f2f6a53..76d393b438b 100644 --- a/source/source_io/module_parameter/read_input_item_md.cpp +++ b/source/source_io/module_parameter/read_input_item_md.cpp @@ -661,6 +661,20 @@ Note: It is a system-dependent empirical parameter. An improper choice might lea read_sync_double(input.mdp.md_damp); this->add_item(item); } + { + Input_Item item("md_csvr_tau"); + item.annotation = "CSVR thermostat characteristic time scale"; + item.category = "Molecular dynamics"; + item.type = "Real"; + item.description = "The characteristic time scale for the CSVR (Canonical Sampling through Velocity " + "Rescaling) thermostat. Larger values give weaker coupling, smaller values give " + "stronger coupling. Recommended value: 100 * md_dt."; + item.default_value = "100.0"; + item.unit = "fs"; + item.availability = "md_thermostat = csvr"; + read_sync_double(input.mdp.md_csvr_tau); + this->add_item(item); + } { Input_Item item("md_tolerance"); item.annotation = "tolerance for velocity rescaling (K)"; diff --git a/source/source_md/verlet.cpp b/source/source_md/verlet.cpp index b2df7a78ced..7189a989321 100644 --- a/source/source_md/verlet.cpp +++ b/source/source_md/verlet.cpp @@ -3,6 +3,8 @@ #include "md_func.h" #include "source_base/timer.h" +#include + Verlet::Verlet(const Parameter& param_in, UnitCell& unit_in) : MD_base(param_in, unit_in) { } @@ -100,6 +102,11 @@ void Verlet::apply_thermostat(void) t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, md_tfirst, md_tlast); thermalize(mdp.md_nraise, t_current, t_target); } + else if (mdp.md_thermostat == "csvr") + { + t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_nstep, md_tfirst, md_tlast); + apply_csvr(t_current, t_target); + } else { ModuleBase::WARNING_QUIT("Verlet", "No such thermostat!"); @@ -126,6 +133,70 @@ void Verlet::thermalize(const int& nraise, const double& current_temp, const dou } +void Verlet::apply_csvr(const double& current_temp, const double& target_temp) +{ + // CSVR thermostat: Canonical Sampling through Velocity Rescaling + // Reference: G. Bussi, D. Donadio, M. Parrinello, J. Chem. Phys. 126, 014101 (2007) + + if (current_temp <= 0.0 || target_temp <= 0.0) + { + return; + } + + // Get degrees of freedom + int ndeg = frozen_freedom_; + + // Calculate kinetic energies + double kin_energy = current_temp * ndeg * 0.5; // in Hartree + double kin_target = target_temp * ndeg * 0.5; // in Hartree + + // Calculate tau parameter (characteristic time scale / dt) + double taut = mdp.md_csvr_tau / mdp.md_dt; + + // Calculate decay factor + double factor = 0.0; + if (taut > 0.1) + { + factor = exp(-1.0 / taut); + } + + // Generate Gaussian random number + std::random_device rd; + std::mt19937 gen(rd()); + if (mdp.md_seed > 0) + { + gen.seed(mdp.md_seed + step_); + } + std::normal_distribution gauss_dist(0.0, 1.0); + double rr = gauss_dist(gen); + + // Calculate sum of squared Gaussian random numbers (ndeg - 1) + double sumnoises = 0.0; + for (int i = 0; i < ndeg - 1; ++i) + { + double r = gauss_dist(gen); + sumnoises += r * r; + } + + // CSVR core formula + double resample = kin_energy + + (1.0 - factor) * (kin_target * (sumnoises + rr * rr) / ndeg - kin_energy) + + 2.0 * rr * sqrt(kin_energy * kin_target / ndeg * (1.0 - factor) * factor); + + // Ensure non-negative + resample = std::max(0.0, resample); + + // Calculate scaling factor + double scale = sqrt(resample / kin_energy); + + // Apply velocity scaling + for (int i = 0; i < ucell.nat; ++i) + { + vel[i] *= scale; + } +} + + void Verlet::print_md(std::ofstream& ofs, const bool& cal_stress) { MD_base::print_md(ofs, cal_stress); diff --git a/source/source_md/verlet.h b/source/source_md/verlet.h index 72e0edcd6e3..9c4ab59d8aa 100644 --- a/source/source_md/verlet.h +++ b/source/source_md/verlet.h @@ -35,6 +35,14 @@ class Verlet : public MD_base * @param target_temp the target temperature */ void thermalize(const int& nraise, const double& current_temp, const double& target_temp); + + /** + * @brief apply CSVR thermostat + * + * @param current_temp the current temperature + * @param target_temp the target temperature + */ + void apply_csvr(const double& current_temp, const double& target_temp); }; #endif \ No newline at end of file From 3c89d41d676992d9d9025b55a26565994cb73094 Mon Sep 17 00:00:00 2001 From: jiebin chen Date: Mon, 8 Jun 2026 17:31:48 +0000 Subject: [PATCH 2/3] docs(md): add CSVR thermostat documentation - Add csvr option to md_thermostat parameter description - Add md_csvr_tau parameter documentation Implements #6941 --- docs/advanced/input_files/input-main.md | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index a91c6dd7530..21088c6c6dd 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -3194,6 +3194,7 @@ - berendsen: Berendsen thermostat, see md_nraise in detail. - rescaling: velocity Rescaling method 1, see md_tolerance in detail. - rescale_v: velocity Rescaling method 2, see md_nraise in detail. + - csvr: Canonical Sampling through Velocity Rescaling, see md_csvr_tau in detail. - **Default**: nhc ### md_tfirst @@ -3445,6 +3446,13 @@ - **Default**: 1.0 - **Unit**: fs +### md_csvr_tau + +- **Type**: Real +- **Description**: The characteristic time scale for the CSVR (Canonical Sampling through Velocity Rescaling) thermostat. Larger values give weaker coupling (longer relaxation time), smaller values give stronger coupling (shorter relaxation time). Recommended value: 100 * md_dt. +- **Default**: 100.0 +- **Unit**: fs + ### md_tolerance - **Type**: Real From 550f7adc1720c16cea183a7b1ffcb4780bdfbf36 Mon Sep 17 00:00:00 2001 From: jiebin chen Date: Mon, 8 Jun 2026 17:40:25 +0000 Subject: [PATCH 3/3] test(md): add unit test for CSVR thermostat Add CSVR thermostat test case to verlet_test.cpp: - Test position update correctness - Verify temperature is in reasonable range Implements #6941 --- source/source_md/test/verlet_test.cpp | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/source/source_md/test/verlet_test.cpp b/source/source_md/test/verlet_test.cpp index 8f2c00f74f3..86dc66a85e1 100644 --- a/source/source_md/test/verlet_test.cpp +++ b/source/source_md/test/verlet_test.cpp @@ -281,6 +281,26 @@ TEST_F(Verlet_test, rescale_v) EXPECT_NEAR(mdrun->vel[3].z, -2.8328663233253657e-05, doublethreshold); } +TEST_F(Verlet_test, CSVR) +{ + mdrun->first_half(GlobalV::ofs_running); + param_in.input.mdp.md_type = "nvt"; + param_in.input.mdp.md_thermostat = "csvr"; + param_in.input.mdp.md_csvr_tau = 100.0; + param_in.input.mdp.md_seed = 12345; + mdrun->second_half(); + + // Check that positions are updated correctly + EXPECT_NEAR(mdrun->pos[0].x, -0.00054545529007222658, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].y, 0.00029590658162135359, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].z, -5.7952328034033513e-05, doublethreshold); + + // Check that temperature is in reasonable range + double temp = mdrun->t_current * ModuleBase::Hartree_to_K; + EXPECT_GT(temp, 0.0); + EXPECT_LT(temp, 1000.0); +} + TEST_F(Verlet_test, write_restart) { mdrun->step_ = 1;