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

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

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