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

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

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