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, sinhc},
8    },
9    physics::molecular::{
10        potential::Potential,
11        single_chain::{
12            Ensemble, Isometric, Isotensional, Legendre, SingleChain, SingleChainError,
13            Thermodynamics,
14        },
15    },
16};
17
18/// The freely-jointed chain model with an arbitrary link potential.[^1]
19/// [^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).
20#[derive(Clone, Debug)]
21pub struct ArbitraryPotentialFreelyJointedChain<T>
22where
23    T: Potential,
24{
25    /// The link potential $`u`$.
26    pub link_potential: T,
27    /// The number of links $`N_b`$.
28    pub number_of_links: u8,
29    /// The thermodynamic ensemble.
30    pub ensemble: Ensemble,
31}
32
33impl<T> ArbitraryPotentialFreelyJointedChain<T>
34where
35    T: Potential,
36{
37    fn correction(&self) -> Scalar {
38        1.0 / (1.0
39            - 0.5
40                * self
41                    .link_potential
42                    .nondimensional_anharmonicity(1.0, self.temperature())
43                / self
44                    .link_potential
45                    .nondimensional_stiffness(1.0, self.temperature()))
46    }
47    fn nondimensional_link_stiffness(&self) -> Scalar {
48        self.link_potential
49            .nondimensional_stiffness(0.0, self.temperature())
50    }
51}
52
53impl<T> SingleChain for ArbitraryPotentialFreelyJointedChain<T>
54where
55    T: Potential,
56{
57    fn link_length(&self) -> Scalar {
58        self.link_potential.rest_length()
59    }
60    fn number_of_links(&self) -> u8 {
61        self.number_of_links
62    }
63}
64
65impl<T> Thermodynamics for ArbitraryPotentialFreelyJointedChain<T>
66where
67    T: Potential,
68{
69    fn ensemble(&self) -> Ensemble {
70        self.ensemble
71    }
72}
73
74impl<T> Isometric for ArbitraryPotentialFreelyJointedChain<T>
75where
76    T: Potential,
77{
78    fn nondimensional_helmholtz_free_energy(
79        &self,
80        _nondimensional_extension: Scalar,
81    ) -> Result<Scalar, SingleChainError> {
82        unimplemented!()
83    }
84    fn nondimensional_force(
85        &self,
86        _nondimensional_extension: Scalar,
87    ) -> Result<Scalar, SingleChainError> {
88        unimplemented!()
89    }
90    fn nondimensional_stiffness(
91        &self,
92        _nondimensional_extension: Scalar,
93    ) -> Result<Scalar, SingleChainError> {
94        unimplemented!()
95    }
96    fn nondimensional_spherical_distribution(
97        &self,
98        _nondimensional_extension: Scalar,
99    ) -> Result<Scalar, SingleChainError> {
100        unimplemented!()
101    }
102}
103
104impl<T> Isotensional for ArbitraryPotentialFreelyJointedChain<T>
105where
106    T: Potential,
107{
108    /// ```math
109    /// \varrho(\eta) = \ln\left[\frac{\eta}{\sinh(\eta)}\right] - \ln\left[1 + \frac{\eta}{c\kappa}\,\coth(\eta)\right] - \beta v(\eta)
110    /// ```
111    fn nondimensional_gibbs_free_energy_per_link(
112        &self,
113        nondimensional_force: Scalar,
114    ) -> Result<Scalar, SingleChainError> {
115        nondimensional_gibbs_free_energy_per_link(
116            nondimensional_force,
117            self.nondimensional_link_stiffness(),
118            self.link_potential
119                .nondimensional_legendre(nondimensional_force, self.temperature()),
120            self.correction(),
121        )
122    }
123    /// ```math
124    /// \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)
125    /// ```
126    fn nondimensional_extension(
127        &self,
128        nondimensional_force: Scalar,
129    ) -> Result<Scalar, SingleChainError> {
130        nondimensional_extension(
131            nondimensional_force,
132            self.nondimensional_link_stiffness(),
133            self.link_potential
134                .nondimensional_extension(nondimensional_force, self.temperature()),
135            self.correction(),
136        )
137    }
138    /// ```math
139    /// \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)
140    /// ```
141    fn nondimensional_compliance(
142        &self,
143        nondimensional_force: Scalar,
144    ) -> Result<Scalar, SingleChainError> {
145        nondimensional_compliance(
146            nondimensional_force,
147            self.nondimensional_link_stiffness(),
148            self.link_potential
149                .nondimensional_compliance(nondimensional_force, self.temperature()),
150            self.correction(),
151        )
152    }
153}
154
155impl<T> Legendre for ArbitraryPotentialFreelyJointedChain<T>
156where
157    T: Potential,
158{
159    fn nondimensional_spherical_distribution(
160        &self,
161        _nondimensional_extension: Scalar,
162    ) -> Result<Scalar, SingleChainError> {
163        unimplemented!()
164    }
165}
166
167pub fn nondimensional_gibbs_free_energy_per_link(
168    eta: Scalar,
169    kappa: Scalar,
170    beta_v: Scalar,
171    c: Scalar,
172) -> Result<Scalar, SingleChainError> {
173    Ok(-((sinhc(eta) * (1.0 + eta / c / kappa / eta.tanh())).ln() - beta_v))
174}
175
176pub fn nondimensional_extension(
177    eta: Scalar,
178    kappa: Scalar,
179    delta_lambda: Scalar,
180    c: Scalar,
181) -> Result<Scalar, SingleChainError> {
182    if eta == 0.0 {
183        Ok(0.0)
184    } else {
185        let eta_coth = 1.0 / eta.tanh();
186        let gamma_0 = langevin(eta);
187        let eta_over_kappa = eta / kappa;
188        Ok(gamma_0
189            + eta_over_kappa * (1.0 - gamma_0 * eta_coth) / (c + eta_over_kappa * eta_coth)
190            + delta_lambda)
191    }
192}
193
194pub fn nondimensional_compliance(
195    eta: Scalar,
196    kappa: Scalar,
197    zeta: Scalar,
198    c: Scalar,
199) -> Result<Scalar, SingleChainError> {
200    if eta == 0.0 {
201        Ok(1.0 / 3.0 + 2.0 / 3.0 / c / kappa + zeta)
202    } else {
203        let eta_tanh = eta.tanh();
204        let eta_coth = 1.0 / eta_tanh;
205        let gamma_0 = langevin(eta);
206        let eta_over_kappa = eta / kappa;
207        let c_0 = langevin_derivative(eta);
208        let g = 1.0 - gamma_0 * eta_coth;
209        let h = c + eta_over_kappa * eta_coth;
210        let dcth = 1.0 - 1.0 / (eta_tanh * eta_tanh);
211        let dg = -(c_0 * eta_coth + gamma_0 * dcth);
212        let dh = eta_coth / kappa + eta_over_kappa * dcth;
213        Ok(c_0 + (g / h) / kappa + eta_over_kappa * (dg * h - g * dh) / (h * h) + zeta)
214    }
215}
conspire/physics/molecular/single_chain/ufjc/mod.rs

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