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

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

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