vegas/
program.rs

1//! Programs to run on samples.
2//!
3//! A program is a sequence of steps that can be run on a system.
4//! Programs implement the `Program` trait, which requires a `run` method. This method takes a
5//! random number generator and a mutable reference to a `Machine`, allowing the program to
6//! manipulate the machine's state and behavior.
7//!
8//! Some of the provided programs include:
9//!
10//! * `Relax` - A program that relaxes the system at a specified temperature for a given number of steps.
11//! * `CoolDown` - A program that gradually cools down the system from a maximum temperature to a minimum temperature over a series of steps.
12//! * `HysteresisLoop` - A program that simulates a hysteresis loop by varying the external magnetic field and measuring the system's response.
13//!
14//! # Example
15//!
16//! ```rust
17//! use rand::SeedableRng;
18//! use rand_pcg::Pcg64;
19//! use vegas::{
20//!    energy::{Hamiltonian, Zeeman},
21//!    integrator::MetropolisIntegrator,
22//!    machine::Machine,
23//!    program::{CoolDown, Program},
24//!    state::{Field, IsingSpin, State, Spin},
25//!    thermostat::Thermostat,
26//! };
27//!
28//! // Define a Hamiltonian (e.g., Zeeman Energy).
29//! let hamiltonian = Zeeman::new();
30//! let program = CoolDown::default()
31//!    .set_relax(10)
32//!    .set_steps(10);
33//! let mut rng = Pcg64::from_rng(&mut rand::rng());
34//! let state = State::<IsingSpin>::rand_with_size(&mut rng, 100);
35//! let integrator = MetropolisIntegrator::new();
36//! let thermostat = Thermostat::new(2.8, Field::zero());
37//! let instruments = Vec::new();
38//! let mut machine = Machine::new(thermostat, hamiltonian, integrator, instruments, state);
39//! let _ = program.run(&mut rng, &mut machine);
40//! ```
41
42use core::f64;
43
44use crate::{
45    energy::Hamiltonian,
46    error::{ProgramError, ProgramResult},
47    integrator::Integrator,
48    machine::Machine,
49    state::{Field, Spin},
50};
51use rand::Rng;
52use serde::{Deserialize, Serialize};
53
54/// A program is a sequence of steps that can be run on a system.
55pub trait Program {
56    /// Run the program on a system returning the last state.
57    fn run<R, I, H, S>(&self, rng: &mut R, machine: &mut Machine<H, I, S>) -> ProgramResult<()>
58    where
59        R: Rng,
60        I: Integrator<S>,
61        H: Hamiltonian<S>,
62        S: Spin;
63}
64
65/// A program that relaxes the system.
66#[derive(Debug, Deserialize, Serialize)]
67pub struct Relax {
68    steps: usize,
69    temperature: f64,
70}
71
72impl Relax {
73    /// Create a new relaxation program.
74    pub fn new(steps: usize, temperature: f64) -> Self {
75        Self { steps, temperature }
76    }
77
78    /// Set the number of steps.
79    pub fn set_steps(mut self, steps: usize) -> Self {
80        self.steps = steps;
81        self
82    }
83
84    /// Set the temperature.
85    pub fn set_temperature(mut self, temperature: f64) -> Self {
86        self.temperature = temperature;
87        self
88    }
89}
90
91impl Default for Relax {
92    fn default() -> Self {
93        Self::new(1000, 3.0)
94    }
95}
96
97impl Program for Relax {
98    fn run<R, I, H, S>(&self, rng: &mut R, machine: &mut Machine<H, I, S>) -> ProgramResult<()>
99    where
100        I: Integrator<S>,
101        H: Hamiltonian<S>,
102        S: Spin,
103        R: Rng,
104    {
105        if self.steps == 0 {
106            return Err(ProgramError::NoSteps);
107        }
108        if self.temperature < f64::EPSILON {
109            return Err(ProgramError::ZeroTemperature);
110        }
111        machine.set_thermostat(machine.thermostat().with_temperature(self.temperature));
112        machine.relax_for(rng, self.steps)?;
113        Ok(())
114    }
115}
116
117/// A program that cools the system to find the Curie temperature.
118#[derive(Debug, Deserialize, Serialize)]
119pub struct CoolDown {
120    max_temperature: f64,
121    min_temperature: f64,
122    cool_rate: f64,
123    relax: usize,
124    steps: usize,
125}
126
127impl CoolDown {
128    /// Create a new Curie temperature program.
129    pub fn new(
130        max_temperature: f64,
131        min_temperature: f64,
132        cool_rate: f64,
133        relax: usize,
134        steps: usize,
135    ) -> Self {
136        Self {
137            max_temperature,
138            min_temperature,
139            cool_rate,
140            relax,
141            steps,
142        }
143    }
144
145    /// Set the maximum temperature.
146    pub fn set_max_temperature(mut self, max_temp: f64) -> Self {
147        self.max_temperature = max_temp;
148        self
149    }
150
151    /// Set the minimum temperature.
152    pub fn set_min_temperature(mut self, min_temp: f64) -> Self {
153        self.min_temperature = min_temp;
154        self
155    }
156
157    /// Set the cooling rate.
158    pub fn set_cool_rate(mut self, cool_rate: f64) -> Self {
159        self.cool_rate = cool_rate;
160        self
161    }
162
163    /// Set the number of relaxation steps.
164    pub fn set_relax(mut self, relax: usize) -> Self {
165        self.relax = relax;
166        self
167    }
168
169    /// Set the number of steps.
170    pub fn set_steps(mut self, steps: usize) -> Self {
171        self.steps = steps;
172        self
173    }
174}
175
176impl Default for CoolDown {
177    fn default() -> Self {
178        Self::new(3.0, 0.05, 0.1, 1000, 20000)
179    }
180}
181
182impl Program for CoolDown {
183    fn run<R, I, H, S>(&self, rng: &mut R, machine: &mut Machine<H, I, S>) -> ProgramResult<()>
184    where
185        I: Integrator<S>,
186        H: Hamiltonian<S>,
187        S: Spin,
188        R: Rng,
189    {
190        if self.max_temperature < self.min_temperature {
191            return Err(ProgramError::TemperatureMaxLessThanMin);
192        }
193        if self.steps == 0 {
194            return Err(ProgramError::NoSteps);
195        }
196        if self.min_temperature < f64::EPSILON {
197            return Err(ProgramError::ZeroTemperature);
198        }
199        if self.cool_rate < f64::EPSILON {
200            return Err(ProgramError::ZeroCoolRate);
201        }
202        let mut temperature = self.max_temperature;
203        loop {
204            machine.set_thermostat(machine.thermostat().with_temperature(temperature));
205            machine.relax_for(rng, self.relax)?;
206            machine.measure_for(rng, self.steps)?;
207            temperature -= self.cool_rate;
208            if temperature < self.min_temperature {
209                break;
210            }
211        }
212        Ok(())
213    }
214}
215
216/// A program that runs a hysteresis loop.
217#[derive(Debug, Deserialize, Serialize)]
218pub struct HysteresisLoop {
219    steps: usize,
220    relax: usize,
221    temperature: f64,
222    max_field: f64,
223    field_step: f64,
224}
225
226impl HysteresisLoop {
227    /// Create a new hysteresis loop program.
228    pub fn new(
229        steps: usize,
230        relax: usize,
231        temperature: f64,
232        max_field: f64,
233        field_step: f64,
234    ) -> Self {
235        Self {
236            steps,
237            relax,
238            temperature,
239            max_field,
240            field_step,
241        }
242    }
243
244    /// Set the number of steps.
245    pub fn set_steps(mut self, steps: usize) -> Self {
246        self.steps = steps;
247        self
248    }
249
250    /// Set the number of relaxation steps.
251    pub fn set_relax(mut self, relax: usize) -> Self {
252        self.relax = relax;
253        self
254    }
255
256    /// Set the temperature.
257    pub fn set_temperature(mut self, temperature: f64) -> Self {
258        self.temperature = temperature;
259        self
260    }
261
262    /// Set the maximum field.
263    pub fn set_max_field(mut self, max_field: f64) -> Self {
264        self.max_field = max_field;
265        self
266    }
267
268    /// Set the field step.
269    pub fn set_field_step(mut self, field_step: f64) -> Self {
270        self.field_step = field_step;
271        self
272    }
273}
274
275impl Default for HysteresisLoop {
276    fn default() -> Self {
277        Self::new(1000, 1000, 3.0, 1.0, 0.1)
278    }
279}
280
281impl Program for HysteresisLoop {
282    fn run<R, I, H, S>(&self, rng: &mut R, machine: &mut Machine<H, I, S>) -> ProgramResult<()>
283    where
284        R: Rng,
285        I: Integrator<S>,
286        H: Hamiltonian<S>,
287        S: Spin,
288    {
289        if self.steps == 0 {
290            return Err(ProgramError::NoSteps);
291        }
292        if self.temperature < f64::EPSILON {
293            return Err(ProgramError::ZeroTemperature);
294        }
295        if self.max_field < f64::EPSILON {
296            return Err(ProgramError::ZeroField);
297        }
298        if self.field_step < f64::EPSILON {
299            return Err(ProgramError::ZeroFieldStep);
300        }
301        machine.set_thermostat(machine.thermostat().with_temperature(self.temperature));
302        let mut magnitude = 0.0;
303        loop {
304            let field = Field::new(S::up(), magnitude);
305            machine.set_thermostat(machine.thermostat().with_field(field));
306            machine.relax_for(rng, self.relax)?;
307            machine.measure_for(rng, self.steps)?;
308            magnitude += self.field_step;
309            if magnitude > self.max_field {
310                break;
311            }
312        }
313        loop {
314            let field = Field::new(S::up(), magnitude);
315            machine.set_thermostat(machine.thermostat().with_field(field));
316            machine.relax_for(rng, self.relax)?;
317            machine.measure_for(rng, self.steps)?;
318            magnitude -= self.field_step;
319            if magnitude < -self.max_field {
320                break;
321            }
322        }
323        loop {
324            let field = Field::new(S::up(), magnitude);
325            machine.set_thermostat(machine.thermostat().with_field(field));
326            machine.relax_for(rng, self.relax)?;
327            machine.measure_for(rng, self.steps)?;
328            magnitude += self.field_step;
329            if magnitude > self.max_field {
330                break;
331            }
332        }
333
334        Ok(())
335    }
336}
vegas/program.rs

vegas/
program.rs
