Merge pull request #179 from sandialabs/wlc-isometric

WLC isometric

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "polymers"
-version = "0.3.3"
+version = "0.3.4"
 edition = "2021"
 description = "Polymers Modeling Library"
 license = "BSD-3-Clause"

Project.toml

@@ -1,7 +1,7 @@
 name = "Polymers"
 uuid = "8aef037c-a721-4e8a-9d81-eb7093daef2c"
 authors = ["mrbuche <[email protected]>"]
-version = "0.3.3"
+version = "0.3.4"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"

docs/source/physics/single_chain.rst

@@ -9,6 +9,7 @@ Single-chain models for polymer physics
    EFJC <single_chain/efjc>
    SWFJC <single_chain/swfjc>
    uFJC <single_chain/ufjc>
+   WLC <single_chain/wlc>
 
 .. automodule:: polymers.physics.single_chain
    :members:

docs/source/physics/single_chain/wlc.rst

@@ -0,0 +1,15 @@
+WLC model
+=========
+
+.. toctree::
+   :maxdepth: 1
+
+   Thermodynamics <wlc/thermodynamics>
+
+.. autoclass:: polymers.physics.single_chain.fjc::WLC(number_of_links, link_length, hinge_mass, persistance_length)
+
+   .. autoattribute:: number_of_links
+   .. autoattribute:: link_length
+   .. autoattribute:: hinge_mass
+   .. autoattribute:: persistance_length
+   .. autoattribute:: thermodynamics

docs/source/physics/single_chain/wlc/thermodynamics.rst

@@ -0,0 +1,15 @@
+WLC model thermodynamics
+========================
+
+.. toctree::
+   :maxdepth: 1
+
+   Isometric <thermodynamics/isometric>
+
+.. autoclass:: polymers.physics.single_chain.wlc.thermodynamics::WLC(number_of_links, link_length, hinge_mass, persistance_length)
+
+   .. autoattribute:: number_of_links
+   .. autoattribute:: link_length
+   .. autoattribute:: hinge_mass
+   .. autoattribute:: persistance_length
+   .. autoattribute:: isometric

docs/source/physics/single_chain/wlc/thermodynamics/isometric.rst

@@ -0,0 +1,36 @@
+WLC model thermodynamics (isometric)
+====================================
+
+.. toctree::
+   :maxdepth: 1
+
+   Legendre <isometric/legendre>
+
+.. autoclass:: polymers.physics.single_chain.wlc.thermodynamics.isometric::WLC(number_of_links, link_length, hinge_mass, persistance_length)
+
+   .. autoattribute:: number_of_links
+   .. autoattribute:: link_length
+   .. autoattribute:: hinge_mass
+   .. autoattribute:: persistance_length
+   .. autoattribute:: legendre
+   .. automethod:: force(end_to_end_length, temperature)
+   .. automethod:: nondimensional_force(nondimensional_end_to_end_length_per_link)
+   .. automethod:: helmholtz_free_energy(end_to_end_length, temperature)
+   .. automethod:: helmholtz_free_energy_per_link(end_to_end_length, temperature)
+   .. automethod:: relative_helmholtz_free_energy(end_to_end_length, temperature)
+   .. automethod:: relative_helmholtz_free_energy_per_link(end_to_end_length, temperature)
+   .. automethod:: nondimensional_helmholtz_free_energy(nondimensional_end_to_end_length_per_link, temperature)
+   .. automethod:: nondimensional_helmholtz_free_energy_per_link(nondimensional_end_to_end_length_per_link, temperature)
+   .. automethod:: nondimensional_relative_helmholtz_free_energy(nondimensional_end_to_end_length_per_link)
+   .. automethod:: nondimensional_relative_helmholtz_free_energy_per_link(nondimensional_end_to_end_length_per_link)
+   .. automethod:: equilibrium_distribution(end_to_end_length)
+   .. automethod:: nondimensional_equilibrium_distribution(nondimensional_end_to_end_length_per_link)
+   .. automethod:: equilibrium_radial_distribution(end_to_end_length)
+   .. automethod:: nondimensional_equilibrium_radial_distribution(nondimensional_end_to_end_length_per_link)
+
+.. raw::
+ html
+
+   <hr>
+
+.. footbibliography::

docs/source/physics/single_chain/wlc/thermodynamics/isometric/legendre.rst

