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

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

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