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

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

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