From 9ba776b781db1ef8baece6aeb6b40a90797b877d Mon Sep 17 00:00:00 2001 From: mrbuche Date: Tue, 5 Dec 2023 18:16:30 -0700 Subject: [PATCH] MCMC FRC --- src/physics/single_chain/frc/mod.rs | 2 + .../frc/thermodynamics/isometric/mod.rs | 2 + .../isometric/monte_carlo/mod.rs | 84 +++++++++++++++++++ .../isometric/monte_carlo/test.rs | 43 ++++++++++ .../single_chain/frc/thermodynamics/mod.rs | 2 + src/physics/single_chain/mod.rs | 3 + tests/frc.rs | 26 ++++++ 7 files changed, 162 insertions(+) create mode 100644 src/physics/single_chain/frc/mod.rs create mode 100644 src/physics/single_chain/frc/thermodynamics/isometric/mod.rs create mode 100644 src/physics/single_chain/frc/thermodynamics/isometric/monte_carlo/mod.rs create mode 100644 src/physics/single_chain/frc/thermodynamics/isometric/monte_carlo/test.rs create mode 100644 src/physics/single_chain/frc/thermodynamics/mod.rs create mode 100644 tests/frc.rs diff --git a/src/physics/single_chain/frc/mod.rs b/src/physics/single_chain/frc/mod.rs new file mode 100644 index 00000000..a3b4212b --- /dev/null +++ b/src/physics/single_chain/frc/mod.rs @@ -0,0 +1,2 @@ +/// The freely-rotating chain (FRC) model thermodynamics. +pub mod thermodynamics; diff --git a/src/physics/single_chain/frc/thermodynamics/isometric/mod.rs b/src/physics/single_chain/frc/thermodynamics/isometric/mod.rs new file mode 100644 index 00000000..7565ab9a --- /dev/null +++ b/src/physics/single_chain/frc/thermodynamics/isometric/mod.rs @@ -0,0 +1,2 @@ +/// The freely-rotating chain (FRC) model thermodynamics in the isometric ensemble calculated using Monte Carlo methods. +pub mod monte_carlo; diff --git a/src/physics/single_chain/frc/thermodynamics/isometric/monte_carlo/mod.rs b/src/physics/single_chain/frc/thermodynamics/isometric/monte_carlo/mod.rs new file mode 100644 index 00000000..bd53e2eb --- /dev/null +++ b/src/physics/single_chain/frc/thermodynamics/isometric/monte_carlo/mod.rs @@ -0,0 +1,84 @@ +mod test; + +use rand::prelude::*; + +use std::f64::consts::PI; +use std::f64::consts::TAU as TWO_PI; + +pub fn cross(u: &[f64; 3], v: &[f64; 3]) -> [f64; 3] +{ + [u[1] * v[2] - u[2] * v[1], + u[2] * v[0] - u[0] * v[2], + u[0] * v[1] - u[1] * v[0]] +} + +pub fn random_configuration(theta: &f64, rng: &mut ThreadRng) -> [[f64; 3]; NUMBER_OF_LINKS] +{ + let mut phi: f64 = 0.0; + let mut phi_cos: f64 = 0.0; + let mut phi_sin: f64 = 0.0; + let theta_cos: f64 = theta.cos(); + let theta_sin: f64 = theta.sin(); + let mut configuration = [[0.0; 3]; NUMBER_OF_LINKS]; + let mut position = [0.0; 3]; + let mut r = [1.0, 0.0, 0.0]; + let mut t = [0.0; 3]; + let mut u = [0.0; 3]; + let mut v = [0.0; 3]; + let mut u_normalization = 0.0; + configuration.iter_mut().for_each(|coordinate|{ + phi = TWO_PI * rng.gen::(); + phi_cos = phi.cos(); + phi_sin = phi.sin(); + t = std::array::from_fn(|_| rng.gen::()); + u = cross(&r, &t); + u_normalization = u.iter().map(|u_i| u_i * u_i).sum::().sqrt(); + u.iter_mut().for_each(|u_i| *u_i /= u_normalization); + v = cross(&r, &u); + r.iter_mut().zip(u.iter().zip(v.iter())).for_each(|(r_i, (u_i, v_i))| + *r_i = (u_i * phi_cos + v_i * phi_sin) * theta_sin + *r_i * theta_cos + ); + position.iter_mut().zip(r.iter()).for_each(|(position_i, r_i)| + *position_i += r_i + ); + coordinate.iter_mut().zip(position.iter()).for_each(|(coordinate_i, position_i)| + *coordinate_i = *position_i + ); + }); + configuration +} + +pub fn random_nondimensional_end_to_end_length(theta: &f64, rng: &mut ThreadRng) -> f64 +{ + random_configuration::(theta, rng)[NUMBER_OF_LINKS - 1].iter().map(|entry| entry * entry).sum::().sqrt() +} + +pub fn nondimensional_equilibrium_radial_distribution(theta: &f64, number_of_samples: usize) -> ([f64; NUMBER_OF_BINS], [f64; NUMBER_OF_BINS]) +{ + let mut rng = rand::thread_rng(); + let number_of_links_f64 = NUMBER_OF_LINKS as f64; + let gamma_max = (2.0 - 2.0*(PI - theta).cos()).sqrt()/2.0; + let mut bin_centers = [0.0_f64; NUMBER_OF_BINS]; + bin_centers.iter_mut().enumerate().for_each(|(bin_index, bin_center)| + *bin_center = gamma_max*((bin_index as f64) + 0.5)/(NUMBER_OF_BINS as f64) + ); + let mut bin_counts = [0_u128; NUMBER_OF_BINS]; + let mut gamma: f64 = 0.0; + (0..number_of_samples).for_each(|_|{ + gamma = random_nondimensional_end_to_end_length::(theta, &mut rng)/number_of_links_f64; + for (bin_center, bin_count) in bin_centers.iter().zip(bin_counts.iter_mut()) + { + if &gamma < bin_center + { + *bin_count += 1; + break + } + } + }); + let normalization = (number_of_samples as f64)/(NUMBER_OF_BINS as f64); + let mut bin_probabilities = [0.0_f64; NUMBER_OF_BINS]; + bin_probabilities.iter_mut().zip(bin_counts.iter()).for_each(|(bin_probability, bin_count)| + *bin_probability = (*bin_count as f64)/normalization + ); + (bin_centers, bin_probabilities) +} \ No newline at end of file diff --git a/src/physics/single_chain/frc/thermodynamics/isometric/monte_carlo/test.rs b/src/physics/single_chain/frc/thermodynamics/isometric/monte_carlo/test.rs new file mode 100644 index 00000000..a24cf799 --- /dev/null +++ b/src/physics/single_chain/frc/thermodynamics/isometric/monte_carlo/test.rs @@ -0,0 +1,43 @@ +#![cfg(test)] + +use super::*; + +use std::f64::consts::PI; + +const NUMBER_OF_BINS: usize = 100; +const NUMBER_OF_LINKS: usize = 8; +const NUMBER_OF_SAMPLES: usize = 100000; +const THETA: f64 = PI/8.0; + +#[test] +fn monte_carlo_random_configuration() +{ + let mut rng = rand::thread_rng(); + let configuration = random_configuration::(&THETA, &mut rng); + assert_eq!(configuration.len(), NUMBER_OF_LINKS); + assert_eq!(configuration[0].len(), 3); +} + +#[test] +fn monte_carlo_random_nondimensional_end_to_end_length() +{ + let mut rng = rand::thread_rng(); + let gamma = random_nondimensional_end_to_end_length::(&THETA, &mut rng)/(NUMBER_OF_LINKS as f64); + let gamma_max = (2.0 - 2.0*(PI - THETA).cos()).sqrt()/2.0; + assert!(gamma >= 0.0); + assert!(gamma <= gamma_max); +} + +#[test] +fn monte_carlo_nondimensional_equilibrium_radial_distribution() +{ + let (gamma, g_eq) = nondimensional_equilibrium_radial_distribution::(&THETA, NUMBER_OF_SAMPLES); + assert_eq!(gamma.len(), NUMBER_OF_BINS); + assert_eq!(g_eq.len(), NUMBER_OF_BINS); + let gamma_max = (2.0 - 2.0*(PI - THETA).cos()).sqrt()/2.0; + gamma.iter().zip(g_eq.iter()).for_each(|(gamma_i, g_eq_i)|{ + assert!(gamma_i > &0.0); + assert!(gamma_i < &gamma_max); + assert!(g_eq_i >= &0.0); + }); +} diff --git a/src/physics/single_chain/frc/thermodynamics/mod.rs b/src/physics/single_chain/frc/thermodynamics/mod.rs new file mode 100644 index 00000000..f0025dcf --- /dev/null +++ b/src/physics/single_chain/frc/thermodynamics/mod.rs @@ -0,0 +1,2 @@ +/// The freely-rotating chain (FRC) model thermodynamics in the isometric ensemble. +pub mod isometric; diff --git a/src/physics/single_chain/mod.rs b/src/physics/single_chain/mod.rs index 1f122f84..3eaecc4c 100644 --- a/src/physics/single_chain/mod.rs +++ b/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 freely-rotating chain (FRC) single-chain model. +pub mod frc; + /// The worm-like chain (WLC) single-chain model. pub mod wlc; diff --git a/tests/frc.rs b/tests/frc.rs new file mode 100644 index 00000000..e60d677f --- /dev/null +++ b/tests/frc.rs @@ -0,0 +1,26 @@ +use polymers::physics::single_chain:: +{ + frc::thermodynamics::isometric::monte_carlo::nondimensional_equilibrium_radial_distribution, + wlc::thermodynamics::isometric::WLC +}; + +const KAPPA: f64 = 5.0/11.0; +const NUMBER_OF_BINS: usize = 1000; +const NUMBER_OF_LINKS: usize = 256; +const NUMBER_OF_SAMPLES: usize = 10000000; +const TOL: f64 = 5e-2; + +#[test] +fn monte_carlo_nondimensional_equilibrium_radial_distribution_wlc_limit() +{ + let theta = (2.0/(NUMBER_OF_LINKS as f64)/KAPPA).sqrt(); + let (gamma, g_eq) = nondimensional_equilibrium_radial_distribution::(&theta, NUMBER_OF_SAMPLES); + let mut check = 0.0; + let mut residual = 0.0; + let wlc = WLC::init(1, 1.0, 1.0, KAPPA); + gamma.iter().zip(g_eq.iter()).for_each(|(gamma_i, g_eq_i)|{ + check = wlc.nondimensional_equilibrium_radial_distribution(&gamma_i); + residual = (g_eq_i - &check).abs(); + assert!(residual < TOL || residual/check < TOL || g_eq_i < &TOL); + }); +}