@@ -0,0 +1,17 @@
+WLC model thermodynamics (isometric/legendre)
+=============================================
+
+.. autoclass:: polymers.physics.single_chain.wlc.thermodynamics.isometric.legendre::WLC(number_of_links, link_length, hinge_mass, persistance_length)
+
+   .. autoattribute:: number_of_links
+   .. autoattribute:: link_length
+   .. autoattribute:: hinge_mass
+   .. autoattribute:: persistance_length
+   .. automethod:: gibbs_free_energy(end_to_end_length, temperature)
+   .. automethod:: gibbs_free_energy_per_link(end_to_end_length, temperature)
+   .. automethod:: relative_gibbs_free_energy(end_to_end_length, temperature)
+   .. automethod:: relative_gibbs_free_energy_per_link(end_to_end_length, temperature)
+   .. automethod:: nondimensional_gibbs_free_energy(nondimensional_end_to_end_length_per_link, temperature)
+   .. automethod:: nondimensional_gibbs_free_energy_per_link(nondimensional_end_to_end_length_per_link, temperature)
+   .. automethod:: nondimensional_relative_gibbs_free_energy(nondimensional_end_to_end_length_per_link)
+   .. automethod:: nondimensional_relative_gibbs_free_energy_per_link(nondimensional_end_to_end_length_per_link)

docs/src/physics/single_chain.md

@@ -5,6 +5,7 @@
   * [Extensible freely-jointed chain (EFJC) model](../efjc)
   * [Square-well freely-jointed chain (SWFJC) model](../swfjc)
   * [Arbitrary link potential freely-jointed chain (uFJC) model](../ufjc)
