hoomd_md/thermostat/nose_hoover_chain.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 `NoséHooverChain`
5
6use rand::Rng;
7use rand_distr::{Distribution, Normal};
8use serde::{Deserialize, Serialize};
9use serde_with::serde_as;
10
11use crate::Thermostat;
12use hoomd_simulation::macrostate::Temperature;
13use hoomd_utility::valid::PositiveReal;
14
15/// Chain of Nosé-Hoover thermostats.
16///
17/// [`NoséHooverChain`] adds new degrees of freedom ($`\eta_i`$)
18/// to a molecular dynamics simulation in such a way that the existing
19/// degrees of freedom sample a constant temperature ensemble. Each
20/// [`NoséHooverChain`] instance stores the $`\eta_i`$ and their momenta,
21/// $`\xi_i`$, internally.
22///
23/// The dynamics of each $`\eta_i, \xi_i`$ are similar to that in
24/// [`MartynaTuckermanTobiasKlein`], but they are also chained together.
25/// See [Martyna et al. 1992] for details.
26///
27/// [`MartynaTuckermanTobiasKlein`]: crate::thermostat::MartynaTuckermanTobiasKlein
28///
29/// # Reference
30///
31/// * [Martyna et al. 1992]
32///
33/// [Martyna et al. 1992]: https://doi.org/10.1063/1.463940
34///
35/// # Example
36///
37/// ```
38/// use hoomd_md::thermostat::NoséHooverChain;
39///
40/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
41/// let thermostat = NoséHooverChain::<3>::zero(0.5.try_into()?);
42/// # Ok(())
43/// # }
44/// ```
45#[serde_as]
46#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
47pub struct NoséHooverChain<const N: usize> {
48 /// Thermostat time constant.
49 tau: PositiveReal,
50
51 /// Chain of thermostat momenta.
52 #[serde_as(as = "[_; N]")]
53 xi: [f64; N],
54
55 /// Chain of thermostat positions.
56 #[serde_as(as = "[_; N]")]
57 eta: [f64; N],
58
59 /// Chain of thermostat accelerations.
60 #[serde_as(as = "[_; N]")]
61 g: [f64; N],
62
63 /// Energy the thermostat contributes to the Hamiltonian.
64 energy: f64,
65}
66
67impl<const N: usize> NoséHooverChain<N> {
68 /// Construct a new `NoséHooverChain` thermostat with the given time constant,
69 /// $` \xi_i = 0 `$, and $` \eta_i = 0 `$ .
70 ///
71 /// This initial condition is likely to be very far from equilibrium which
72 /// will result in wild kinetic energy oscillations for the first hundred to
73 /// thousand time steps. Use [`thermalized`] to choose the initial position
74 /// and momentum from a thermal distribution.
75 ///
76 /// [`thermalized`]: Self::thermalized
77 ///
78 /// # Example
79 ///
80 /// ```
81 /// use hoomd_md::thermostat::NoséHooverChain;
82 ///
83 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
84 /// let thermostat = NoséHooverChain::<3>::zero(0.5.try_into()?);
85 /// # Ok(())
86 /// # }
87 /// ```
88 #[inline]
89 pub fn zero(tau: PositiveReal) -> Self {
90 Self {
91 tau,
92 xi: [0.0; N],
93 eta: [0.0; N],
94 g: [0.0; N],
95 energy: 0.0,
96 }
97 }
98
99 /// Construct a new `NoséHooverChain` thermostat with random $` \xi_i `$
100 /// values drawn from a thermal distribution.
101 ///
102 /// # Panics
103 ///
104 /// This method will panic when `degrees_of_freedom` is 0.
105 ///
106 /// # Example
107 ///
108 /// ```
109 /// use hoomd_md::{TranslationalKineticEnergy, thermostat::NoséHooverChain};
110 /// use hoomd_microstate::{
111 /// Body, Microstate,
112 /// property::{DynamicPoint, Point},
113 /// };
114 /// use hoomd_simulation::macrostate::Isothermal;
115 /// use hoomd_vector::Cartesian;
116 ///
117 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
118 /// let mut microstate = Microstate::builder()
119 /// .bodies([
120 /// Body::single_site(
121 /// DynamicPoint {
122 /// position: Cartesian::from([1.0, 2.0]),
123 /// ..Default::default()
124 /// },
125 /// Point::default(),
126 /// ),
127 /// Body::single_site(
128 /// DynamicPoint {
129 /// position: Cartesian::from([-2.0, 3.0]),
130 /// ..Default::default()
131 /// },
132 /// Point::default(),
133 /// ),
134 /// ])
135 /// .try_build()?;
136 ///
137 /// let macrostate = Isothermal { temperature: 1.5 };
138 /// let mut rng = microstate.counter().make_rng();
139 /// let translational_thermostat = NoséHooverChain::<3>::thermalized(
140 /// &mut rng,
141 /// 0.5.try_into()?,
142 /// ¯ostate,
143 /// microstate.translational_kinetic_energy().1,
144 /// );
145 /// microstate.increment_substep();
146 /// # Ok(())
147 /// # }
148 /// ```
149 #[inline]
150 pub fn thermalized<M, R: Rng + ?Sized>(
151 rng: &mut R,
152 tau: PositiveReal,
153 macrostate: &M,
154 degrees_of_freedom: usize,
155 ) -> Self
156 where
157 M: Temperature,
158 {
159 let sigma_0 = 1.0 / (degrees_of_freedom as f64).sqrt() / tau.get();
160 let sigma_other = 1.0 / tau.get();
161
162 let mut xi = [0.0; N];
163
164 xi[0] = Normal::new(0.0, sigma_0)
165 .expect("Normal distribution should be valid")
166 .sample(rng);
167
168 for xi_i in xi.iter_mut().skip(1) {
169 *xi_i = Normal::new(0.0, sigma_other)
170 .expect("Normal distribution should be valid")
171 .sample(rng);
172 }
173
174 let mut result = Self {
175 tau,
176 xi,
177 eta: [0.0; N],
178 g: [0.0; N],
179 energy: 0.0,
180 };
181
182 let q = result.q(*macrostate.temperature(), degrees_of_freedom);
183 result.energy = result.thermostat_energy(*macrostate.temperature(), degrees_of_freedom, &q);
184
185 result
186 }
187
188 /// Calculate q.
189 #[inline]
190 fn q(&self, temperature: f64, degrees_of_freedom: usize) -> [f64; N] {
191 let n_k_t = (degrees_of_freedom as f64) * temperature;
192 let mut result = [temperature * self.tau.get().powi(2); N];
193
194 result[0] = n_k_t * self.tau.get().powi(2);
195
196 result
197 }
198
199 /// Calculate thermostat energy.
200 #[inline]
201 fn thermostat_energy(
202 &self,
203 temperature_set_point: f64,
204 degrees_of_freedom: usize,
205 q: &[f64; N],
206 ) -> f64 {
207 let mut energy = 0.0;
208 energy += (degrees_of_freedom as f64) * temperature_set_point * self.eta[0]
209 + 0.5 * q[0] * (self.xi[0]).powi(2);
210
211 for (eta_i, (q_i, xi_i)) in self.eta.iter().zip(q.iter().zip(self.xi)) {
212 energy += temperature_set_point * eta_i + 0.5 * q_i * (xi_i).powi(2);
213 }
214 energy
215 }
216
217 /// The total energy of the thermostat.
218 ///
219 /// # Example
220 ///
221 /// ```
222 /// use hoomd_md::thermostat::NoséHooverChain;
223 ///
224 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
225 /// let thermostat = NoséHooverChain::<3>::zero(0.5.try_into()?);
226 ///
227 /// let energy = thermostat.energy();
228 /// # Ok(())
229 /// # }
230 /// ```
231 #[inline]
232 pub fn energy(&self) -> f64 {
233 self.energy
234 }
235
236 /// The chain of thermostat positions.
237 ///
238 /// # Example
239 ///
240 /// ```
241 /// use hoomd_md::thermostat::NoséHooverChain;
242 ///
243 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
244 /// let thermostat = NoséHooverChain::<3>::zero(0.5.try_into()?);
245 ///
246 /// let eta = thermostat.eta();
247 /// # Ok(())
248 /// # }
249 /// ```
250 #[inline]
251 pub fn eta(&self) -> &[f64; N] {
252 &self.eta
253 }
254
255 /// The chain of thermostat momenta.
256 ///
257 /// # Example
258 ///
259 /// ```
260 /// use hoomd_md::thermostat::NoséHooverChain;
261 ///
262 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
263 /// let thermostat = NoséHooverChain::<3>::zero(0.5.try_into()?);
264 ///
265 /// let xi = thermostat.xi();
266 /// # Ok(())
267 /// # }
268 /// ```
269 #[inline]
270 pub fn xi(&self) -> &[f64; N] {
271 &self.xi
272 }
273}
274
275impl<const N: usize, M> Thermostat<M> for NoséHooverChain<N>
276where
277 M: Temperature,
278{
279 #[inline]
280 fn integrate_half_step_one<R: Rng + ?Sized>(
281 &mut self,
282 _rng: &mut R,
283 macrostate: &M,
284 delta_t: f64,
285 kinetic_energy: f64,
286 degrees_of_freedom: usize,
287 ) -> f64 {
288 // Integrate extra degrees-of-freedom and
289 // return the velocity rescaling factor, following
290 // Tuckerman's work <https://doi.org/10.1088/0305-4470/39/19/S18>.
291
292 let n_k_t = (degrees_of_freedom as f64) * *macrostate.temperature();
293 let q = self.q(*macrostate.temperature(), degrees_of_freedom);
294
295 // Update the thermostat acceleration coupled to the real system
296 self.g[0] = (2.0 * kinetic_energy - n_k_t) / q[0];
297
298 // Update the chain of velocity
299 // start from the last one
300 self.xi[N - 1] += 0.25 * delta_t * self.g[N - 1];
301 // update the rest
302 for idx in (0..N - 1).rev() {
303 let xi_rescaling_factor = (-0.125 * delta_t * self.xi[idx + 1]).exp();
304 self.xi[idx] *= xi_rescaling_factor;
305 self.xi[idx] += 0.25 * delta_t * self.g[idx];
306 self.xi[idx] *= xi_rescaling_factor;
307 }
308
309 // calculate real velocity rescaling factor
310 let rescaling_factor = (-0.5 * delta_t * self.xi[0]).exp();
311
312 // Expected temperature update
313 let kinetic_energy_new = kinetic_energy * rescaling_factor.powi(2);
314
315 // Update the thermostat acceleration coupled to the real system
316 self.g[0] = (2.0 * kinetic_energy_new - n_k_t) / q[0];
317
318 // Update the chain of position
319 for idx in 0..N {
320 self.eta[idx] += 0.5 * delta_t * self.xi[idx];
321 }
322
323 // Update the chain of velocity
324 // start from the first one
325 if N > 1 {
326 let xi_rescaling_factor = (-0.125 * delta_t * self.xi[1]).exp();
327 self.xi[0] *= xi_rescaling_factor;
328 self.xi[0] += 0.25 * delta_t * self.g[0];
329 self.xi[0] *= xi_rescaling_factor;
330 } else {
331 self.xi[0] += 0.25 * delta_t * self.g[0];
332 }
333 // update the rest
334 // the chain of acceleration need to be updated here (have done the first one)
335 for idx in 1..N - 1 {
336 let xi_rescaling_factor = (-0.125 * delta_t * self.xi[idx + 1]).exp();
337 self.xi[idx] *= xi_rescaling_factor;
338 self.g[idx] =
339 (q[idx - 1] * (self.xi[idx - 1]).powi(2) - *macrostate.temperature()) / q[idx];
340 self.xi[idx] += 0.25 * delta_t * self.g[idx];
341 self.xi[idx] *= xi_rescaling_factor;
342 }
343 // special for the last one
344 if N > 1 {
345 self.g[N - 1] =
346 (q[N - 2] * (self.xi[N - 2]).powi(2) - *macrostate.temperature()) / q[N - 1];
347 self.xi[N - 1] += 0.25 * delta_t * self.g[N - 1];
348 }
349
350 self.energy = self.thermostat_energy(*macrostate.temperature(), degrees_of_freedom, &q);
351 rescaling_factor
352 }
353
354 #[inline]
355 fn integrate_half_step_two<R: Rng + ?Sized>(
356 &mut self,
357 rng: &mut R,
358 macrostate: &M,
359 delta_t: f64,
360 kinetic_energy: f64,
361 degrees_of_freedom: usize,
362 ) -> f64 {
363 self.integrate_half_step_one(rng, macrostate, delta_t, kinetic_energy, degrees_of_freedom)
364 }
365}
366
367#[cfg(test)]
368mod tests {
369 use super::*;
370 use assert2::check;
371
372 use crate::TranslationalKineticEnergy;
373 use hoomd_microstate::{
374 Body, Microstate,
375 property::{DynamicPoint, Point},
376 };
377 use hoomd_simulation::macrostate::Isothermal;
378 use hoomd_vector::Cartesian;
379
380 #[test]
381 fn test_zero() -> anyhow::Result<()> {
382 let thermostat = NoséHooverChain::<10>::zero(0.5.try_into()?);
383
384 check!(thermostat.tau.get() == 0.5);
385 check!(thermostat.xi() == &[0.0; 10]);
386 check!(thermostat.eta() == &[0.0; 10]);
387 check!(thermostat.energy() == 0.0);
388
389 Ok(())
390 }
391
392 #[test]
393 fn test_thermalized() -> anyhow::Result<()> {
394 let microstate = Microstate::builder()
395 .bodies([
396 Body {
397 properties: DynamicPoint {
398 position: Cartesian::from([1.0, 2.0]),
399 ..Default::default()
400 },
401 sites: vec![Point::default()],
402 },
403 Body {
404 properties: DynamicPoint {
405 position: Cartesian::from([-2.0, 3.0]),
406 ..Default::default()
407 },
408 sites: vec![Point::default()],
409 },
410 ])
411 .try_build()?;
412
413 let macrostate = Isothermal { temperature: 1.5 };
414 let mut rng = microstate.counter().make_rng();
415 let thermostat = NoséHooverChain::<10>::thermalized(
416 &mut rng,
417 0.5.try_into()?,
418 ¯ostate,
419 microstate.translational_kinetic_energy().1,
420 );
421
422 check!(thermostat.tau.get() == 0.5);
423 check!(thermostat.xi() != &[0.0; 10]);
424 check!(thermostat.eta() == &[0.0; 10]);
425 check!(thermostat.energy() != 0.0);
426
427 Ok(())
428 }
429}