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

1#[cfg(test)]
2mod test;
3
4use crate::{
5    math::{
6        Quantity, Scalar,
7        special::{langevin, langevin_derivative},
8    },
9    physics::molecular::{
10        potential::{Harmonic, Potential},
11        single_chain::{
12            Ensemble, Extensible, Isometric, Isotensional, IsotensionalExtensible, Legendre,
13            SingleChain, SingleChainError, Thermodynamics, ThermodynamicsExtensible,
14        },
15    },
16    units::Length,
17};
18use std::f64::consts::TAU;
19
20/// The freely-jointed chain model with an arbitrary link potential.[^1]
21/// [^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).
22#[derive(Clone, Debug)]
23pub struct ArbitraryPotentialFreelyJointedChain<T>
24where
25    T: Potential,
26{
27    /// The link potential $`u`$.
28    pub link_potential: T,
29    /// The number of links $`N_b`$.
30    pub number_of_links: u8,
31    /// The thermodynamic ensemble.
32    pub ensemble: Ensemble,
33}
34
35impl<T> ArbitraryPotentialFreelyJointedChain<T>
36where
37    T: Potential,
38{
39    fn correction(&self) -> Scalar {
40        1.0 / (1.0
41            - 0.5
42                * self
43                    .link_potential
44                    .nondimensional_anharmonicity(1.0, self.temperature())
45                / self
46                    .link_potential
47                    .nondimensional_stiffness(1.0, self.temperature()))
48    }
49    fn nondimensional_link_stiffness(&self) -> Scalar {
50        self.link_potential
51            .nondimensional_stiffness(1.0, self.temperature())
52    }
53}
54
55impl<T> SingleChain for ArbitraryPotentialFreelyJointedChain<T>
56where
57    T: Potential,
58{
59    fn link_length(&self) -> Quantity<Length> {
60        self.link_potential.rest_length()
61    }
62    fn number_of_links(&self) -> u8 {
63        self.number_of_links
64    }
65}
66
67impl<T> Extensible for ArbitraryPotentialFreelyJointedChain<T> where T: Potential {}
68
69impl<T> Thermodynamics for ArbitraryPotentialFreelyJointedChain<T>
70where
71    T: Potential,
72{
73    fn ensemble(&self) -> Ensemble {
74        self.ensemble
75    }
76}
77
78impl ThermodynamicsExtensible for ArbitraryPotentialFreelyJointedChain<Harmonic> {}
79
80impl<T> Isometric for ArbitraryPotentialFreelyJointedChain<T>
81where
82    T: Potential,
83{
84    fn nondimensional_helmholtz_free_energy(
85        &self,
86        _nondimensional_extension: Scalar,
87    ) -> Result<Scalar, SingleChainError> {
88        unimplemented!()
89    }
90    fn nondimensional_force(
91        &self,
92        _nondimensional_extension: Scalar,
93    ) -> Result<Scalar, SingleChainError> {
94        unimplemented!()
95    }
96    fn nondimensional_stiffness(
97        &self,
98        _nondimensional_extension: Scalar,
99    ) -> Result<Scalar, SingleChainError> {
100        unimplemented!()
101    }
102    fn nondimensional_spherical_distribution(
103        &self,
104        _nondimensional_extension: Scalar,
105    ) -> Result<Scalar, SingleChainError> {
106        unimplemented!()
107    }
108}
109
110impl<T> Isotensional for ArbitraryPotentialFreelyJointedChain<T>
111where
112    T: Potential,
113{
114    /// ```math
115    /// \varrho(\eta) = \ln\left[\frac{\eta}{\sinh(\eta)}\right] - \ln\left[1 + \frac{\eta}{c\kappa}\,\coth(\eta)\right] - \nu(\eta)
116    /// ```
117    fn nondimensional_gibbs_free_energy_per_link(
118        &self,
119        nondimensional_force: Scalar,
120    ) -> Result<Scalar, SingleChainError> {
121        nondimensional_gibbs_free_energy_per_link(
122            nondimensional_force,
123            self.nondimensional_link_stiffness(),
124            self.link_potential
125                .nondimensional_legendre(nondimensional_force, self.temperature()),
126            self.correction(),
127        )
128    }
129    /// ```math
130    /// \gamma(\eta) = \mathcal{L}(\eta) + \frac{\eta}{\kappa}\left[\frac{1 - \mathcal{L}(\eta)\coth(\eta)}{c + (\eta/\kappa)\coth(\eta)}\right] + \Delta\lambda(\eta)
131    /// ```
132    fn nondimensional_extension(
133        &self,
134        nondimensional_force: Scalar,
135    ) -> Result<Scalar, SingleChainError> {
136        nondimensional_extension(
137            nondimensional_force,
138            self.nondimensional_link_stiffness(),
139            self.link_potential
140                .nondimensional_extension(nondimensional_force, self.temperature()),
141            self.correction(),
142        )
143    }
144    /// ```math
145    /// \zeta(\eta) = \mathcal{L}'(\eta) + \frac{\partial}{\partial\eta}\left\{\frac{\eta}{\kappa}\left[\frac{1 - \mathcal{L}(\eta)\coth(\eta)}{c + (\eta/\kappa)\coth(\eta)}\right]\right\} + \zeta(\eta)
146    /// ```
147    fn nondimensional_compliance(
148        &self,
149        nondimensional_force: Scalar,
150    ) -> Result<Scalar, SingleChainError> {
151        nondimensional_compliance(
152            nondimensional_force,
153            self.nondimensional_link_stiffness(),
154            self.link_potential
155                .nondimensional_compliance(nondimensional_force, self.temperature()),
156            self.correction(),
157        )
158    }
159}
160
161impl IsotensionalExtensible for ArbitraryPotentialFreelyJointedChain<Harmonic> {
162    /// ```math
163    /// \langle\upsilon\rangle = \frac{1}{2} + \frac{(\eta/\kappa)\coth(\eta)}{c + (\eta/\kappa)\coth(\eta)} + \upsilon[\lambda(\eta)]
164    /// ```
165    fn nondimensional_link_energy_average(
166        &self,
167        nondimensional_force: Scalar,
168    ) -> Result<Scalar, SingleChainError> {
169        Ok(0.5
170            + helper(
171                nondimensional_force,
172                self.nondimensional_link_stiffness(),
173                self.correction(),
174            )
175            + self
176                .link_potential
177                .nondimensional_energy_at_nondimensional_force(
178                    nondimensional_force,
179                    self.temperature(),
180                ))
181    }
182    /// ```math
183    /// \sigma_\upsilon^2 = \frac{1}{2} + \frac{(\eta/\kappa)\coth(\eta)}{c + (\eta/\kappa)\coth(\eta)}\left[2 - \frac{(\eta/\kappa)\coth(\eta)}{c + (\eta/\kappa)\coth(\eta)}\right] + 2\upsilon[\lambda(\eta)]
184    /// ```
185    fn nondimensional_link_energy_variance(
186        &self,
187        nondimensional_force: Scalar,
188    ) -> Result<Scalar, SingleChainError> {
189        let hlpr = helper(
190            nondimensional_force,
191            self.nondimensional_link_stiffness(),
192            self.correction(),
193        );
194        Ok(0.5
195            + hlpr * (2.0 - hlpr)
196            + 2.0
197                * self
198                    .link_potential
199                    .nondimensional_energy_at_nondimensional_force(
200                        nondimensional_force,
201                        self.temperature(),
202                    ))
203    }
204    /// ```math
205    /// p(\upsilon\,|\,\eta) = \left|\frac{\partial\upsilon}{\partial\lambda}\right|^{-1} \Big[p(\lambda_+\,|\,\eta) + p(\lambda_-\,|\,\eta)\Big]
206    /// ```
207    fn nondimensional_link_energy_probability(
208        &self,
209        nondimensional_energy: Scalar,
210        nondimensional_force: Scalar,
211    ) -> Result<Scalar, SingleChainError> {
212        self.link_potential
213            .nondimensional_forces_at_nondimensional_energy(
214                nondimensional_energy,
215                self.temperature(),
216            )
217            .into_iter()
218            .zip(
219                self.link_potential
220                    .nondimensional_lengths_at_nondimensional_energy(
221                        nondimensional_energy,
222                        self.temperature(),
223                    ),
224            )
225            .map(|(eta, nondimensional_length)| {
226                Ok(
227                    IsotensionalExtensible::nondimensional_link_length_probability(
228                        self,
229                        nondimensional_length,
230                        nondimensional_force,
231                    )? / eta.abs(),
232                )
233            })
234            .sum()
235    }
236    /// ```math
237    /// \langle\lambda\rangle = 1 + \frac{1/\kappa + (\eta/\kappa)(1 - \eta/\kappa)(\coth(\eta) - 1)}{1 + (\eta/\kappa)\coth(\eta)} + \Delta\lambda(\eta)
238    /// ```
239    fn nondimensional_link_length_average(
240        &self,
241        nondimensional_force: Scalar,
242    ) -> Result<Scalar, SingleChainError> {
243        let eta = nondimensional_force;
244        let kappa = self.nondimensional_link_stiffness();
245        if eta == 0.0 {
246            Ok(1.0 + 2.0 / (1.0 * kappa + 1.0))
247        } else {
248            let eta_coth = 1.0 / eta.tanh();
249            let eta_over_kappa = eta / kappa;
250            Ok(1.0
251                + (1.0 / kappa + eta_over_kappa * (1.0 - eta_over_kappa) * (eta_coth - 1.0))
252                    / (1.0 + eta_over_kappa * eta_coth)
253                + eta_over_kappa)
254        }
255    }
256    /// ```math
257    /// \sigma_\lambda^2 = 1 + \frac{3/\kappa + 2\eta^2/\kappa^2 + (3/\kappa + 2)(\eta/\kappa)\coth(\eta)}{1 + (\eta/\kappa)\coth(\eta)} + \Delta\lambda^2(\eta) - \langle\lambda\rangle^2
258    /// ```
259    fn nondimensional_link_length_variance(
260        &self,
261        nondimensional_force: Scalar,
262    ) -> Result<Scalar, SingleChainError> {
263        let eta = nondimensional_force;
264        let kappa = self.nondimensional_link_stiffness();
265        let mean_squared =
266            ThermodynamicsExtensible::nondimensional_link_length_average(self, eta)?.powi(2);
267        if eta == 0.0 {
268            Ok(1.0 + 3.0 / kappa + 2.0 / (kappa + 1.0) - mean_squared)
269        } else {
270            let eta_coth = 1.0 / eta.tanh();
271            let eta_over_kappa = eta / kappa;
272            let eta_over_kappa_coth = eta_over_kappa * eta_coth;
273            Ok(1.0
274                + (3.0 / kappa
275                    + 2.0 * eta_over_kappa.powi(2)
276                    + (3.0 / kappa + 2.0) * eta_over_kappa_coth)
277                    / (1.0 + eta_over_kappa_coth)
278                + eta_over_kappa.powi(2)
279                - mean_squared)
280        }
281    }
282    /// ```math
283    /// p(\lambda\,|\,\eta) = \sqrt{\frac{\kappa}{2\pi}}\,\frac{\mathrm{sinhc}(\eta\lambda)}{\mathrm{sinhc}(\eta)}\,\frac{\lambda^2\,e^{-\upsilon(\lambda)}\,e^{-\eta^2/2\kappa}}{1 + (\eta/c\kappa)\coth(\eta)}
284    /// ```
285    fn nondimensional_link_length_probability(
286        &self,
287        nondimensional_length: Scalar,
288        nondimensional_force: Scalar,
289    ) -> Result<Scalar, SingleChainError> {
290        let eta = nondimensional_force;
291        let lambda = nondimensional_length;
292        let kappa = self.nondimensional_link_stiffness();
293        let upsilon_twice = self
294            .link_potential
295            .nondimensional_energy(nondimensional_length, self.temperature())
296            + eta.powi(2) / 2.0 / kappa;
297        Ok((kappa / TAU).sqrt()
298            * lambda
299            * ((eta * (lambda - 1.0) - upsilon_twice).exp()
300                - (-eta * (lambda + 1.0) - upsilon_twice).exp())
301            / (1.0 - (-2.0 * eta).exp())
302            / (1.0 + eta / kappa / self.correction() / eta.tanh()))
303    }
304}
305
306pub(super) fn nondimensional_gibbs_free_energy_per_link(
307    eta: Scalar,
308    kappa: Scalar,
309    nu: Scalar,
310    c: Scalar,
311) -> Result<Scalar, SingleChainError> {
312    Ok(nu
313        - eta
314        - (0.5 - 0.5 * (-2.0 * eta).exp()).ln()
315        - (1.0 / eta + 1.0 / c / kappa / eta.tanh()).ln())
316}
317
318pub(super) fn nondimensional_extension(
319    eta: Scalar,
320    kappa: Scalar,
321    delta_lambda: Scalar,
322    c: Scalar,
323) -> Result<Scalar, SingleChainError> {
324    if eta == 0.0 {
325        Ok(0.0)
326    } else {
327        let eta_coth = 1.0 / eta.tanh();
328        let gamma_0 = langevin(eta);
329        let eta_over_kappa = eta / kappa;
330        Ok(gamma_0
331            + eta_over_kappa * (1.0 - gamma_0 * eta_coth) / (c + eta_over_kappa * eta_coth)
332            + delta_lambda)
333    }
334}
335
336pub(super) fn nondimensional_compliance(
337    eta: Scalar,
338    kappa: Scalar,
339    zeta: Scalar,
340    c: Scalar,
341) -> Result<Scalar, SingleChainError> {
342    if eta == 0.0 {
343        Ok(1.0 / 3.0 + 2.0 / 3.0 / c / kappa + zeta)
344    } else {
345        let eta_tanh = eta.tanh();
346        let eta_coth = 1.0 / eta_tanh;
347        let gamma_0 = langevin(eta);
348        let eta_over_kappa = eta / kappa;
349        let c_0 = langevin_derivative(eta);
350        let g = 1.0 - gamma_0 * eta_coth;
351        let h = c + eta_over_kappa * eta_coth;
352        let dcth = 1.0 - 1.0 / (eta_tanh * eta_tanh);
353        let dg = -(c_0 * eta_coth + gamma_0 * dcth);
354        let dh = eta_coth / kappa + eta_over_kappa * dcth;
355        Ok(c_0 + (g / h) / kappa + eta_over_kappa * (dg * h - g * dh) / (h * h) + zeta)
356    }
357}
358
359fn helper(
360    nondimensional_force: Scalar,
361    nondimensional_stiffness: Scalar,
362    correction: Scalar,
363) -> Scalar {
364    let eta_over_kappa = nondimensional_force / nondimensional_stiffness;
365    eta_over_kappa / (eta_over_kappa + correction * nondimensional_force.tanh())
366}
367
368impl<T> Legendre for ArbitraryPotentialFreelyJointedChain<T>
369where
370    T: Potential,
371{
372    fn nondimensional_spherical_distribution(
373        &self,
374        _nondimensional_extension: Scalar,
375    ) -> Result<Scalar, SingleChainError> {
376        unimplemented!()
377    }
378}
conspire/physics/molecular/single_chain/ufjc/mod.rs

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