feat(md): add CSVR thermostat for NVT molecular dynamics (#7461)

* 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

* docs(md): add CSVR thermostat documentation

- Add csvr option to md_thermostat parameter description
- Add md_csvr_tau parameter documentation

Implements #6941

* 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

* fix: address review comments for CSVR thermostat

- Fix degrees of freedom: use 3N - frozen_freedom instead of frozen_freedom
- Use MD_func::gaussrand() instead of std::random for consistency
- Simplify CSVR core formula with factor2 variable
- Remove unused #include <random>
Co-authored-by: monkeycode-ai <monkeycode-ai@chaitin.com>
Co-authored-by: monkeycode-ai <monkeycode-ai@chaitin.com>

---------

Co-authored-by: monkeycode-ai <monkeycode-ai@chaitin.com>

Author: mintleaf84

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

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
 

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)";

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;

source/source_md/verlet.cpp

@@ -100,6 +100,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 +131,62 @@ 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 (3N - frozen)
+    int ndeg = 3 * ucell.nat - 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 numbers using MD_func
+    double rr = MD_func::gaussrand();
+
+    // Calculate sum of squared Gaussian random numbers (ndeg - 1)
+    double sumnoises = 0.0;
+    for (int i = 0; i < ndeg - 1; ++i)
+    {
+        double r = MD_func::gaussrand();
+        sumnoises += r * r;
+    }
+
+    // CSVR core formula (simplified)
+    double factor2 = (1.0 - factor) * kin_target / kin_energy / ndeg;
+    double resample = factor + factor2 * (rr * rr + sumnoises) + 2.0 * rr * sqrt(factor * factor2);
+
+    // Ensure non-negative
+    resample = std::max(0.0, resample);
+
+    // Calculate scaling factor
+    double scale = sqrt(resample);
+
+    // 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);

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