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

1#[cfg(test)]
2mod test;
3
4use crate::{
5    math::{
6        Scalar, random_uniform, random_x2_normal,
7        special::{erf, erfc},
8    },
9    mechanics::CurrentCoordinate,
10    physics::{
11        BOLTZMANN_CONSTANT,
12        molecular::single_chain::{
13            Configuration, Ensemble, Extensible, Isometric, Isotensional, IsotensionalExtensible,
14            Legendre, MonteCarlo, SingleChain, SingleChainError, Thermodynamics,
15            ThermodynamicsExtensible,
16            ufjc::{
17                // nondimensional_compliance as nondimensional_compliance_asymptotic,
18                nondimensional_extension as nondimensional_extension_asymptotic,
19                nondimensional_gibbs_free_energy_per_link as nondimensional_gibbs_free_energy_per_link_asymptotic,
20            },
21        },
22    },
23};
24use std::f64::consts::{PI, TAU};
25
26/// The extensible freely-jointed chain model.[^1]<sup>,</sup>[^2]
27/// [^1]: N. Balabaev and T. Khazanovich, [Russian Journal of Physical Chemistry B  **3**, 242 (2009)](https://doi.org/10.1134/S1990793109020109).
28/// [^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).
29#[derive(Clone, Debug)]
30pub struct ExtensibleFreelyJointedChain {
31    /// The link length $`\ell_b`$.
32    pub link_length: Scalar,
33    /// The link stiffness $`k_b`$.
34    pub link_stiffness: Scalar,
35    /// The number of links $`N_b`$.
36    pub number_of_links: u8,
37    /// The thermodynamic ensemble.
38    pub ensemble: Ensemble,
39}
40
41impl ExtensibleFreelyJointedChain {
42    fn nondimensional_link_stiffness(&self) -> Scalar {
43        self.link_stiffness * self.link_length().powi(2) / BOLTZMANN_CONSTANT / self.temperature()
44    }
45}
46
47impl SingleChain for ExtensibleFreelyJointedChain {
48    fn link_length(&self) -> Scalar {
49        self.link_length
50    }
51    fn number_of_links(&self) -> u8 {
52        self.number_of_links
53    }
54}
55
56impl Extensible for ExtensibleFreelyJointedChain {}
57
58impl Thermodynamics for ExtensibleFreelyJointedChain {
59    fn ensemble(&self) -> Ensemble {
60        self.ensemble
61    }
62}
63
64impl ThermodynamicsExtensible for ExtensibleFreelyJointedChain {}
65
66impl Isometric for ExtensibleFreelyJointedChain {
67    fn nondimensional_helmholtz_free_energy(
68        &self,
69        _nondimensional_extension: Scalar,
70    ) -> Result<Scalar, SingleChainError> {
71        unimplemented!()
72    }
73    fn nondimensional_force(
74        &self,
75        _nondimensional_extension: Scalar,
76    ) -> Result<Scalar, SingleChainError> {
77        unimplemented!()
78    }
79    fn nondimensional_stiffness(
80        &self,
81        _nondimensional_extension: Scalar,
82    ) -> Result<Scalar, SingleChainError> {
83        unimplemented!()
84    }
85    fn nondimensional_spherical_distribution(
86        &self,
87        _nondimensional_extension: Scalar,
88    ) -> Result<Scalar, SingleChainError> {
89        unimplemented!()
90    }
91}
92
93impl Isotensional for ExtensibleFreelyJointedChain {
94    /// ```math
95    /// \varrho(\eta) = \ln\left[\frac{\eta}{\sinh(\eta)}\right] - \ln\left[1 + \frac{\eta}{\kappa}\,\coth(\eta)\right] - \frac{\eta^2}{2\kappa} - \ln\left[1 + g(\eta)\right]
96    /// ```
97    fn nondimensional_gibbs_free_energy_per_link(
98        &self,
99        nondimensional_force: Scalar,
100    ) -> Result<Scalar, SingleChainError> {
101        let eta = nondimensional_force;
102        let kappa = self.nondimensional_link_stiffness();
103        let eta_over_kappa = eta / kappa;
104        let neg_2_eta_exp = (-2.0 * eta).exp();
105        Ok(nondimensional_gibbs_free_energy_per_link_asymptotic(
106            eta,
107            kappa,
108            -0.5 * eta.powi(2) / kappa,
109            1.0,
110        )? - (0.5
111            + ((eta_over_kappa + 1.0) * erf((eta + kappa) / (2.0 * kappa).sqrt())
112                - (eta_over_kappa - 1.0)
113                    * neg_2_eta_exp
114                    * erf((eta - kappa) / (2.0 * kappa).sqrt()))
115                / (2.0 * (1.0 - neg_2_eta_exp) * (1.0 + eta / eta.tanh() / kappa)))
116            .ln())
117    }
118    /// ```math
119    /// \gamma(\eta) = \mathcal{L}(\eta) + \frac{\eta}{\kappa}\left[\frac{1 - \mathcal{L}(\eta)\coth(\eta)}{1 + (\eta/\kappa)\coth(\eta)}\right] + \frac{\eta}{\kappa} + \frac{g'(\eta)}{1 + g(\eta)}
120    /// ```
121    fn nondimensional_extension(
122        &self,
123        nondimensional_force: Scalar,
124    ) -> Result<Scalar, SingleChainError> {
125        let eta = nondimensional_force;
126        let kappa = self.nondimensional_link_stiffness();
127        let eta_over_kappa = eta / kappa;
128        let neg_2_eta_exp = (-2.0 * eta).exp();
129        let denominator = 2.0 * (1.0 - neg_2_eta_exp) * (1.0 + eta / eta.tanh() / kappa);
130        let fraction = ((eta_over_kappa + 1.0) * erf((eta + kappa) / (2.0 * kappa).sqrt())
131            - (eta_over_kappa - 1.0) * neg_2_eta_exp * erf((eta - kappa) / (2.0 * kappa).sqrt()))
132            / denominator;
133        Ok(
134            nondimensional_extension_asymptotic(eta, kappa, eta_over_kappa, 1.0)?
135                + (((2.0 / PI / kappa).sqrt()
136                    * (eta_over_kappa + 1.0)
137                    * (-(eta + kappa).powi(2) / 2.0 / kappa).exp()
138                    + (1.0 + (1.0 + eta) / kappa))
139                    - 1.0
140                        * neg_2_eta_exp
141                        * ((2.0 / PI / kappa).sqrt()
142                            * (eta_over_kappa - 1.0)
143                            * (-(eta - kappa).powi(2) / 2.0 / kappa).exp()
144                            + (1.0 + (1.0 - eta) / kappa)
145                                * erf((eta - kappa) / (2.0 * kappa).sqrt()))
146                    - fraction
147                        * (2.0
148                            * ((1.0 + neg_2_eta_exp) * (1.0 + (1.0 + eta / eta.tanh()) / kappa)
149                                - 4.0 * eta_over_kappa / (1.0 / neg_2_eta_exp - 1.0))))
150                    / denominator
151                    / (1.0 + fraction),
152        )
153    }
154    /// ```math
155    /// \zeta(\eta) = \mathcal{L}'(\eta) + \frac{\partial}{\partial\eta}\left\{\frac{\eta}{\kappa}\left[\frac{1 - \mathcal{L}(\eta)\coth(\eta)}{1 + (\eta/\kappa)\coth(\eta)}\right]\right\} + \frac{1}{\kappa} + \frac{g''(\eta)}{1 + g(\eta)} - \left[\frac{g'(\eta)}{1 + g(\eta)}\right]^2
156    /// ```
157    fn nondimensional_compliance(
158        &self,
159        _nondimensional_force: Scalar,
160    ) -> Result<Scalar, SingleChainError> {
161        unimplemented!()
162    }
163}
164
165impl IsotensionalExtensible for ExtensibleFreelyJointedChain {
166    /// ```math
167    /// \langle\upsilon\rangle = \frac{\kappa}{2}\Big(\langle\lambda^2\rangle - 2\langle\lambda\rangle + 1\Big)
168    /// ```
169    fn nondimensional_link_energy_average(
170        &self,
171        nondimensional_force: Scalar,
172    ) -> Result<Scalar, SingleChainError> {
173        Ok(0.5
174            * self.nondimensional_link_stiffness()
175            * (nondimensional_link_length_squared_average(
176                self.nondimensional_link_stiffness(),
177                nondimensional_force,
178            )? - 2.0
179                * ThermodynamicsExtensible::nondimensional_link_length_average(
180                    self,
181                    nondimensional_force,
182                )?
183                + 1.0))
184    }
185    /// ```math
186    /// \sigma_\upsilon^2 = \frac{\kappa^2}{4}\Big(\langle\lambda^4\rangle - 4\langle\lambda^3\rangle + 6\langle\lambda^2\rangle - 4\langle\lambda\rangle + 1\Big) - \langle\upsilon\rangle^2
187    /// ```
188    fn nondimensional_link_energy_variance(
189        &self,
190        nondimensional_force: Scalar,
191    ) -> Result<Scalar, SingleChainError> {
192        Ok(0.25
193            * self.nondimensional_link_stiffness().powi(2)
194            * (nondimensional_link_length_quad_average(
195                self.nondimensional_link_stiffness(),
196                nondimensional_force,
197            )? - 4.0
198                * nondimensional_link_length_cubed_average(
199                    self.nondimensional_link_stiffness(),
200                    nondimensional_force,
201                )?
202                + 6.0
203                    * nondimensional_link_length_squared_average(
204                        self.nondimensional_link_stiffness(),
205                        nondimensional_force,
206                    )?
207                - 4.0
208                    * ThermodynamicsExtensible::nondimensional_link_length_average(
209                        self,
210                        nondimensional_force,
211                    )?
212                + 1.0)
213            - ThermodynamicsExtensible::nondimensional_link_energy_average(
214                self,
215                nondimensional_force,
216            )?
217            .powi(2))
218    }
219    /// ```math
220    /// p(\upsilon\,|\,\eta) = \left|\frac{\partial\upsilon}{\partial\lambda}\right|^{-1} \Big[p(\lambda_+\,|\,\eta) + p(\lambda_-\,|\,\eta)\Big]
221    /// ```
222    fn nondimensional_link_energy_probability(
223        &self,
224        nondimensional_energy: Scalar,
225        nondimensional_force: Scalar,
226    ) -> Result<Scalar, SingleChainError> {
227        let kappa = self.nondimensional_link_stiffness();
228        let eta = (2.0 * kappa * nondimensional_energy).sqrt();
229        let delta_lambda = (2.0 * nondimensional_energy / kappa).sqrt();
230        [eta, -eta]
231            .into_iter()
232            .zip([1.0 + delta_lambda, 1.0 - delta_lambda])
233            .map(|(eta, nondimensional_length)| {
234                Ok(
235                    IsotensionalExtensible::nondimensional_link_length_probability(
236                        self,
237                        nondimensional_length,
238                        nondimensional_force,
239                    )? / eta.abs(),
240                )
241            })
242            .sum()
243    }
244    /// ```math
245    /// \langle\lambda\rangle = \frac{\mu_1^+(\kappa,\eta) - \mu_1^-(\kappa,\eta) + \nu_1(\kappa,\eta)}{\mu_0^+(\kappa,\eta) - \mu_0^-(\kappa,\eta)}
246    /// ```
247    fn nondimensional_link_length_average(
248        &self,
249        nondimensional_force: Scalar,
250    ) -> Result<Scalar, SingleChainError> {
251        let eta = nondimensional_force;
252        let kappa = self.nondimensional_link_stiffness();
253        let eta_over_kappa = eta / kappa;
254        let erfd_p = 1.0 + erf((eta + kappa) / (2.0 * kappa).sqrt());
255        let exp_n2_eta_erfc_m = (-2.0 * eta).exp() * erfc((eta - kappa) / (2.0 * kappa).sqrt());
256        Ok(
257            (4.0 * (-0.5 * (eta.powi(2) / kappa + kappa) - eta).exp() / (TAU * kappa).sqrt()
258                * eta_over_kappa
259                + (1.0 / kappa + (eta_over_kappa + 1.0).powi(2)) * erfd_p
260                - (1.0 / kappa + (eta_over_kappa - 1.0).powi(2)) * exp_n2_eta_erfc_m)
261                / ((eta / kappa + 1.0) * erfd_p + (eta / kappa - 1.0) * exp_n2_eta_erfc_m),
262        )
263    }
264    /// ```math
265    /// \sigma_\lambda^2 = \frac{\mu_2^+(\kappa,\eta) - \mu_2^-(\kappa,\eta) + \nu_2^+(\kappa,\eta) - \nu_2^-(\kappa,\eta)}{\mu_0^+(\kappa,\eta) - \mu_0^-(\kappa,\eta)} - \langle\lambda\rangle^2
266    /// ```
267    fn nondimensional_link_length_variance(
268        &self,
269        nondimensional_force: Scalar,
270    ) -> Result<Scalar, SingleChainError> {
271        Ok(nondimensional_link_length_squared_average(
272            self.nondimensional_link_stiffness(),
273            nondimensional_force,
274        )? - ThermodynamicsExtensible::nondimensional_link_length_average(
275            self,
276            nondimensional_force,
277        )?
278        .powi(2))
279    }
280    /// ```math
281    /// p(\lambda\,|\,\eta) = \sqrt{\frac{\kappa}{2\pi}}\,\frac{4\lambda\sinh(\eta\lambda)\,e^{-\upsilon(\lambda)}\,e^{-\eta^2/2\kappa}}{\mu_0^+(\kappa,\eta) - \mu_0^-(\kappa,\eta)}
282    /// ```
283    fn nondimensional_link_length_probability(
284        &self,
285        nondimensional_length: Scalar,
286        nondimensional_force: Scalar,
287    ) -> Result<Scalar, SingleChainError> {
288        let eta = nondimensional_force;
289        let lambda = nondimensional_length;
290        let kappa = self.nondimensional_link_stiffness();
291        let eta_over_kappa = eta / kappa;
292        let upsilon_twice = 0.5 * kappa * ((lambda - 1.0).powi(2) + eta_over_kappa.powi(2));
293        Ok((kappa / TAU).sqrt()
294            * 2.0
295            * lambda
296            * ((eta * (lambda - 1.0) - upsilon_twice).exp()
297                - (-eta * (lambda + 1.0) - upsilon_twice).exp())
298            / ((1.0 + eta_over_kappa) * (1.0 + erf((eta + kappa) / (2.0 * kappa).sqrt()))
299                - (1.0 - eta_over_kappa)
300                    * (-2.0 * eta).exp()
301                    * erfc((eta - kappa) / (2.0 * kappa).sqrt())))
302    }
303}
304
305impl Legendre for ExtensibleFreelyJointedChain {
306    fn nondimensional_spherical_distribution(
307        &self,
308        _nondimensional_extension: Scalar,
309    ) -> Result<Scalar, SingleChainError> {
310        unimplemented!()
311    }
312}
313
314impl MonteCarlo for ExtensibleFreelyJointedChain {
315    fn random_nondimensional_link_vectors(&self, nondimensional_force: Scalar) -> Configuration {
316        let sigma = 1.0 / self.nondimensional_link_stiffness().sqrt();
317        (0..self.number_of_links())
318            .map(|_| {
319                let cos_theta = if nondimensional_force == 0.0 {
320                    2.0 * random_uniform() - 1.0
321                } else {
322                    todo!("Force biases the link stretch too.")
323                };
324                let sin_theta = (1.0 - cos_theta * cos_theta).sqrt();
325                let phi = TAU * random_uniform();
326                let (sin_phi, cos_phi) = phi.sin_cos();
327                let lambda = random_x2_normal(1.0, sigma);
328                CurrentCoordinate::from([
329                    lambda * sin_theta * cos_phi,
330                    lambda * sin_theta * sin_phi,
331                    lambda * cos_theta,
332                ])
333            })
334            .collect()
335    }
336}
337
338fn nondimensional_link_length_squared_average(
339    kappa: Scalar,
340    eta: Scalar,
341) -> Result<Scalar, SingleChainError> {
342    let eta_over_kappa = eta / kappa;
343    let erfd_p_pre = (eta / kappa + 1.0) * (1.0 + erf((eta + kappa) / (2.0 * kappa).sqrt()));
344    let exp_n2_eta_erfc_m_pre =
345        (eta / kappa - 1.0) * (-2.0 * eta).exp() * erfc((eta - kappa) / (2.0 * kappa).sqrt());
346    Ok(
347        (2.0 * (-0.5 * (eta.powi(2) / kappa + kappa) - eta).exp() / (TAU * kappa).sqrt()
348            * ((2.0 / kappa + (eta / kappa + 1.0).powi(2))
349                - (2.0 / kappa + (eta / kappa - 1.0).powi(2)))
350            + (3.0 / kappa + (eta_over_kappa + 1.0).powi(2)) * erfd_p_pre
351            + (3.0 / kappa + (eta_over_kappa - 1.0).powi(2)) * exp_n2_eta_erfc_m_pre)
352            / (erfd_p_pre + exp_n2_eta_erfc_m_pre),
353    )
354}
355
356fn nondimensional_link_length_cubed_average(
357    kappa: Scalar,
358    eta: Scalar,
359) -> Result<Scalar, SingleChainError> {
360    let eta_over_kappa = eta / kappa;
361    let x_p = (eta + kappa) / (2.0 * kappa).sqrt();
362    let x_m = (eta - kappa) / (2.0 * kappa).sqrt();
363    let one_plus_erf_p = 1.0 + erf(x_p);
364    let erfc_m = erfc(x_m);
365    let exp_n2_eta = (-2.0 * eta).exp();
366    let denominator =
367        (eta_over_kappa + 1.0) * one_plus_erf_p + exp_n2_eta * (eta_over_kappa - 1.0) * erfc_m;
368    let p_p = eta.powi(4)
369        + 4.0 * eta.powi(3) * kappa
370        + 6.0 * eta.powi(2) * kappa * (1.0 + kappa)
371        + 4.0 * eta * kappa.powi(2) * (3.0 + kappa)
372        + kappa.powi(2) * (3.0 + 6.0 * kappa + kappa.powi(2));
373    let p_m = eta.powi(4) - 4.0 * eta.powi(3) * kappa + 6.0 * eta.powi(2) * kappa * (1.0 + kappa)
374        - 4.0 * eta * kappa.powi(2) * (3.0 + kappa)
375        + kappa.powi(2) * (3.0 + 6.0 * kappa + kappa.powi(2));
376    let boundary =
377        2.0 * eta * (eta.powi(2) + 5.0 * kappa + 3.0 * kappa.powi(2)) * (2.0 / PI).sqrt()
378            / kappa.powf(3.5)
379            * (-(eta.powi(2) / (2.0 * kappa) + eta + 0.5 * kappa)).exp();
380    let branch_terms = (p_p * one_plus_erf_p - exp_n2_eta * p_m * erfc_m) / kappa.powi(4);
381    Ok((boundary + branch_terms) / denominator)
382}
383
384fn nondimensional_link_length_quad_average(
385    kappa: Scalar,
386    eta: Scalar,
387) -> Result<Scalar, SingleChainError> {
388    let sqrt_kappa = kappa.sqrt();
389    let sqrt_2 = 2.0_f64.sqrt();
390    let sqrt_pi = PI.sqrt();
391    let sqrt_2_pi = (2.0 * PI).sqrt();
392    let x_p = (eta + kappa) / (2.0 * kappa).sqrt();
393    let x_m = (eta - kappa) / (2.0 * kappa).sqrt();
394    let erf_p = erf(x_p);
395    let erfc_m = erfc(x_m);
396    let exp_p = ((eta + kappa).powi(2) / (2.0 * kappa)).exp();
397    let exp_m = ((eta - kappa).powi(2) / (2.0 * kappa)).exp();
398    let exp_n2_eta = (-2.0 * eta).exp();
399    let denominator =
400        (eta / kappa + 1.0) * (1.0 + erf_p) + exp_n2_eta * (eta / kappa - 1.0) * erfc_m;
401    let poly_p = eta.powi(5)
402        + 5.0 * eta.powi(4) * kappa
403        + 10.0 * eta.powi(3) * kappa * (1.0 + kappa)
404        + 10.0 * eta.powi(2) * kappa.powi(2) * (3.0 + kappa)
405        + 5.0 * eta * kappa.powi(2) * (3.0 + 6.0 * kappa + kappa.powi(2))
406        + kappa.powi(3) * (15.0 + 10.0 * kappa + kappa.powi(2));
407    let inner = -2.0 * (eta - kappa).powi(4) * sqrt_kappa
408        - 18.0 * (eta - kappa).powi(2) * kappa.powf(1.5)
409        + 18.0 * kappa.powf(3.5)
410        + 2.0 * kappa.powf(4.5)
411        + 2.0
412            * sqrt_2
413            * eta.powi(3)
414            * kappa
415            * (2.0 * sqrt_2 * sqrt_kappa + 5.0 * exp_p * sqrt_pi + 5.0 * exp_p * kappa * sqrt_pi)
416        + eta.powi(5) * exp_p * sqrt_2_pi
417        + 15.0 * exp_p * kappa.powi(3) * sqrt_2_pi
418        + 10.0 * exp_p * kappa.powi(4) * sqrt_2_pi
419        + exp_p * kappa.powi(5) * sqrt_2_pi
420        + eta.powi(4) * (2.0 * sqrt_kappa + 5.0 * exp_p * kappa * sqrt_2_pi)
421        + 2.0
422            * eta.powi(2)
423            * kappa.powf(1.5)
424            * (9.0
425                + 6.0 * kappa
426                + 15.0 * exp_p * sqrt_kappa * sqrt_2_pi
427                + 5.0 * exp_p * kappa.powf(1.5) * sqrt_2_pi)
428        + eta
429            * kappa.powi(2)
430            * (36.0 * sqrt_kappa
431                + 8.0 * kappa.powf(1.5)
432                + 15.0 * exp_p * sqrt_2_pi
433                + 30.0 * exp_p * kappa * sqrt_2_pi
434                + 5.0 * exp_p * kappa.powi(2) * sqrt_2_pi)
435        + exp_m * (eta - kappa).powi(5) * sqrt_2_pi * erfc_m
436        + 10.0 * exp_m * (eta - kappa).powi(3) * kappa * sqrt_2_pi * erfc_m
437        + 15.0 * exp_m * (eta - kappa) * kappa.powi(2) * sqrt_2_pi * erfc_m
438        + exp_p * poly_p * sqrt_2_pi * erf_p;
439    let numerator = (-(eta.powi(2) / (2.0 * kappa) + eta + 0.5 * kappa)).exp()
440        / (TAU.sqrt() * kappa.powi(5))
441        * inner;
442    Ok(numerator / denominator)
443}
conspire/physics/molecular/single_chain/efjc/mod.rs

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