Skip to main content

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

1#[cfg(test)]
2mod test;
3
4use crate::{
5    math::{Scalar, Tensor, random_uniform},
6    mechanics::CurrentCoordinate,
7    physics::molecular::single_chain::{
8        Configuration, Ensemble, Inextensible, MonteCarlo, SingleChain,
9    },
10};
11use std::f64::consts::TAU;
12
13/// The freely-rotating chain model.
14#[derive(Clone, Debug)]
15pub struct FreelyRotatingChain {
16    /// The link angle $`\theta_b`$.
17    pub link_angle: Scalar,
18    /// The link length $`\ell_b`$.
19    pub link_length: Scalar,
20    /// The number of links $`N_b`$.
21    pub number_of_links: u8,
22    /// The thermodynamic ensemble.
23    pub ensemble: Ensemble,
24}
25
26impl SingleChain for FreelyRotatingChain {
27    fn link_length(&self) -> Scalar {
28        self.link_length
29    }
30    fn number_of_links(&self) -> u8 {
31        self.number_of_links
32    }
33}
34
35impl Inextensible for FreelyRotatingChain {
36    fn maximum_nondimensional_extension(&self) -> Scalar {
37        1.0
38    }
39}
40
41impl MonteCarlo for FreelyRotatingChain {
42    fn random_nondimensional_link_vectors(&self, nondimensional_force: Scalar) -> Configuration {
43        if nondimensional_force != 0.0 {
44            todo!()
45        }
46        let cos_theta = 2.0 * random_uniform() - 1.0;
47        let sin_theta = (1.0 - cos_theta * cos_theta).sqrt();
48        let phi = TAU * random_uniform();
49        let (sin_phi, cos_phi) = phi.sin_cos();
50        const AY: CurrentCoordinate = CurrentCoordinate::const_from([0.0, 1.0, 0.0]);
51        const AZ: CurrentCoordinate = CurrentCoordinate::const_from([0.0, 0.0, 1.0]);
52        let mut a = AY;
53        let mut b =
54            CurrentCoordinate::const_from([sin_theta * cos_phi, sin_theta * sin_phi, cos_theta]);
55        let (sin_theta, cos_theta) = self.link_angle.sin_cos();
56        (0..self.number_of_links())
57            .map(|link| {
58                if link > 0 {
59                    a = if b[1].abs() < 0.9 { AY } else { AZ };
60                    let u = a.cross(&b).normalized();
61                    let v = b.cross(&u);
62                    let phi = TAU * random_uniform();
63                    let (sin_phi, cos_phi) = phi.sin_cos();
64                    b = &b * cos_theta + (&u * cos_phi + &v * sin_phi) * sin_theta;
65                }
66                b.clone()
67            })
68            .collect()
69    }
70}