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

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

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