+  * [Worm-like chain (WLC) model](../wlc)
 
 ```@autodocs
 Modules = [Polymers.Physics.SingleChain]

docs/src/physics/single_chain/wlc.md

@@ -0,0 +1,7 @@
+# Worm-like chain (WLC) model
+
+  * [WLC model thermodynamics](../../thermodynamics)
+
+```@autodocs
+Modules = [Polymers.Physics.SingleChain.WLC]
+```

docs/src/physics/single_chain/wlc/thermodynamics.md

@@ -0,0 +1,7 @@
+# WLC model thermodynamics
+
+  * [WLC model thermodynamics (isometric)](../../../isometric)
+
+```@autodocs
+Modules = [Polymers.Physics.SingleChain.WLC.Thermodynamics]
+```

docs/src/physics/single_chain/wlc/thermodynamics/isometric.md

@@ -0,0 +1,7 @@
+# WLC model thermodynamics (isometric)
+
+  * [WLC model thermodynamics (isometric/legendre)](../../../../legendre)
+
+```@autodocs
+Modules = [Polymers.Physics.SingleChain.WLC.Thermodynamics.Isometric]
+```

docs/src/physics/single_chain/wlc/thermodynamics/isometric/legendre.md

@@ -0,0 +1,5 @@
+# WLC model thermodynamics (isometric/legendre)
+
+```@autodocs
+Modules = [Polymers.Physics.SingleChain.WLC.Thermodynamics.Isometric.Legendre]
+```

src/math/mod.rs

@@ -92,6 +92,157 @@ pub fn lambert_w(x: &f64) -> f64
     w
 }
 
+pub fn bessel_i(nu: &u8, x: &f64) -> f64
+{
+    if nu == &0
+    {
+        bessel_i0(x)
+    }
+    else if nu == &1
+    {
+        bessel_i1(x)
+    }
+    else
+    {
+        -1.0
+    }
+}
+
+fn bessel_i0(x: &f64) -> f64
+{
+    if x < &7.75
+    {
+        let coefficients = vec![
+            1.0,
+            2.499_999_999_999_999e-1,
+            2.777_777_777_777_822e-2,
+            1.736_111_111_110_237e-3,
+            6.944_444_444_533_525e-5,
+            1.929_012_345_132_199e-6,
+            3.936_759_911_025_107e-8,
+            6.151_186_727_044_392e-10,
+            7.594_070_020_589_734e-12,
+            7.593_897_933_698_363e-14,
+            6.277_677_736_362_926e-16,
+            4.347_097_041_532_722e-18,
+            2.634_177_426_901_091e-20,
+            1.139_430_377_448_228e-22,
+            9.079_269_200_856_248e-25
+        ];
+        let t = 0.25*x.powi(2);
+        1.0 + t*coefficients.iter().enumerate().map(|(i, c)| c*t.powi(i.try_into().unwrap())).sum::<f64>()
+    }
+    else if x < &500.0
+    {
+        let coefficients = vec![
+             3.989_422_804_014_25e-1,
+             4.986_778_506_049_619e-2,
+             2.805_062_339_283_126e-2,
+             2.922_112_251_660_478e-2,
+             4.442_072_994_936_595e-2,
+             1.309_705_746_058_567e-1,
+            -3.350_522_802_317_27,
+             2.330_257_115_835_147e2,
+            -1.133_663_506_971_723e4,
+             4.240_576_743_178_673e5,
+            -1.231_570_285_956_987e7,
+             2.802_319_381_552_675e8,
+            -5.018_839_997_137_779e9,
+             7.080_292_430_151_091e10,
+            -7.842_610_821_248_111e11,
+             6.768_257_378_540_965e12,
+            -4.490_348_496_961_38e13,
+             2.241_552_399_669_589e14,
+            -8.134_264_678_656_593e14,
+             2.023_910_973_916_877e15,
+            -3.086_757_152_953_708e15,
+             2.175_875_438_638_19e15
+        ];
+        x.exp()/x.sqrt()*coefficients.iter().enumerate().map(|(i, c)| c/x.powi(i.try_into().unwrap())).sum::<f64>()
+    }
+    else
+    {
+        let coefficients = vec![
+            3.989_422_804_014_329e-1,
+            4.986_778_504_914_345e-2,
+            2.805_063_089_165_061e-2,
+            2.921_790_968_539_151e-2,
+            4.533_712_087_625_794e-2
+        ];
+        let expf = (0.5*x).exp();
+        (expf/x.sqrt()*coefficients.iter().enumerate().map(|(i, c)| c/x.powi(i.try_into().unwrap())).sum::<f64>())*expf
+    }
+}
+
+fn bessel_i1(x: &f64) -> f64
+{
+    if x < &7.75
+    {
+        let coefficients = vec![
+            8.333_333_333_333_333e-2,
+            6.944_444_444_444_341e-3,
+            3.472_222_222_225_921e-4,
+            1.157_407_407_354_987e-5,
+            2.755_731_926_254_79e-7,
+            4.920_949_692_800_671e-9,
+            6.834_657_311_305_621e-11,
+            7.593_969_849_687_574e-13,
+            6.904_822_652_741_917e-15,
+            5.220_157_095_351_373e-17,
+            3.410_720_494_727_771e-19,
+            1.625_212_890_947_171e-21,
+            1.332_898_928_162_29e-23
+        ];
+        let t = 0.25*x.powi(2);
+        let more_coefficients = vec![
+            1.0,
+            0.5,
+            coefficients.iter().enumerate().map(|(i, c)| c*t.powi(i.try_into().unwrap())).sum::<f64>()
+        ];
+        0.5*x*more_coefficients.iter().enumerate().map(|(i, c)| c*t.powi(i.try_into().unwrap())).sum::<f64>()
+    }
+    else if x < &500.0
+    {
+        let coefficients = vec![
+             3.989_422_804_014_406e-1,
+            -1.496_033_551_613_111e-1,
+            -4.675_104_253_598_537e-2,
+            -4.090_895_951_581_637e-2,
+            -5.719_036_414_430_205e-2,
+            -1.528_189_554_374_492e-1,
+             3.458_284_470_977_172e0,
+            -2.426_181_371_595_021e2,
+             1.178_785_865_993_44e4,
+            -4.404_655_582_443_487e5,
+             1.277_677_779_341_446e7,
+            -2.903_390_398_236_656e8,
+             5.192_386_898_222_206e9,
+            -7.313_784_438_967_834e10,
+             8.087_824_484_994_859e11,
+            -6.967_602_516_005_787e12,
+             4.614_040_809_616_582e13,
+            -2.298_849_639_457_172e14,
+             8.325_554_073_334_618e14,
+            -2.067_285_045_778_906e15,
+             3.146_401_654_361_325e15,
+            -2.213_318_202_179_221e15
+        ];
+        x.exp()/x.sqrt()*coefficients.iter().enumerate().map(|(i, c)| c/x.powi(i.try_into().unwrap())).sum::<f64>()
+    }
+    else
+    {
+        let coefficients = vec![
+             3.989_422_804_014_314e-1,
+            -1.496_033_551_467_584e-1,
+            -4.675_105_322_571_775e-2,
+            -4.090_421_597_376_992e-2,
+            -5.843_630_344_778_927e-2
+        ];
+        let expf = (0.5*x).exp();
+        (expf/x.sqrt()*coefficients.iter().enumerate().map(|(i, c)| c/x.powi(i.try_into().unwrap())).sum::<f64>())*expf
+    }
+}
+
 pub fn erf(x: &f64) -> f64
 {
     1.0 - erfc(x)

src/physics/single_chain/fjc/thermodynamics/isometric/test.jl

@@ -1557,4 +1557,5 @@ end
               abs(residual_rel) <= 1.0 / sqrt(number_of_links)
     end
 end
+
 end

src/physics/single_chain/mod.jl

@@ -28,6 +28,8 @@ struct Parameters
     hinge_mass_scale::Float64
     link_length_reference::Float64
     link_length_scale::Float64
+    persistance_length_reference::Float64
+    persistance_length_scale::Float64
     number_of_links_minimum::UInt8
     number_of_links_maximum::UInt8
     link_stiffness_reference::Float64
@@ -67,9 +69,11 @@ parameters = Parameters(
     12e-1,
     8,
     1e0,
+    1e-1,
     1e0,
-    1e0,
-    1e0,
+    1e-2,
+    25e-1,
+    49e-1,
     0x08,
     0x19,
     5e5,
@@ -106,5 +110,6 @@ include("fjc/mod.jl")
 include("efjc/mod.jl")
 include("swfjc/mod.jl")
 include("ufjc/mod.jl")
+include("wlc/mod.jl")
 
 end

src/physics/single_chain/mod.rs

@@ -18,6 +18,9 @@ pub mod swfjc;
 /// The arbitrary link potential freely-jointed chain (uFJC) single-chain model.
 pub mod ufjc;
 
+/// The worm-like chain (WLC) single-chain model.
+pub mod wlc;
+
 static ONE: f64 = 1.0;
 static ZERO: f64 = 1e-6;
 static POINTS: u128 = 64;

src/physics/single_chain/py.rs

@@ -8,6 +8,7 @@ pub fn register_module(py: Python<'_>, parent_module: &PyModule) -> PyResult<()>
     super::efjc::py::register_module(py, single_chain)?;
     super::swfjc::py::register_module(py, single_chain)?;
     super::ufjc::py::register_module(py, single_chain)?;
+    super::wlc::py::register_module(py, single_chain)?;
     parent_module.add_submodule(single_chain)?;
     Ok(())
 }

src/physics/single_chain/test.py

@@ -28,9 +28,11 @@ def __init__(self):
         self.log_log_scale = 12e-1
         self.number_of_loops = 8
         self.hinge_mass_reference = 1e0
-        self.hinge_mass_scale = 1e0
+        self.hinge_mass_scale = 1e-1
         self.link_length_reference = 1e0
-        self.link_length_scale = 1e0
+        self.link_length_scale = 1e-2
+        self.persistance_length_reference = 25e-1
+        self.persistance_length_scale = 49e-1
         self.number_of_links_minimum = 5
         self.number_of_links_maximum = 25
         self.link_stiffness_reference = 5e5

src/physics/single_chain/test.rs

@@ -18,6 +18,8 @@ pub struct Parameters
     pub hinge_mass_scale: f64,
     pub link_length_reference: f64,
     pub link_length_scale: f64,
+    pub persistance_length_reference: f64,
+    pub persistance_length_scale: f64,
     pub number_of_links_minimum: u8,
     pub number_of_links_maximum: u8,
     pub link_stiffness_reference: f64,
@@ -61,9 +63,11 @@ impl Default for Parameters
             log_log_scale: 12e-1,
             number_of_loops: 8,
             hinge_mass_reference: 1e0,
-            hinge_mass_scale: 1e0,
+            hinge_mass_scale: 1e-1,
             link_length_reference: 1e0,
-            link_length_scale: 1e0,
+            link_length_scale: 1e-2,
+            persistance_length_reference: 25e-1,
+            persistance_length_scale: 49e-1,
             number_of_links_minimum: 5,
             number_of_links_maximum: 25,
             link_stiffness_reference: 5e5,