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

1#[cfg(test)]
2mod test;
3
4use crate::math::Current;
5use crate::{
6    math::{
7        Quantity, Scalar,
8        random::random_uniform,
9        special::{inverse_langevin, langevin, langevin_derivative, sinhc},
10    },
11    mechanics::Vector,
12    physics::molecular::single_chain::{
13        Configuration, Ensemble, Inextensible, Isometric, Isotensional, Legendre, MonteCarlo,
14        SingleChain, SingleChainError, Thermodynamics,
15    },
16    units::Length,
17};
18use std::f64::consts::{PI, TAU};
19
20/// The freely-jointed chain model.
21#[derive(Clone, Debug)]
22pub struct FreelyJointedChain {
23    /// The link length $`\ell_b`$.
24    pub link_length: Scalar,
25    /// The number of links $`N_b`$.
26    pub number_of_links: u8,
27    /// The thermodynamic ensemble.
28    pub ensemble: Ensemble,
29}
30
31impl SingleChain for FreelyJointedChain {
32    fn link_length(&self) -> Quantity<Length> {
33        Quantity::new(self.link_length)
34    }
35    fn number_of_links(&self) -> u8 {
36        self.number_of_links
37    }
38}
39
40impl Inextensible for FreelyJointedChain {
41    /// ```math
42    /// \lim_{\eta\to\infty}\gamma(\eta) = 1
43    /// ```
44    fn maximum_nondimensional_extension(&self) -> Scalar {
45        1.0
46    }
47}
48
49impl Thermodynamics for FreelyJointedChain {
50    fn ensemble(&self) -> Ensemble {
51        self.ensemble
52    }
53}
54
55impl Isometric for FreelyJointedChain {
56    fn nondimensional_helmholtz_free_energy(
57        &self,
58        nondimensional_extension: Scalar,
59    ) -> Result<Scalar, SingleChainError> {
60        self.nondimensional_extension_check(nondimensional_extension)?;
61        if nondimensional_extension == 0.0 {
62            Ok(0.0)
63        } else {
64            let [s0, _, _] = treloar_sums(self.number_of_links(), nondimensional_extension);
65            Ok(nondimensional_extension.abs().ln() - s0.ln())
66        }
67    }
68    /// ```math
69    /// \eta(\gamma) = \frac{1}{N_b\gamma} + \left(\frac{1}{2} - \frac{1}{N_b}\right)\frac{\sum_{s=0}^{s_\mathrm{max}}(-1)^s\binom{N_b}{s}\left(m - \frac{s}{N_b}\right)^{N_b - 3}}{\sum_{s=0}^{s_\mathrm{max}}(-1)^s\binom{N_b}{s}\left(m - \frac{s}{N_b}\right)^{N_b - 2}}
70    /// ```
71    fn nondimensional_force(
72        &self,
73        nondimensional_extension: Scalar,
74    ) -> Result<Scalar, SingleChainError> {
75        self.nondimensional_extension_check(nondimensional_extension)?;
76        if nondimensional_extension == 0.0 {
77            Ok(0.0)
78        } else {
79            let [s0, s1, _] = treloar_sums(self.number_of_links(), nondimensional_extension);
80            let n = self.number_of_links() as Scalar;
81            Ok((1.0 / nondimensional_extension + (0.5 * n - 1.0) * s1 / s0) / n)
82        }
83    }
84    /// ```math
85    /// \kappa(\gamma) = \frac{\partial\eta}{\partial\gamma}
86    /// ```
87    fn nondimensional_stiffness(
88        &self,
89        nondimensional_extension: Scalar,
90    ) -> Result<Scalar, SingleChainError> {
91        self.nondimensional_extension_check(nondimensional_extension)?;
92        if nondimensional_extension == 0.0 {
93            Ok(Scalar::NAN)
94        } else {
95            let [s0, s1, s2] = treloar_sums(self.number_of_links(), nondimensional_extension);
96            if !s0.is_finite() || s0 == 0.0 {
97                return Ok(Scalar::NAN);
98            }
99            let n = self.number_of_links() as Scalar;
100            let p = n - 2.0;
101            let b = (0.5 * n - 1.0) / n;
102            let ds0dx = -(p / 2.0) * s1;
103            let ds1dx = -((p - 1.0) / 2.0) * s2;
104            let d_ratio_dx = (ds1dx * s0 - s1 * ds0dx) / (s0 * s0);
105            Ok(-1.0 / (n * nondimensional_extension * nondimensional_extension) + b * d_ratio_dx)
106        }
107    }
108    /// ```math
109    /// \mathcal{P}(\gamma) = \frac{1}{8\pi\gamma}\frac{N_b^{N_b}}{(N_b - 2)!}\sum_{s=0}^{s_\mathrm{max}}(-1)^s\binom{N_b}{s}\left(m - \frac{s}{N_b}\right)^{N_b - 2}
110    /// ```
111    fn nondimensional_spherical_distribution(
112        &self,
113        nondimensional_extension: Scalar,
114    ) -> Result<Scalar, SingleChainError> {
115        self.nondimensional_extension_check(nondimensional_extension)?;
116        if nondimensional_extension <= 0.0 || nondimensional_extension >= 1.0 {
117            Ok(0.0)
118        } else {
119            let number_of_links = self.number_of_links();
120            let [s0, _, _] = treloar_sums(number_of_links, nondimensional_extension);
121            let n = number_of_links as Scalar;
122            let factorial_n_minus_2 = (1..=(number_of_links - 2))
123                .map(|i| i as Scalar)
124                .product::<Scalar>();
125            Ok((n.powf(n) / (8.0 * PI * nondimensional_extension * factorial_n_minus_2)) * s0)
126        }
127    }
128}
129
130impl Isotensional for FreelyJointedChain {
131    /// ```math
132    /// \varrho(\eta) = N_b\ln\left[\frac{\eta}{\sinh(\eta)}\right]
133    /// ```
134    fn nondimensional_gibbs_free_energy_per_link(
135        &self,
136        nondimensional_force: Scalar,
137    ) -> Result<Scalar, SingleChainError> {
138        Ok(-sinhc(nondimensional_force).ln())
139    }
140    /// ```math
141    /// \gamma(\eta) = \mathcal{L}(\eta)
142    /// ```
143    fn nondimensional_extension(
144        &self,
145        nondimensional_force: Scalar,
146    ) -> Result<Scalar, SingleChainError> {
147        Ok(langevin(nondimensional_force))
148    }
149    /// ```math
150    /// \zeta(\eta) = \mathcal{L}'(\eta)
151    /// ```
152    fn nondimensional_compliance(
153        &self,
154        nondimensional_force: Scalar,
155    ) -> Result<Scalar, SingleChainError> {
156        Ok(langevin_derivative(nondimensional_force))
157    }
158}
159
160impl Legendre for FreelyJointedChain {
161    /// ```math
162    /// \eta(\gamma) = \mathcal{L}^{-1}(\gamma)
163    /// ```
164    fn nondimensional_force(
165        &self,
166        nondimensional_extension: Scalar,
167    ) -> Result<Scalar, SingleChainError> {
168        self.nondimensional_extension_check(nondimensional_extension)?;
169        Ok(inverse_langevin(nondimensional_extension))
170    }
171    /// ```math
172    /// \mathcal{P}(\gamma) \propto \left\{\frac{\sinh[\eta(\gamma)]}{\eta(\gamma)\exp[\eta(\gamma)\gamma]}\right\}^{N_b}
173    /// ```
174    fn nondimensional_spherical_distribution(
175        &self,
176        nondimensional_extension: Scalar,
177    ) -> Result<Scalar, SingleChainError> {
178        let nondimensional_force = Legendre::nondimensional_force(self, nondimensional_extension)?;
179        Ok(
180            (((nondimensional_force * (1.0 - nondimensional_extension)).exp()
181                - (-nondimensional_force * (1.0 + nondimensional_extension)).exp())
182                / 2.0
183                / nondimensional_force)
184                .powi(self.number_of_links() as i32)
185                / normalization(self.number_of_links()),
186        )
187    }
188}
189
190impl MonteCarlo for FreelyJointedChain {
191    fn random_nondimensional_link_vectors(&self, nondimensional_force: Scalar) -> Configuration {
192        let eta = nondimensional_force;
193        let eta_exp = eta.exp();
194        let eta_nexp = 1.0 / eta_exp;
195        (0..self.number_of_links())
196            .map(|_| {
197                let cos_theta = if eta == 0.0 {
198                    2.0 * random_uniform() - 1.0
199                } else {
200                    (eta_nexp + random_uniform() * (eta_exp - eta_nexp)).ln() / eta
201                };
202                let sin_theta = (1.0 - cos_theta * cos_theta).sqrt();
203                let phi = TAU * random_uniform();
204                let (sin_phi, cos_phi) = phi.sin_cos();
205                Vector::<Current>::from([sin_theta * cos_phi, sin_theta * sin_phi, cos_theta])
206            })
207            .collect()
208    }
209}
210
211fn treloar_sums(number_of_links: u8, x: Scalar) -> [Scalar; 3] {
212    if number_of_links <= 2 {
213        return [Scalar::NAN; 3];
214    }
215    let n = number_of_links as Scalar;
216    let p = (number_of_links - 2) as i32;
217    let m = 0.5 * (1.0 - x);
218    let k = ((n * m).ceil() as usize)
219        .saturating_sub(1)
220        .min(number_of_links as usize);
221    let k_float = n * m;
222    if (k_float - k_float.round()).abs() == 0.0 {
223        return [Scalar::NAN; 3];
224    }
225    let mut binom = 1.0;
226    let mut s0 = 0.0;
227    let mut s1 = 0.0;
228    let mut s2 = 0.0;
229    for s in 0..=k {
230        let sign = if s % 2 == 0 { 1.0 } else { -1.0 };
231        let t = m - (s as Scalar) / n;
232        let t0 = if p >= 0 {
233            t.powi(p)
234        } else if t == 0.0 {
235            0.0
236        } else {
237            t.powi(p)
238        };
239        let t1 = if p > 0 {
240            t.powi(p - 1)
241        } else if t == 0.0 {
242            0.0
243        } else {
244            t.powi(p - 1)
245        };
246        let t2 = if p > 1 {
247            t.powi(p - 2)
248        } else if t == 0.0 {
249            0.0
250        } else {
251            t.powi(p - 2)
252        };
253        s0 += sign * binom * t0;
254        s1 += sign * binom * t1;
255        s2 += sign * binom * t2;
256        let sf = s as Scalar;
257        binom *= (n - sf) / (sf + 1.0);
258    }
259    [s0, s1, s2]
260}
261
262fn normalization(number_of_links: u8) -> Scalar {
263    match number_of_links {
264        0 => Scalar::NAN,
265        1 => 1.389_063_303_837_301_3,
266        2 => 0.714_480_944_477_587_6,
267        3 => 0.446_182_225_454_993_8,
268        4 => 0.310_582_574_239_989_03,
269        5 => 0.231_583_731_936_937_35,
270        6 => 0.181_026_390_997_248_38,
271        7 => 0.146_444_713_993_307_1,
272        8 => 0.121_590_329_098_661_26,
273        9 => 0.103_031_548_251_807_95,
274        10 => 0.088_746_746_615_799_2,
275        11 => 0.077_477_021_054_147_71,
276        12 => 0.068_402_348_281_970_37,
277        13 => 0.060_968_329_153_341_53,
278        14 => 0.054_788_235_109_506_34,
279        15 => 0.049_585_008_268_986_38,
280        16 => 0.045_155_543_187_723_454,
281        17 => 0.041_347_934_607_350_624,
282        18 => 0.038_046_552_454_682_33,
283        19 => 0.035_161_995_757_729_23,
284        20 => 0.032_624_174_659_916_245,
285        21 => 0.030_377_448_781_842_47,
286        22 => 0.028_377_147_909_992_65,
287        23 => 0.026_587_040_753_179_49,
288        24 => 0.024_977_465_826_024_475,
289        25 => 0.023_523_932_434_435_773,
290        26 => 0.022_206_060_476_346_62,
291        27 => 0.021_006_767_816_333_316,
292        28 => 0.019_911_640_864_367_884,
293        29 => 0.018_908_442_314_802_338,
294        30 => 0.017_986_722_687_273_728,
295        31 => 0.017_137_511_214_426_318,
296        32 => 0.016_353_067_950_360_97,
297        33 => 0.015_626_683_526_651_468,
298        34 => 0.014_952_516_294_544_284,
299        35 => 0.014_325_459_026_026_574,
300        36 => 0.013_741_029_152_869_741,
301        37 => 0.013_195_277_875_651_702,
302        38 => 0.012_684_714_496_711_45,
303        39 => 0.012_206_243_109_212_051,
304        40 => 0.011_757_109_371_641_886,
305        41 => 0.011_334_855_558_612_255,
306        42 => 0.010_937_282_437_959_286,
307        43 => 0.010_562_416_805_450_496,
308        44 => 0.010_208_483_730_069_894,
309        45 => 0.009_873_882_738_568_507,
310        46 => 0.009_557_167_308_025_053,
311        47 => 0.009_257_027_147_393_918,
312        48 => 0.008_972_272_839_407_26,
313        49 => 0.008_701_822_487_349_997,
314        50 => 0.008_444_690_070_700_62,
315        51 => 0.008_199_975_262_200_628,
316        52 => 0.007_966_854_498_747_187,
317        53 => 0.007_744_573_131_302_788,
318        54 => 0.007_532_438_506_130_219,
319        55 => 0.007_329_813_852_159_101,
320        56 => 0.007_136_112_868_025_883,
321        57 => 0.006_950_794_917_986_018,
322        58 => 0.006_773_360_759_023_209,
323        59 => 0.006_603_348_732_523_006,
324        60 => 0.006_440_331_363_193_850_5,
325        61 => 0.006_283_912_315_803_12,
326        62 => 0.006_133_723_666_987_055,
327        63 => 0.005_989_423_455_088_637,
328        64 => 0.005_850_693_475_837_429,
329        65 => 0.005_717_237_295_843_889,
330        66 => 0.005_588_778_459_447_456,
331        67 => 0.005_465_058_867_524_899,
332        68 => 0.005_345_837_309_508_891,
333        69 => 0.005_230_888_132_150_507,
334        70 => 0.005_120_000_030_536_583,
335        71 => 0.005_012_974_948_588_448,
336        72 => 0.004_909_627_077_760_272,
337        73 => 0.004_809_781_943_954_87,
338        74 => 0.004_713_275_573_809_406_5,
339        75 => 0.004_619_953_732_495_746,
340        76 => 0.004_529_671_226_049_816,
341        77 => 0.004_442_291_262_007_673_5,
342        78 => 0.004_357_684_862_797_343,
343        79 => 0.004_275_730_326_926_852,
344        80 => 0.004_196_312_733_530_766,
345        81 => 0.004_119_323_486_298_74,
346        82 => 0.004_044_659_893_217_922,
347        83 => 0.003_972_224_778_923_046,
348        84 => 0.003_901_926_126_769_502,
349        85 => 0.003_833_676_748_030_511_5,
350        86 => 0.003_767_393_975_874_093_7,
351        87 => 0.003_702_999_382_002_534_4,
352        88 => 0.003_640_418_514_039_760_3,
353        89 => 0.003_579_580_651_933_304,
354        90 => 0.003_520_418_581_799_841,
355        91 => 0.003_462_868_385_788_78,
356        92 => 0.003_406_869_246_668_994,
357        93 => 0.003_352_363_265_961_134_8,
358        94 => 0.003_299_295_294_543_597_2,
359        95 => 0.003_247_612_774_755_336_3,
360        96 => 0.003_197_265_593_104_502_5,
361        97 => 0.003_148_205_942_769_372,
362        98 => 0.003_100_388_195_147_996,
363        99 => 0.003_053_768_779_776_389_4,
364        100 => 0.003_008_306_071_992_423_4,
365        101 => 0.002_963_960_287_774_609_6,
366        102 => 0.002_920_693_385_232_166,
367        103 => 0.002_878_468_972_265_675,
368        104 => 0.002_837_252_219_956_577_4,
369        105 => 0.002_797_009_781_279_318,
370        106 => 0.002_757_709_714_762_243_4,
371        107 => 0.002_719_321_412_752_858,
372        108 => 0.002_681_815_533_969_964_8,
373        109 => 0.002_645_163_940_049_785,
374        110 => 0.002_609_339_635_815_636,
375        111 => 0.002_574_316_713_021_305,
376        112 => 0.002_540_070_297_337_108_2,
377        113 => 0.002_506_576_498_364_846,
378        114 => 0.002_473_812_362_483_738,
379        115 => 0.002_441_755_828_343_935,
380        116 => 0.002_410_385_684_837_54,
381        117 => 0.002_379_681_531_389_384_6,
382        118 => 0.002_349_623_740_421_053,
383        119 => 0.002_320_193_421_852_075,
384        120 => 0.002_291_372_389_511_768,
385        121 => 0.002_263_143_129_344_052_4,
386        122 => 0.002_235_488_769_295_701_7,
387        123 => 0.002_208_393_050_786_021,
388        124 => 0.002_181_840_301_662_884_3,
389        125 => 0.002_155_815_410_556_500_4,
390        126 => 0.002_130_303_802_548_212_4,
391        127 => 0.002_105_291_416_077_128,
392        128 => 0.002_080_764_681_012_503,
393        129 => 0.002_056_710_497_824_479_8,
394        130 => 0.002_033_116_217_790_202_3,
395        131 => 0.002_009_969_624_176_372,
396        132 => 0.001_987_258_914_343_074,
397        133 => 0.001_964_972_682_717_241_4,
398        134 => 0.001_943_099_904_587_340_1,
399        135 => 0.001_921_629_920_673_915_4,
400        136 => 0.001_900_552_422_433_467_6,
401        137 => 0.001_879_857_438_055_702_6,
402        138 => 0.001_859_535_319_116_718_1,
403        139 => 0.001_839_576_727_852_897_4,
404        140 => 0.001_819_972_625_022_462_4,
405        141 => 0.001_800_714_258_323_577,
406        142 => 0.001_781_793_151_339_777_3,
407        143 => 0.001_763_201_092_985_212,
408        144 => 0.001_744_930_127_423_801_3,
409        145 => 0.001_726_972_544_437_928,
410        146 => 0.001_709_320_870_223_689,
411        147 => 0.001_691_967_858_591_060_8,
412        148 => 0.001_674_906_482_548_559_4,
413        149 => 0.001_658_129_926_253_140_5,
414        150 => 0.001_641_631_577_307_177_5,
415        151 => 0.001_625_405_019_385_355_6,
416        152 => 0.001_609_444_025_175_287_7,
417        153 => 0.001_593_742_549_616_558,
418        154 => 0.001_578_294_723_423_718_3,
419        155 => 0.001_563_094_846_879_572_4,
420        156 => 0.001_548_137_383_885_819,
421        157 => 0.001_533_416_956_258_811_4,
422        158 => 0.001_518_928_338_258_862_9,
423        159 => 0.001_504_666_451_342_125,
424        160 => 0.001_490_626_359_124_665_4,
425        161 => 0.001_476_803_262_548_898,
426        162 => 0.001_463_192_495_243_045,
427        163 => 0.001_449_789_519_064_795,
428        164 => 0.001_436_589_919_820_758_2,
429        165 => 0.001_423_589_403_153_785_4,
430        166 => 0.001_410_783_790_590_583,
431        167 => 0.001_398_169_015_742_466_7,
432        168 => 0.001_385_741_120_652_447_7,
433        169 => 0.001_373_496_252_282_184_2,
434        170 => 0.001_361_430_659_132_658_5,
435        171 => 0.001_349_540_687_992_734_2,
436        172 => 0.001_337_822_780_810_052_7,
437        173 => 0.001_326_273_471_678_944_4,
438        174 => 0.001_314_889_383_940_535_3,
439        175 => 0.001_303_667_227_389_772_3,
440        176 => 0.001_292_603_795_585_479_2,
441        177 => 0.001_281_695_963_258_609_4,
442        178 => 0.001_270_940_683_814_773_4,
443        179 => 0.001_260_334_986_927_068_6,
444        180 => 0.001_249_875_976_215_474_5,
445        181 => 0.001_239_560_827_009_231_7,
446        182 => 0.001_229_386_784_188_803_6,
447        183 => 0.001_219_351_160_104_175_2,
448        184 => 0.001_209_451_332_566_384_2,
449        185 => 0.001_199_684_742_909_333_8,
450        186 => 0.001_190_048_894_119_060_9,
451        187 => 0.001_180_541_349_027_766_8,
452        188 => 0.001_171_159_728_570_035,
453        189 => 0.001_161_901_710_098_784_4,
454        190 => 0.001_152_765_025_758_600_2,
455        191 => 0.001_143_747_460_914_206_6,
456        192 => 0.001_134_846_852_631_931_8,
457        193 => 0.001_126_061_088_212_113_2,
458        194 => 0.001_117_388_103_770_484_4,
459        195 => 0.001_108_825_882_866_665,
460        196 => 0.001_100_372_455_177_960_5,
461        197 => 0.001_092_025_895_216_747_5,
462        198 => 0.001_083_784_321_089_808_8,
463        199 => 0.001_075_645_893_298_036_3,
464        200 => 0.001_067_608_813_574_995_4,
465        201 => 0.001_059_671_323_762_909_2,
466        202 => 0.001_051_831_704_724_672_7,
467        203 => 0.001_044_088_275_290_578_8,
468        204 => 0.001_036_439_391_238_473_6,
469        205 => 0.001_028_883_444_306_137,
470        206 => 0.001_021_418_861_234_703,
471        207 => 0.001_014_044_102_842_014_4,
472        208 => 0.001_006_757_663_124_827_5,
473        209 => 0.000_999_558_068_388_836_5,
474        210 => 0.000_992_443_876_405_530_2,
475        211 => 0.000_985_413_675_594_929_3,
476        212 => 0.000_978_466_084_233_29,
477        213 => 0.000_971_599_749_684_902_6,
478        214 => 0.000_964_813_347_657_139_1,
479        215 => 0.000_958_105_581_477_943_6,
480        216 => 0.000_951_475_181_394_988_1,
481        217 => 0.000_944_920_903_895_750_9,
482        218 => 0.000_938_441_531_047_793_8,
483        219 => 0.000_932_035_869_858_555_4,
484        220 => 0.000_925_702_751_653_993_6,
485        221 => 0.000_919_441_031_475_440_3,
486        222 => 0.000_913_249_587_494_056_1,
487        223 => 0.000_907_127_320_442_295_6,
488        224 => 0.000_901_073_153_061_812_5,
489        225 => 0.000_895_086_029_567_262_1,
490        226 => 0.000_889_164_915_125_473_2,
491        227 => 0.000_883_308_795_349_485_9,
492        228 => 0.000_877_516_675_806_960_8,
493        229 => 0.000_871_787_581_542_502_7,
494        230 => 0.000_866_120_556_613_433,
495        231 => 0.000_860_514_663_638_586_4,
496        232 => 0.000_854_968_983_359_706_2,
497        233 => 0.000_849_482_614_215_034_8,
498        234 => 0.000_844_054_671_924_712_5,
499        235 => 0.000_838_684_289_087_606_4,
500        236 => 0.000_833_370_614_789_207_2,
501        237 => 0.000_828_112_814_220_247_1,
502        238 => 0.000_822_910_068_305_699,
503        239 => 0.000_817_761_573_343_835,
504        240 => 0.000_812_666_540_655_029_1,
505        241 => 0.000_807_624_196_240_001_8,
506        242 => 0.000_802_633_780_447_217_2,
507        243 => 0.000_797_694_547_649_146_8,
508        244 => 0.000_792_805_765_927_132_3,
509        245 => 0.000_787_966_716_764_583_6,
510        246 => 0.000_783_176_694_748_257,
511        247 => 0.000_778_435_007_277_372_5,
512        248 => 0.000_773_740_974_280_329_3,
513        249 => 0.000_769_093_927_938_796_6,
514        250 => 0.000_764_493_212_418_953_9,
515        251 => 0.000_759_938_183_609_671_9,
516        252 => 0.000_755_428_208_867_465_4,
517        253 => 0.000_750_962_666_767_783,
518        254 => 0.000_746_540_946_863_032_5,
519        255 => 0.000_742_162_449_446_367_6,
520    }
521}
conspire/physics/molecular/single_chain/fjc/mod.rs

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