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

1#[cfg(test)]
2mod test;
3
4use crate::{
5    math::{Scalar, TensorArray, random_uniform},
6    mechanics::CurrentCoordinate,
7    physics::molecular::single_chain::{
8        Configuration, Ensemble, Inextensible, Isometric, Isotensional, Legendre, MonteCarlo,
9        SingleChain, SingleChainError, Thermodynamics,
10    },
11};
12use std::f64::consts::TAU;
13
14/// The square-well freely-jointed chain model.[^1]
15/// [^1]: M.R. Buche, M.N. Silberstein, and S.J. Grutzik, [Physical Review E **106**, 024502 (2022)](https://doi.org/10.1103/PhysRevE.106.024502).
16#[derive(Clone, Debug)]
17pub struct SquareWellFreelyJointedChain {
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 well width $`w_b`$.
23    pub well_width: Scalar,
24    /// The thermodynamic ensemble.
25    pub ensemble: Ensemble,
26}
27
28impl SingleChain for SquareWellFreelyJointedChain {
29    fn link_length(&self) -> Scalar {
30        self.link_length
31    }
32    fn number_of_links(&self) -> u8 {
33        self.number_of_links
34    }
35}
36
37impl Inextensible for SquareWellFreelyJointedChain {
38    /// ```math
39    /// \lim_{\eta\to\infty}\gamma(\eta) = 1 + \frac{w_b}{\ell_b} = \varsigma
40    /// ```
41    fn maximum_nondimensional_extension(&self) -> Scalar {
42        1.0 + self.well_width / self.link_length
43    }
44}
45
46impl Thermodynamics for SquareWellFreelyJointedChain {
47    fn ensemble(&self) -> Ensemble {
48        self.ensemble
49    }
50}
51
52impl Isometric for SquareWellFreelyJointedChain {
53    fn nondimensional_helmholtz_free_energy(
54        &self,
55        _nondimensional_extension: Scalar,
56    ) -> Result<Scalar, SingleChainError> {
57        unimplemented!()
58    }
59    fn nondimensional_force(
60        &self,
61        _nondimensional_extension: Scalar,
62    ) -> Result<Scalar, SingleChainError> {
63        unimplemented!()
64    }
65    fn nondimensional_stiffness(
66        &self,
67        _nondimensional_extension: Scalar,
68    ) -> Result<Scalar, SingleChainError> {
69        unimplemented!()
70    }
71    fn nondimensional_spherical_distribution(
72        &self,
73        _nondimensional_extension: Scalar,
74    ) -> Result<Scalar, SingleChainError> {
75        unimplemented!()
76    }
77}
78
79impl Isotensional for SquareWellFreelyJointedChain {
80    /// ```math
81    /// \beta\varphi(\eta) = N_b\ln\left[\frac{\eta^3}{\varsigma\eta\cosh(\varsigma\eta) - \sinh(\varsigma\eta) - \eta\cosh(\eta) + \sinh(\eta)}\right]
82    /// ```
83    fn nondimensional_gibbs_free_energy(
84        &self,
85        nondimensional_force: Scalar,
86    ) -> Result<Scalar, SingleChainError> {
87        let varsigma = self.maximum_nondimensional_extension();
88        let varsigma_eta = varsigma * nondimensional_force;
89        Ok(self.number_of_links() as Scalar
90            * (nondimensional_force.powi(3)
91                / (varsigma_eta * varsigma_eta.cosh()
92                    - varsigma_eta.sinh()
93                    - nondimensional_force * nondimensional_force.cosh()
94                    + nondimensional_force.sinh()))
95            .ln())
96    }
97    /// ```math
98    /// \gamma(\eta) = \frac{\varsigma^2\eta\sinh(\varsigma\eta) - \eta\sinh(\eta)}{\varsigma\eta\cosh(\varsigma\eta) - \sinh(\varsigma\eta) - \eta\cosh(\eta) + \sinh(\eta)} - \frac{3}{\eta}
99    /// ```
100    fn nondimensional_extension(
101        &self,
102        nondimensional_force: Scalar,
103    ) -> Result<Scalar, SingleChainError> {
104        if nondimensional_force == 0.0 {
105            Ok(0.0)
106        } else {
107            let eta = nondimensional_force;
108            let eta_sinh = eta.sinh();
109            let varsigma = self.maximum_nondimensional_extension();
110            let varsigma_eta = varsigma * eta;
111            let varsigma_eta_sinh = varsigma_eta.sinh();
112            Ok(eta * (varsigma.powi(2) * varsigma_eta_sinh - eta_sinh)
113                / (varsigma_eta * varsigma_eta.cosh() - varsigma_eta_sinh - eta * eta.cosh()
114                    + eta_sinh)
115                - 3.0 / eta)
116        }
117    }
118    /// ```math
119    /// \zeta(\eta) = \frac{\left(\varsigma^2\sinh(\varsigma\eta)+\varsigma^3\eta\cosh(\varsigma\eta)-\sinh(\eta)-\eta\cosh(\eta)\right)\left(\varsigma\eta\cosh(\varsigma\eta)-\sinh(\varsigma\eta)-\eta\cosh(\eta)+\sinh(\eta)\right)-\left(\varsigma^2\eta\sinh(\varsigma\eta)-\eta\sinh(\eta)\right)^2}{\left(\varsigma\eta\cosh(\varsigma\eta)-\sinh(\varsigma\eta)-\eta\cosh(\eta)+\sinh(\eta)\right)^2}+\frac{3}{\eta^2}
120    /// ```
121    fn nondimensional_compliance(
122        &self,
123        nondimensional_force: Scalar,
124    ) -> Result<Scalar, SingleChainError> {
125        if nondimensional_force == 0.0 {
126            Ok(Scalar::NAN)
127        } else {
128            let eta = nondimensional_force;
129            let eta_sinh = eta.sinh();
130            let eta_cosh = eta.cosh();
131            let varsigma = self.maximum_nondimensional_extension();
132            let varsigma_eta = varsigma * nondimensional_force;
133            let varsigma_eta_sinh = varsigma_eta.sinh();
134            let varsigma_eta_cosh = varsigma_eta.cosh();
135            let a = eta * (varsigma * varsigma * varsigma_eta_sinh - eta_sinh);
136            let b =
137                varsigma_eta * varsigma_eta_cosh - varsigma_eta_sinh - eta * eta_cosh + eta_sinh;
138            let a_prime = (varsigma * varsigma * varsigma_eta_sinh - eta_sinh)
139                + eta * (varsigma * varsigma * varsigma * varsigma_eta_cosh - eta_cosh);
140            Ok((a_prime * b - a * a) / (b * b) + 3.0 / (eta * eta))
141        }
142    }
143}
144
145impl Legendre for SquareWellFreelyJointedChain {
146    fn nondimensional_spherical_distribution(
147        &self,
148        _nondimensional_extension: Scalar,
149    ) -> Result<Scalar, SingleChainError> {
150        unimplemented!()
151    }
152}
153
154impl MonteCarlo for SquareWellFreelyJointedChain {
155    fn random_configuration(&self) -> Configuration {
156        let mut position = CurrentCoordinate::zero();
157        let max_strain = self.maximum_nondimensional_extension() - 1.0;
158        (0..self.number_of_links())
159            .map(|_| {
160                let cos_theta = 2.0 * random_uniform() - 1.0;
161                let sin_theta = (1.0 - cos_theta * cos_theta).sqrt();
162                let phi = TAU * random_uniform();
163                let (sin_phi, cos_phi) = phi.sin_cos();
164                let lambda = 1.0 + max_strain * random_uniform();
165                position[0] += lambda * sin_theta * cos_phi;
166                position[1] += lambda * sin_theta * sin_phi;
167                position[2] += lambda * cos_theta;
168                position.clone()
169            })
170            .collect()
171    }
172}
conspire/physics/molecular/single_chain/swfjc/mod.rs

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