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

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

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