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

1#[cfg(test)]
2mod test;
3
4use crate::{
5    math::{Scalar, TensorArray, random_uniform, random_x2_normal, special::erf},
6    mechanics::CurrentCoordinate,
7    physics::{
8        BOLTZMANN_CONSTANT,
9        molecular::single_chain::{
10            Configuration, Ensemble, Extensible, Isometric, Isotensional, Legendre, MonteCarlo,
11            SingleChain, SingleChainError, Thermodynamics,
12            ufjc::{
13                // nondimensional_compliance as nondimensional_compliance_asymptotic,
14                nondimensional_extension as nondimensional_extension_asymptotic,
15                nondimensional_gibbs_free_energy_per_link as nondimensional_gibbs_free_energy_per_link_asymptotic,
16            },
17        },
18    },
19};
20use std::f64::consts::{PI, TAU};
21
22/// The extensible freely-jointed chain model.[^1]<sup>,</sup>[^2]
23/// [^1]: N. Balabaev and T. Khazanovich, [Russian Journal of Physical Chemistry B  **3**, 242 (2009)](https://doi.org/10.1134/S1990793109020109).
24/// [^2]: M.R. Buche, M.N. Silberstein, and S.J. Grutzik, [Physical Review E **106**, 024502 (2022)](https://doi.org/10.1103/PhysRevE.106.024502).
25#[derive(Clone, Debug)]
26pub struct ExtensibleFreelyJointedChain {
27    /// The link length $`\ell_b`$.
28    pub link_length: Scalar,
29    /// The link stiffness $`k_b`$.
30    pub link_stiffness: Scalar,
31    /// The number of links $`N_b`$.
32    pub number_of_links: u8,
33    /// The thermodynamic ensemble.
34    pub ensemble: Ensemble,
35}
36
37impl ExtensibleFreelyJointedChain {
38    fn nondimensional_link_stiffness(&self) -> Scalar {
39        self.link_stiffness * self.link_length().powi(2) / BOLTZMANN_CONSTANT / self.temperature()
40    }
41}
42
43impl SingleChain for ExtensibleFreelyJointedChain {
44    fn link_length(&self) -> Scalar {
45        self.link_length
46    }
47    fn number_of_links(&self) -> u8 {
48        self.number_of_links
49    }
50}
51
52impl Extensible for ExtensibleFreelyJointedChain {}
53
54impl Thermodynamics for ExtensibleFreelyJointedChain {
55    fn ensemble(&self) -> Ensemble {
56        self.ensemble
57    }
58}
59
60impl Isometric for ExtensibleFreelyJointedChain {
61    fn nondimensional_helmholtz_free_energy(
62        &self,
63        _nondimensional_extension: Scalar,
64    ) -> Result<Scalar, SingleChainError> {
65        unimplemented!()
66    }
67    fn nondimensional_force(
68        &self,
69        _nondimensional_extension: Scalar,
70    ) -> Result<Scalar, SingleChainError> {
71        unimplemented!()
72    }
73    fn nondimensional_stiffness(
74        &self,
75        _nondimensional_extension: Scalar,
76    ) -> Result<Scalar, SingleChainError> {
77        unimplemented!()
78    }
79    fn nondimensional_spherical_distribution(
80        &self,
81        _nondimensional_extension: Scalar,
82    ) -> Result<Scalar, SingleChainError> {
83        unimplemented!()
84    }
85}
86
87impl Isotensional for ExtensibleFreelyJointedChain {
88    /// ```math
89    /// \varrho(\eta) = \varrho_a(\eta) - \ln\left[1 + g(\eta)\right]
90    /// ```
91    fn nondimensional_gibbs_free_energy_per_link(
92        &self,
93        nondimensional_force: Scalar,
94    ) -> Result<Scalar, SingleChainError> {
95        let eta = nondimensional_force;
96        let kappa = self.nondimensional_link_stiffness();
97        let eta_over_kappa = eta / kappa;
98        let eta_exp = eta.exp();
99        Ok(nondimensional_gibbs_free_energy_per_link_asymptotic(
100            eta,
101            kappa,
102            -0.5 * eta.powi(2) / kappa,
103            1.0,
104        )? - (0.5
105            + ((eta_over_kappa + 1.0) * eta_exp * erf(&((eta + kappa) / (2.0 * kappa).sqrt()))
106                - (eta_over_kappa - 1.0) / eta_exp * erf(&((eta - kappa) / (2.0 * kappa).sqrt())))
107                / (4.0 * eta.sinh() * (1.0 + eta / eta.tanh() / kappa)))
108            .ln())
109    }
110    /// ```math
111    /// \gamma(\eta) = \gamma_a(\eta) + \frac{g'(\eta)}{1 + g(\eta)}
112    /// ```
113    fn nondimensional_extension(
114        &self,
115        nondimensional_force: Scalar,
116    ) -> Result<Scalar, SingleChainError> {
117        let eta = nondimensional_force;
118        let kappa = self.nondimensional_link_stiffness();
119        let eta_over_kappa = eta / kappa;
120        let denominator = 4.0 * eta.sinh() * (1.0 + eta / eta.tanh() / kappa);
121        let fraction = ((eta_over_kappa + 1.0)
122            * eta.exp()
123            * erf(&((eta + kappa) / (2.0 * kappa).sqrt()))
124            - (eta_over_kappa - 1.0) / eta.exp() * erf(&((eta - kappa) / (2.0 * kappa).sqrt())))
125            / denominator;
126        Ok(
127            nondimensional_extension_asymptotic(eta, kappa, eta_over_kappa, 1.0)?
128                + (eta.exp()
129                    * ((2.0 / PI / kappa).sqrt()
130                        * (eta_over_kappa + 1.0)
131                        * (-(eta + kappa).powi(2) / 2.0 / kappa).exp()
132                        + (1.0 + (1.0 + eta) / kappa))
133                    - 1.0 / eta.exp()
134                        * ((2.0 / PI / kappa).sqrt()
135                            * (eta_over_kappa - 1.0)
136                            * (-(eta - kappa).powi(2) / 2.0 / kappa).exp()
137                            + (1.0 + (1.0 - eta) / kappa)
138                                * erf(&((eta - kappa) / (2.0 * kappa).sqrt())))
139                    - fraction
140                        * (4.0
141                            * (eta.cosh() * (1.0 + (1.0 + eta / eta.tanh()) / kappa)
142                                - eta_over_kappa / eta.sinh())))
143                    / denominator
144                    / (1.0 + fraction),
145        )
146    }
147    /// ```math
148    /// \zeta(\eta) = \zeta_a(\eta) + \frac{g''(\eta)}{1 + g(\eta)} - \left[\frac{g'(\eta)}{1 + g(\eta)}\right]^2
149    /// ```
150    fn nondimensional_compliance(
151        &self,
152        _nondimensional_force: Scalar,
153    ) -> Result<Scalar, SingleChainError> {
154        unimplemented!()
155    }
156}
157
158impl Legendre for ExtensibleFreelyJointedChain {
159    fn nondimensional_spherical_distribution(
160        &self,
161        _nondimensional_extension: Scalar,
162    ) -> Result<Scalar, SingleChainError> {
163        unimplemented!()
164    }
165}
166
167impl MonteCarlo for ExtensibleFreelyJointedChain {
168    fn random_configuration(&self) -> Configuration {
169        let mut position = CurrentCoordinate::zero();
170        let std = 1.0 / self.nondimensional_link_stiffness().sqrt();
171        (0..self.number_of_links())
172            .map(|_| {
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                let link_stretch = random_x2_normal(1.0, std);
178                position[0] += link_stretch * sin_theta * cos_phi;
179                position[1] += link_stretch * sin_theta * sin_phi;
180                position[2] += link_stretch * cos_theta;
181                position.clone()
182            })
183            .collect()
184    }
185}
conspire/physics/molecular/single_chain/efjc/mod.rs

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