conspire/physics/molecular/single_chain/efrc/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::{
5    math::{Scalar, Tensor, random_uniform, random_x2_normal},
6    mechanics::CurrentCoordinate,
7    physics::{
8        BOLTZMANN_CONSTANT,
9        molecular::single_chain::{
10            Configuration, Ensemble, Extensible, Isometric, Isotensional, Legendre, MonteCarlo,
11            SingleChain, SingleChainError, Thermodynamics,
12        },
13    },
14};
15use std::f64::consts::TAU;
16
17/// The extensible freely-rotating chain model.
18#[derive(Clone, Debug)]
19pub struct ExtensibleFreelyRotatingChain {
20    /// The link angle $`\theta_b`$.
21    pub link_angle: Scalar,
22    /// The link length $`\ell_b`$.
23    pub link_length: Scalar,
24    /// The link stiffness $`k_b`$.
25    pub link_stiffness: Scalar,
26    /// The number of links $`N_b`$.
27    pub number_of_links: u8,
28    /// The thermodynamic ensemble.
29    pub ensemble: Ensemble,
30}
31
32impl ExtensibleFreelyRotatingChain {
33    fn nondimensional_link_stiffness(&self) -> Scalar {
34        self.link_stiffness * self.link_length().powi(2) / BOLTZMANN_CONSTANT / self.temperature()
35    }
36}
37
38impl SingleChain for ExtensibleFreelyRotatingChain {
39    fn link_length(&self) -> Scalar {
40        self.link_length
41    }
42    fn number_of_links(&self) -> u8 {
43        self.number_of_links
44    }
45}
46
47impl Extensible for ExtensibleFreelyRotatingChain {}
48
49impl Thermodynamics for ExtensibleFreelyRotatingChain {
50    fn ensemble(&self) -> Ensemble {
51        self.ensemble
52    }
53}
54
55impl Isometric for ExtensibleFreelyRotatingChain {
56    fn nondimensional_helmholtz_free_energy(
57        &self,
58        _nondimensional_extension: Scalar,
59    ) -> Result<Scalar, SingleChainError> {
60        unimplemented!()
61    }
62    fn nondimensional_force(
63        &self,
64        _nondimensional_extension: Scalar,
65    ) -> Result<Scalar, SingleChainError> {
66        unimplemented!()
67    }
68    fn nondimensional_stiffness(
69        &self,
70        _nondimensional_extension: Scalar,
71    ) -> Result<Scalar, SingleChainError> {
72        unimplemented!()
73    }
74    fn nondimensional_spherical_distribution(
75        &self,
76        _nondimensional_extension: Scalar,
77    ) -> Result<Scalar, SingleChainError> {
78        unimplemented!()
79    }
80}
81
82impl Isotensional for ExtensibleFreelyRotatingChain {
83    fn nondimensional_gibbs_free_energy_per_link(
84        &self,
85        _nondimensional_force: Scalar,
86    ) -> Result<Scalar, SingleChainError> {
87        unimplemented!()
88    }
89    fn nondimensional_extension(
90        &self,
91        _nondimensional_force: Scalar,
92    ) -> Result<Scalar, SingleChainError> {
93        unimplemented!()
94    }
95    fn nondimensional_compliance(
96        &self,
97        _nondimensional_force: Scalar,
98    ) -> Result<Scalar, SingleChainError> {
99        unimplemented!()
100    }
101}
102
103impl Legendre for ExtensibleFreelyRotatingChain {
104    fn nondimensional_spherical_distribution(
105        &self,
106        _nondimensional_extension: Scalar,
107    ) -> Result<Scalar, SingleChainError> {
108        unimplemented!()
109    }
110}
111
112impl MonteCarlo for ExtensibleFreelyRotatingChain {
113    fn nondimensional_longitudinal_extension(
114        &self,
115        nondimensional_force: Scalar,
116        number_of_samples: usize,
117        number_of_threads: usize,
118    ) -> Scalar {
119        nondimensional_extension_reweighted_biased_stretch(
120            self,
121            nondimensional_force,
122            nondimensional_force,
123            number_of_samples,
124            number_of_threads,
125        )
126    }
127    fn random_nondimensional_link_vectors(&self, nondimensional_force: Scalar) -> Configuration {
128        if nondimensional_force != 0.0 {
129            unimplemented!()
130        }
131        let std = 1.0 / self.nondimensional_link_stiffness().sqrt();
132        let cos_theta = 2.0 * random_uniform() - 1.0;
133        let sin_theta = (1.0 - cos_theta * cos_theta).sqrt();
134        let phi = TAU * random_uniform();
135        let (sin_phi, cos_phi) = phi.sin_cos();
136        const AY: CurrentCoordinate = CurrentCoordinate::const_from([0.0, 1.0, 0.0]);
137        const AZ: CurrentCoordinate = CurrentCoordinate::const_from([0.0, 0.0, 1.0]);
138        let mut a = AY;
139        let mut b =
140            CurrentCoordinate::const_from([sin_theta * cos_phi, sin_theta * sin_phi, cos_theta]);
141        let (sin_theta, cos_theta) = self.link_angle.sin_cos();
142        (0..self.number_of_links())
143            .map(|link| {
144                if link > 0 {
145                    a = if b[1].abs() < 0.9 { AY } else { AZ };
146                    let u = a.cross(&b).normalized();
147                    let v = b.cross(&u);
148                    let phi = TAU * random_uniform();
149                    let (sin_phi, cos_phi) = phi.sin_cos();
150                    b = &b * cos_theta + (&u * cos_phi + &v * sin_phi) * sin_theta;
151                }
152                &b * random_x2_normal(1.0, std)
153            })
154            .collect()
155    }
156}
157
158fn random_nondimensional_link_vectors_biased_stretch(
159    model: &ExtensibleFreelyRotatingChain,
160    nondimensional_stretch_bias: Scalar,
161) -> Configuration {
162    let kappa = model.nondimensional_link_stiffness();
163    let std = 1.0 / kappa.sqrt();
164    let mean = 1.0 + nondimensional_stretch_bias / kappa;
165
166    let cos_theta = 2.0 * random_uniform() - 1.0;
167    let sin_theta = (1.0 - cos_theta * cos_theta).sqrt();
168    let phi = TAU * random_uniform();
169    let (sin_phi, cos_phi) = phi.sin_cos();
170
171    const AY: CurrentCoordinate = CurrentCoordinate::const_from([0.0, 1.0, 0.0]);
172    const AZ: CurrentCoordinate = CurrentCoordinate::const_from([0.0, 0.0, 1.0]);
173
174    let mut a = AY;
175    let mut b =
176        CurrentCoordinate::const_from([sin_theta * cos_phi, sin_theta * sin_phi, cos_theta]);
177
178    let (sin_theta, cos_theta) = model.link_angle.sin_cos();
179
180    (0..model.number_of_links())
181        .map(|link| {
182            if link > 0 {
183                a = if b[1].abs() < 0.9 { AY } else { AZ };
184                let u = a.cross(&b).normalized();
185                let v = b.cross(&u);
186                let phi = TAU * random_uniform();
187                let (sin_phi, cos_phi) = phi.sin_cos();
188                b = &b * cos_theta + (&u * cos_phi + &v * sin_phi) * sin_theta;
189            }
190            &b * random_x2_normal(mean, std)
191        })
192        .collect()
193}
194
195use std::thread::scope;
196
197fn nondimensional_extension_reweighted_biased_stretch(
198    model: &ExtensibleFreelyRotatingChain,
199    nondimensional_force: Scalar,
200    nondimensional_stretch_bias: Scalar,
201    number_of_samples: usize,
202    number_of_threads: usize,
203) -> Scalar {
204    let base = number_of_samples / number_of_threads;
205    let remainder = number_of_samples % number_of_threads;
206
207    scope(|s| {
208        (0..number_of_threads)
209            .map(|t| {
210                s.spawn(move || {
211                    nondimensional_extension_reweighted_biased_stretch_inner(
212                        model,
213                        nondimensional_force,
214                        nondimensional_stretch_bias,
215                        base + usize::from(t < remainder),
216                    )
217                })
218            })
219            .collect::<Vec<_>>()
220            .into_iter()
221            .map(|handle| handle.join().unwrap())
222            .reduce(|mut acc, (x_max, z_scaled, ext_scaled)| {
223                let x_max_new = acc.0.max(x_max);
224                let scale_acc = (acc.0 - x_max_new).exp();
225                let scale_new = (x_max - x_max_new).exp();
226
227                acc.1 = acc.1 * scale_acc + z_scaled * scale_new;
228                acc.2 = acc.2 * scale_acc + ext_scaled * scale_new;
229                acc.0 = x_max_new;
230                acc
231            })
232            .map(|(_x_max, z_scaled, ext_scaled)| {
233                ext_scaled / z_scaled / model.number_of_links() as Scalar
234            })
235            .unwrap()
236    })
237}
238
239fn nondimensional_extension_reweighted_biased_stretch_inner(
240    model: &ExtensibleFreelyRotatingChain,
241    nondimensional_force: Scalar,
242    nondimensional_stretch_bias: Scalar,
243    number_of_samples: usize,
244) -> (Scalar, Scalar, Scalar) {
245    let mut x_max = Scalar::NEG_INFINITY;
246    let mut z_scaled = 0.0;
247    let mut ext_scaled = 0.0;
248
249    for _ in 0..number_of_samples {
250        let links =
251            random_nondimensional_link_vectors_biased_stretch(model, nondimensional_stretch_bias);
252
253        let extension_sum: Scalar = links.iter().map(|link| link[2]).sum();
254        let stretch_sum: Scalar = links.iter().map(|link| link.norm()).sum();
255
256        let x = nondimensional_force * extension_sum - nondimensional_stretch_bias * stretch_sum;
257
258        if x > x_max {
259            let scale = if x_max.is_finite() {
260                (x_max - x).exp()
261            } else {
262                0.0
263            };
264            z_scaled *= scale;
265            ext_scaled *= scale;
266            x_max = x;
267        }
268
269        let w = (x - x_max).exp();
270        z_scaled += w;
271        ext_scaled += extension_sum * w;
272    }
273
274    (x_max, z_scaled, ext_scaled)
275}
conspire/physics/molecular/single_chain/efrc/mod.rs

conspire/physics/molecular/single_chain/efrc/
mod.rs
