From a12dc90fb45a9410456ff2b72059d59b8eb1e564 Mon Sep 17 00:00:00 2001 From: mrbuche Date: Tue, 12 Dec 2023 22:55:18 -0700 Subject: [PATCH] seems almost perfect, why do they tend to seem translated to the right? --- Cargo.toml | 1 + .../efjc/thermodynamics/isometric/mod.rs | 3 + .../isometric/monte_carlo/mod.rs | 63 +++++++++++++++++++ .../isometric/monte_carlo/test.rs | 40 ++++++++++++ .../isometric/monte_carlo/mod.rs | 4 +- tests/efjc.rs | 16 +++++ tests/frc.rs | 2 +- 7 files changed, 126 insertions(+), 3 deletions(-) create mode 100644 src/physics/single_chain/efjc/thermodynamics/isometric/monte_carlo/mod.rs create mode 100644 src/physics/single_chain/efjc/thermodynamics/isometric/monte_carlo/test.rs create mode 100644 tests/efjc.rs diff --git a/Cargo.toml b/Cargo.toml index ec598d54..42c97252 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -17,4 +17,5 @@ python = ["dep:numpy", "dep:pyo3"] numpy = {version = "0.19", optional = true} pyo3 = {version = "0.19", features = ["extension-module"], optional = true} rand = "0.8.5" +rand_distr = "0.4.3" diff --git a/src/physics/single_chain/efjc/thermodynamics/isometric/mod.rs b/src/physics/single_chain/efjc/thermodynamics/isometric/mod.rs index b826a9b4..639530b0 100644 --- a/src/physics/single_chain/efjc/thermodynamics/isometric/mod.rs +++ b/src/physics/single_chain/efjc/thermodynamics/isometric/mod.rs @@ -6,6 +6,9 @@ mod test; /// The extensible freely-jointed chain (EFJC) model thermodynamics in the isometric ensemble approximated using an asymptotic approach. pub mod asymptotic; +/// The extensible freely-jointed chain (EFJC) model thermodynamics in the isometric ensemble calculated using Monte Carlo methods. +pub mod monte_carlo; + /// The structure of the thermodynamics of the EFJC model in the isometric ensemble. pub struct EFJC { diff --git a/src/physics/single_chain/efjc/thermodynamics/isometric/monte_carlo/mod.rs b/src/physics/single_chain/efjc/thermodynamics/isometric/monte_carlo/mod.rs new file mode 100644 index 00000000..bef697c3 --- /dev/null +++ b/src/physics/single_chain/efjc/thermodynamics/isometric/monte_carlo/mod.rs @@ -0,0 +1,63 @@ +mod test; + +use rand::prelude::*; + +use rand_distr::{Normal, Distribution}; + +use std::f64::consts::TAU as TWO_PI; + +pub fn random_configuration(dist: Normal, rng: &mut ThreadRng) -> [[f64; 3]; NUMBER_OF_LINKS] +{ + let mut lambda: f64 = 0.0; + let mut phi: f64 = 0.0; + let mut theta: f64 = 0.0; + let mut position = [0.0; 3]; + let mut configuration = [[0.0; 3]; NUMBER_OF_LINKS]; + configuration.iter_mut().for_each(|coordinate|{ + lambda = dist.sample(rng); + phi = TWO_PI * rng.gen::(); + theta = (1.0 - 2.0 * rng.gen::()).acos(); + position[0] += lambda * (theta.sin() * phi.cos()); + position[1] += lambda * (theta.sin() * phi.sin()); + position[2] += lambda * theta.cos(); + 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(dist: Normal, rng: &mut ThreadRng) -> f64 +{ + random_configuration::(dist, rng)[NUMBER_OF_LINKS - 1].iter().map(|entry| entry * entry).sum::().sqrt() +} + +pub fn nondimensional_equilibrium_radial_distribution(gamma_max: &f64, kappa: &f64, number_of_samples: usize) -> ([f64; NUMBER_OF_BINS], [f64; NUMBER_OF_BINS]) +{ + let mut rng = rand::thread_rng(); + let dist = Normal::new(1.0, 1.0/kappa.sqrt()).unwrap(); + let number_of_links_f64 = NUMBER_OF_LINKS as f64; + 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::(dist, &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 = gamma_max * (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) +} diff --git a/src/physics/single_chain/efjc/thermodynamics/isometric/monte_carlo/test.rs b/src/physics/single_chain/efjc/thermodynamics/isometric/monte_carlo/test.rs new file mode 100644 index 00000000..9854d4f5 --- /dev/null +++ b/src/physics/single_chain/efjc/thermodynamics/isometric/monte_carlo/test.rs @@ -0,0 +1,40 @@ +#![cfg(test)] + +use super::*; + +const GAMMA_MAX: f64 = 1.5; +const KAPPA: f64 = 50.0; +const NUMBER_OF_BINS: usize = 100; +const NUMBER_OF_LINKS: usize = 3; +const NUMBER_OF_SAMPLES: usize = 100000; + +#[test] +fn monte_carlo_random_configuration() +{ + let mut rng = rand::thread_rng(); + let dist = Normal::new(1.0, 1.0/KAPPA.sqrt()).unwrap(); + let configuration = random_configuration::(dist, &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 dist = Normal::new(1.0, 1.0/KAPPA.sqrt()).unwrap(); + let gamma = random_nondimensional_end_to_end_length::(dist, &mut rng)/(NUMBER_OF_LINKS as f64); + assert!(gamma >= 0.0); +} + +#[test] +fn monte_carlo_nondimensional_equilibrium_radial_distribution() +{ + let (gamma, g_eq) = nondimensional_equilibrium_radial_distribution::(&GAMMA_MAX, &KAPPA, NUMBER_OF_SAMPLES); + assert_eq!(gamma.len(), NUMBER_OF_BINS); + assert_eq!(g_eq.len(), NUMBER_OF_BINS); + gamma.iter().zip(g_eq.iter()).for_each(|(gamma_i, g_eq_i)|{ + assert!(gamma_i > &0.0); + assert!(g_eq_i >= &0.0); + }); +} \ No newline at end of file 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 index bd53e2eb..f3f3ae15 100644 --- a/src/physics/single_chain/frc/thermodynamics/isometric/monte_carlo/mod.rs +++ b/src/physics/single_chain/frc/thermodynamics/isometric/monte_carlo/mod.rs @@ -60,7 +60,7 @@ pub fn nondimensional_equilibrium_radial_distribution(&GAMMA_MAX, &KAPPA, NUMBER_OF_SAMPLES); + gamma.iter().zip(g_eq.iter()).for_each(|output|{ + println!("{:?}", output) + }); +} diff --git a/tests/frc.rs b/tests/frc.rs index e60d677f..a6b6eea3 100644 --- a/tests/frc.rs +++ b/tests/frc.rs @@ -8,7 +8,7 @@ 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; +const TOL: f64 = 8e-2; #[test] fn monte_carlo_nondimensional_equilibrium_radial_distribution_wlc_limit()