lattice_qcd_rs/simulation/monte_carlo/
heat_bath.rs

1//! Pseudo heat bath methods
2//!
3//! # Example
4//! see [`HeatBathSweep`].
5
6use na::ComplexField;
7#[cfg(feature = "serde-serialize")]
8use serde::{Deserialize, Serialize};
9
10use super::{
11    super::{
12        super::{
13            error::Never, lattice::LatticeLinkCanonical, statistics::HeatBathDistribution, su2,
14            su3, CMatrix2, Complex,
15        },
16        state::{LatticeState, LatticeStateDefault},
17    },
18    staple, MonteCarlo,
19};
20
21/// Pseudo heat bath algorithm
22///
23/// # Example
24/// ```
25/// # use std::error::Error;
26/// #
27/// # fn main() -> Result<(), Box<dyn Error>> {
28/// use lattice_qcd_rs::simulation::{HeatBathSweep, LatticeState, LatticeStateDefault};
29/// use rand::SeedableRng;
30///
31/// let rng = rand::rngs::StdRng::seed_from_u64(0); // change with your seed
32/// let mut heat_bath = HeatBathSweep::new(rng);
33///
34/// let mut state = LatticeStateDefault::<3>::new_cold(1_f64, 6_f64, 4)?;
35/// for _ in 0..2 {
36///     state = state.monte_carlo_step(&mut heat_bath)?;
37///     // operation to track the progress or the evolution
38/// }
39/// // operation at the end of the simulation
40/// #     Ok(())
41/// # }
42/// ```
43#[derive(Clone, Debug, PartialEq, Eq, Hash)]
44#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
45pub struct HeatBathSweep<Rng: rand::Rng> {
46    rng: Rng,
47}
48
49impl<Rng: rand::Rng> HeatBathSweep<Rng> {
50    /// Create a new Self form a rng.
51    pub const fn new(rng: Rng) -> Self {
52        Self { rng }
53    }
54
55    /// Absorbed self and return the RNG as owned. It essentially deconstruct the structure.
56    #[allow(clippy::missing_const_for_fn)] // false positive
57    pub fn rng_owned(self) -> Rng {
58        self.rng
59    }
60
61    /// Get a mutable reference to the rng.
62    pub fn rng_mut(&mut self) -> &mut Rng {
63        &mut self.rng
64    }
65
66    /// Get a reference to the rng.
67    pub const fn rng(&self) -> &Rng {
68        &self.rng
69    }
70
71    /// Apply te SU2 heat bath method.
72    #[inline]
73    fn heat_bath_su2(&mut self, staple: CMatrix2, beta: f64) -> CMatrix2 {
74        let staple_coefficient = staple.determinant().real().sqrt();
75        if staple_coefficient.is_normal() {
76            let v_r: CMatrix2 = staple.adjoint() / Complex::from(staple_coefficient);
77            let d_heat_bath_r = HeatBathDistribution::new(beta * staple_coefficient).unwrap();
78            let rand_m_r: CMatrix2 = self.rng.sample(d_heat_bath_r);
79            rand_m_r * v_r
80        }
81        else {
82            // just return a random matrix as all matrices
83            // have the same selection probability
84            su2::random_su2(&mut self.rng)
85        }
86    }
87
88    /// Apply the pseudo heat bath methods and return the new link.
89    #[inline]
90    fn get_modif<const D: usize>(
91        &mut self,
92        state: &LatticeStateDefault<D>,
93        link: &LatticeLinkCanonical<D>,
94    ) -> na::Matrix3<Complex> {
95        let link_matrix = state
96            .link_matrix()
97            .matrix(&link.into(), state.lattice())
98            .unwrap();
99        let a = staple(state.link_matrix(), state.lattice(), link);
100
101        let r = su3::get_r(self.heat_bath_su2(su3::get_su2_r_unorm(link_matrix * a), state.beta()));
102        let s =
103            su3::get_s(self.heat_bath_su2(su3::get_su2_s_unorm(r * link_matrix * a), state.beta()));
104        let t = su3::get_t(
105            self.heat_bath_su2(su3::get_su2_t_unorm(s * r * link_matrix * a), state.beta()),
106        );
107
108        t * s * r * link_matrix
109    }
110
111    #[inline]
112    // TODO improve error handling
113    fn next_element_default<const D: usize>(
114        &mut self,
115        mut state: LatticeStateDefault<D>,
116    ) -> LatticeStateDefault<D> {
117        let lattice = state.lattice().clone();
118        lattice.get_links().for_each(|link| {
119            let potential_modif = self.get_modif(&state, &link);
120            *state.link_mut(&link).unwrap() = potential_modif;
121        });
122        state
123    }
124}
125
126impl<Rng: rand::Rng> AsRef<Rng> for HeatBathSweep<Rng> {
127    fn as_ref(&self) -> &Rng {
128        self.rng()
129    }
130}
131
132impl<Rng: rand::Rng> AsMut<Rng> for HeatBathSweep<Rng> {
133    fn as_mut(&mut self) -> &mut Rng {
134        self.rng_mut()
135    }
136}
137
138impl<Rng: rand::Rng + Default> Default for HeatBathSweep<Rng> {
139    fn default() -> Self {
140        Self::new(Rng::default())
141    }
142}
143
144impl<Rng, const D: usize> MonteCarlo<LatticeStateDefault<D>, D> for HeatBathSweep<Rng>
145where
146    Rng: rand::Rng,
147{
148    type Error = Never;
149
150    #[inline]
151    fn next_element(
152        &mut self,
153        state: LatticeStateDefault<D>,
154    ) -> Result<LatticeStateDefault<D>, Self::Error> {
155        Ok(self.next_element_default(state))
156    }
157}
lattice_qcd_rs/simulation/monte_carlo/heat_bath.rs

lattice_qcd_rs/simulation/monte_carlo/
heat_bath.rs
