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

1use crate::math::Current;
2use crate::{
3    math::{Quantity, Scalar, random::random_uniform},
4    mechanics::Vector,
5    physics::molecular::{
6        potential::{Harmonic, Potential},
7        single_chain::{
8            Configuration, Ensemble, Isometric, Isotensional, Legendre, MonteCarlo, SingleChain,
9            SingleChainError, Thermodynamics,
10        },
11    },
12    units::Length,
13};
14use std::f64::consts::TAU;
15
16/// Options for arbitrary discrete potentials.
17#[derive(Clone, Debug)]
18pub enum ArbitraryDiscretePotential<U>
19where
20    U: Potential,
21{
22    Free,
23    Rigid(Scalar),
24    Strong(U),
25    Weak(U),
26}
27
28/// The arbitrary discrete single-chain model.
29#[derive(Clone, Debug)]
30pub struct ArbitraryDiscrete {
31    /// The number of links $`N_b`$.
32    pub number_of_links: u8,
33    /// The link potential $`u_b`$.
34    pub link_potential: ArbitraryDiscretePotential<Harmonic>,
35    /// The angular potential $`u_\theta`$.
36    pub angular_potential: ArbitraryDiscretePotential<Harmonic>,
37    /// The torsional potential $`u_\phi`$.
38    pub torsional_potential: ArbitraryDiscretePotential<Harmonic>,
39    /// The thermodynamic ensemble.
40    pub ensemble: Ensemble,
41}
42
43impl SingleChain for ArbitraryDiscrete {
44    fn link_length(&self) -> Quantity<Length> {
45        match &self.link_potential {
46            ArbitraryDiscretePotential::Free => panic!(),
47            ArbitraryDiscretePotential::Rigid(link_length) => Quantity::new(*link_length),
48            ArbitraryDiscretePotential::Strong(link_potential) => link_potential.rest_length(),
49            ArbitraryDiscretePotential::Weak(link_potential) => link_potential.rest_length(),
50        }
51    }
52    fn number_of_links(&self) -> u8 {
53        self.number_of_links
54    }
55}
56
57impl Thermodynamics for ArbitraryDiscrete {
58    fn ensemble(&self) -> Ensemble {
59        self.ensemble
60    }
61}
62
63impl Isometric for ArbitraryDiscrete {
64    fn nondimensional_helmholtz_free_energy(
65        &self,
66        _nondimensional_force: Scalar,
67    ) -> Result<Scalar, SingleChainError> {
68        unimplemented!()
69    }
70    fn nondimensional_force(
71        &self,
72        _nondimensional_force: Scalar,
73    ) -> Result<Scalar, SingleChainError> {
74        unimplemented!()
75    }
76    fn nondimensional_stiffness(
77        &self,
78        _nondimensional_force: Scalar,
79    ) -> Result<Scalar, SingleChainError> {
80        unimplemented!()
81    }
82    fn nondimensional_spherical_distribution(
83        &self,
84        _nondimensional_force: Scalar,
85    ) -> Result<Scalar, SingleChainError> {
86        unimplemented!()
87    }
88}
89
90impl Isotensional for ArbitraryDiscrete {
91    fn nondimensional_gibbs_free_energy_per_link(
92        &self,
93        _nondimensional_force: Scalar,
94    ) -> Result<Scalar, SingleChainError> {
95        unimplemented!()
96    }
97    fn nondimensional_extension(
98        &self,
99        _nondimensional_force: Scalar,
100    ) -> Result<Scalar, SingleChainError> {
101        unimplemented!()
102    }
103    fn nondimensional_compliance(
104        &self,
105        _nondimensional_force: Scalar,
106    ) -> Result<Scalar, SingleChainError> {
107        unimplemented!()
108    }
109}
110
111impl Legendre for ArbitraryDiscrete {}
112
113impl MonteCarlo for ArbitraryDiscrete {
114    fn random_nondimensional_link_vectors(&self, nondimensional_force: Scalar) -> Configuration {
115        //
116        // Need to add cases and get them right.
117        //
118        let eta = nondimensional_force;
119        let eta_exp = eta.exp();
120        let eta_nexp = 1.0 / eta_exp;
121        (0..self.number_of_links())
122            .map(|_| {
123                let cos_theta = if eta == 0.0 {
124                    2.0 * random_uniform() - 1.0
125                } else {
126                    (eta_nexp + random_uniform() * (eta_exp - eta_nexp)).ln() / eta
127                };
128                let sin_theta = (1.0 - cos_theta * cos_theta).sqrt();
129                let phi = TAU * random_uniform();
130                let (sin_phi, cos_phi) = phi.sin_cos();
131                Vector::<Current>::from([sin_theta * cos_phi, sin_theta * sin_phi, cos_theta])
132            })
133            .collect()
134    }
135}
conspire/physics/molecular/single_chain/arbitrary/mod.rs

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