hoomd_md/thermostat/
bussi.rs

1// Copyright (c) 2024-2026 The Regents of the University of Michigan.
2// Part of hoomd-rs, released under the BSD 3-Clause License.
3
4//! Implement `Bussi`
5
6use rand::Rng;
7use rand_distr::{Distribution, Gamma, Normal};
8use serde::{Deserialize, Serialize};
9
10use crate::Thermostat;
11use hoomd_simulation::macrostate::Temperature;
12
13/// Stochastic momentum rescaling.
14///
15/// The time constant $` \tau `$ sets how long the kinetic energy of the system
16/// takes to decay to equilibrium.
17///
18/// When $`\tau`$ is 0, the rescaling factor $` \alpha `$ is:
19/// ```math
20///  \alpha = \sqrt{\frac{g kT}{K}}
21/// ```
22/// where $`K`$ is the instantaneous kinetic energy of the corresponding
23/// translational or rotational degrees of freedom, $`N`$ is the number of
24/// degrees of freedom, and $`g`$ is a random value sampled from the gamma
25/// distribution $`\Gamma(N, 1)`$ with the probability density function:
26/// ```math
27///    f_N(g) = \frac{1}{\Gamma{(N)}} g^{N-1} e^{-g}
28/// ```
29///
30/// When $`\tau`$ is non-zero, $` \alpha `$ is given by:
31/// ```math
32/// \alpha = \sqrt{e^{-\delta t / \tau}
33///      + (1 - e^{-\delta t / \tau}) \frac{(2 g + n^2) kT}{2 K}
34///      + 2 n \sqrt{e^{-\delta t / \tau} (1-e^{-\delta t / \tau})
35///         \frac{kT}{2 K}}}
36/// ```
37/// where $`\delta t`$ is the integration time step size and $`n`$ is a random
38/// value sampled from the standard normal distribution $`\mathcal{N}(0, 1)`$.
39///
40/// # Reference
41///
42/// * [Bussi et al. 2007](https://doi.org/10.1063/1.2408420)
43///
44/// # Example
45///
46/// ```
47/// use hoomd_md::thermostat::Bussi;
48///
49/// let bussi = Bussi::new(0.5);
50/// ```
51#[derive(Clone, Debug, Default, PartialEq, Serialize, Deserialize)]
52pub struct Bussi {
53    /// Thermostat time constant $`[ \mathrm{time} ]`$.
54    tau: f64,
55
56    /// Cumulative energy absorbed by the thermostat $`[ \mathrm{energy} ]`$
57    energy: f64,
58}
59
60impl Bussi {
61    /// Construct a new `Bussi` thermostat with the given time constant.
62    ///
63    /// # Example
64    ///
65    /// ```
66    /// use hoomd_md::thermostat::Bussi;
67    ///
68    /// let bussi = Bussi::new(0.5);
69    /// ```
70    #[inline]
71    pub fn new(tau: f64) -> Self {
72        Self { tau, energy: 0.0 }
73    }
74    /// Calculate the energy drift on one time step.
75    #[inline]
76    fn energy_drift(kinetic_energy_old: f64, rescaling_factor: f64) -> f64 {
77        kinetic_energy_old * (1.0 - rescaling_factor.powi(2))
78    }
79
80    /// The total energy of the thermostat.
81    /// # Example
82    ///
83    /// ```
84    /// use hoomd_md::thermostat::Bussi;
85    ///
86    /// let bussi = Bussi::new(0.1);
87    /// let energy = bussi.energy();
88    /// ```
89    #[inline]
90    pub fn energy(&self) -> f64 {
91        self.energy
92    }
93}
94
95impl<M> Thermostat<M> for Bussi
96where
97    M: Temperature,
98{
99    #[inline]
100    fn integrate_half_step_one<R: Rng + ?Sized>(
101        &mut self,
102        rng: &mut R,
103        macrostate: &M,
104        delta_t: f64,
105        kinetic_energy: f64,
106        degrees_of_freedom: usize,
107    ) -> f64 {
108        // Calculate velocity rescaling factor following
109        // the Appendix in https://doi.org/10.1063/1.2408420.
110        if degrees_of_freedom == 0 {
111            return 1.0;
112        }
113
114        assert!(
115            kinetic_energy != 0.0,
116            "The Bussi thermostat requires non-zero kinetic energy."
117        );
118
119        let time_decay_factor = if self.tau == 0.0 {
120            0.0
121        } else {
122            (-delta_t / self.tau).exp()
123        };
124
125        let random_normal_one: f64 = Normal::new(0.0, 1.0)
126            .expect("normal distribution should be valid")
127            .sample(rng);
128
129        let random_gamma = if degrees_of_freedom > 0 {
130            2.0 * Gamma::new((degrees_of_freedom as f64 - 1.0) / 2.0, 1.0)
131                .expect("gamma distribution should be valid")
132                .sample(rng)
133        } else {
134            0.0
135        };
136
137        let v = macrostate.temperature() / 2.0 / kinetic_energy;
138        let term1 = v * (1.0 - time_decay_factor) * (random_gamma + random_normal_one.powi(2));
139        let term2 =
140            2.0 * random_normal_one * (v * (1.0 - time_decay_factor) * time_decay_factor).sqrt();
141        let rescaling_factor = (time_decay_factor + term1 + term2).sqrt();
142
143        self.energy += Self::energy_drift(kinetic_energy, rescaling_factor);
144        rescaling_factor
145    }
146
147    #[inline]
148    fn integrate_half_step_two<R: Rng + ?Sized>(
149        &mut self,
150        _rng: &mut R,
151        _macrostate: &M,
152        _delta_t: f64,
153        _kinetic_energy: f64,
154        _degrees_of_freedom: usize,
155    ) -> f64 {
156        1.0
157    }
158}
159
160#[cfg(test)]
161mod tests {
162    use super::*;
163    use assert2::check;
164    use rand::{SeedableRng, rngs::StdRng};
165
166    use hoomd_simulation::macrostate::Isothermal;
167
168    #[test]
169    fn test_init() {
170        let tau = 1.0;
171        let bussi = Bussi::new(tau);
172
173        check!(bussi.tau == tau);
174        check!(bussi.energy() == 0.0);
175    }
176
177    #[test]
178    fn test_scale_down() {
179        let tau = 0.0;
180        let mut bussi = Bussi::new(tau);
181
182        let mut rng = StdRng::seed_from_u64(0);
183        let macrostate = Isothermal { temperature: 0.4 };
184        let alpha = bussi.integrate_half_step_one(&mut rng, &macrostate, 0.01, 1000.0, 100);
185
186        check!(alpha < 1.0);
187        check!(bussi.energy() > 0.0);
188    }
189
190    #[test]
191    fn test_scale_up() {
192        let tau = 0.0;
193        let mut bussi = Bussi::new(tau);
194
195        let mut rng = StdRng::seed_from_u64(0);
196        let macrostate = Isothermal { temperature: 0.4 };
197        let alpha = bussi.integrate_half_step_one(&mut rng, &macrostate, 0.01, 1000.0, 10_000);
198
199        check!(alpha > 1.0);
200        check!(bussi.energy() < 0.0);
201    }
202}
hoomd_md/thermostat/bussi.rs

hoomd_md/thermostat/
bussi.rs
