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

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

@@ -3,6 +3,8 @@
 #include "md_func.h"
 #include "source_base/timer.h"
 
+#include <random>
+
 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<double> 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);

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