From f5d8e527c883d7b9afda35ebebd7724b514bd93a Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Tue, 1 Nov 2016 20:50:00 +0100 Subject: [PATCH 1/3] Implementation of adjustments of MC moves ranges. --- src/core/src/sim/mc/mod.rs | 2 +- src/core/src/sim/mc/monte_carlo.rs | 124 +++++++++++++++++++++++-- src/core/src/sim/mc/moves/mod.rs | 20 ++-- src/core/src/sim/mc/moves/rotate.rs | 19 +++- src/core/src/sim/mc/moves/translate.rs | 39 ++++---- 5 files changed, 168 insertions(+), 36 deletions(-) diff --git a/src/core/src/sim/mc/mod.rs b/src/core/src/sim/mc/mod.rs index 575ac2488..9ea968b07 100644 --- a/src/core/src/sim/mc/mod.rs +++ b/src/core/src/sim/mc/mod.rs @@ -3,7 +3,7 @@ //! Monte-Carlo Metropolis algorithms mod monte_carlo; -pub use self::monte_carlo::MonteCarlo; +pub use self::monte_carlo::{MonteCarlo, MoveCounter}; mod moves; pub use self::moves::MCMove; diff --git a/src/core/src/sim/mc/monte_carlo.rs b/src/core/src/sim/mc/monte_carlo.rs index 5b9773ca1..ff22d2096 100644 --- a/src/core/src/sim/mc/monte_carlo.rs +++ b/src/core/src/sim/mc/monte_carlo.rs @@ -15,9 +15,12 @@ pub struct MonteCarlo { /// Boltzmann factor: beta = 1/(kB * T) beta: f64, /// List of possible Monte-Carlo moves - moves: Vec>, - /// Cumulative frequencies of the Monte-Carlo moves + moves: Vec<(Box, MoveCounter)>, + /// Cummulative frequencies of the Monte-Carlo moves frequencies: Vec, + /// Specifies the number of moves after which an update of a move's + /// amplitude is performed. + update_frequency: u64, /// Random number generator for the simulation. All random state will be /// taken from this. rng: Box, @@ -44,13 +47,14 @@ impl MonteCarlo { beta: 1.0 / (K_BOLTZMANN * temperature), moves: Vec::new(), frequencies: Vec::new(), + update_frequency: 0, rng: rng, cache: EnergyCache::new(), initialized: false, } } - /// Add a the `mcmove` Monte-Carlo move to this propagator, with frequency + /// Add the `mcmove` Monte-Carlo move to this propagator, with frequency /// `frequency`. All calls to this function should happen before any /// simulation run. /// @@ -64,10 +68,41 @@ impl MonteCarlo { we can not add new moves." ); } - self.moves.push(mcmove); + self.moves.push((mcmove, MoveCounter::new(None))); self.frequencies.push(frequency); } + /// Add the `mcmove` Monte-Carlo move to the propagator. + /// `frequency` describes how frequent a move is called, `target_acceptance` + /// is the desired acceptance ratio of the move. + /// + /// # Panics + /// + /// If called after a simulation run. + /// If `target_acceptance` is either negative or larger than one. + pub fn add_move_with_acceptance(&mut self, mcmove: Box, frequency: f64, target_acceptance: f64) { + if self.initialized { + fatal_error!( + "Monte-Carlo simulation has already been initialized, \ + we can not add new moves." + ); + } + if target_acceptance >= 1.0 || target_acceptance < 0.0 { + fatal_error!( + "The target acceptance ratio of the move has to be a positive \ + value smaller equal than 1." + ); + } + self.moves.push((mcmove, MoveCounter::new(Some(target_acceptance)))); + self.frequencies.push(frequency); + } + + /// Set the number of times a move has to be called before its amplitude + /// is updated. This value is applied to all moves. + pub fn set_amplitude_update_frequency(&mut self, frequency: u64) { + self.update_frequency = frequency; + } + /// Get the temperature of the simulation pub fn temperature(&self) -> f64 { 1.0 / (self.beta * K_BOLTZMANN) @@ -128,27 +163,98 @@ impl Propagator for MonteCarlo { .expect("Could not find a move in MonteCarlo moves list"); &mut self.moves[i] }; - trace!("Selected move is '{}'", mcmove.describe()); + trace!("Selected move is '{}'", mcmove.0.describe()); - if !mcmove.prepare(system, &mut self.rng) { + if !mcmove.0.prepare(system, &mut self.rng) { trace!(" --> Can not perform the move"); return; } - let cost = mcmove.cost(system, self.beta, &mut self.cache); + // attempt the move: increase counter + mcmove.1.ncalled += 1; + mcmove.1.nattempted += 1; + + // compute cost + let cost = mcmove.0.cost(system, self.beta, &mut self.cache); trace!(" --> Move cost is {}", cost); + // apply metropolis criterion let accepted = cost <= 0.0 || self.rng.next_f64() < f64::exp(-cost); if accepted { trace!(" --> Move was accepted"); - mcmove.apply(system); + mcmove.0.apply(system); self.cache.update(system); + mcmove.1.naccepted += 1 } else { trace!(" --> Move was rejected"); - mcmove.restore(system); + mcmove.0.restore(system); + } + + // Do the adjustments for the selected move as needed + if mcmove.1.nattempted == self.update_frequency { + mcmove.0.update_amplitude(mcmove.1.compute_scaling_factor()); + // Set `nattempted` and `naccepted` to zero. + // This way, only the running acceptance since the last update is considered. + mcmove.1.naccepted = 0; + mcmove.1.nattempted = 0; + } + } +} + + +/// This struct keeps track of the number of times a move was called +/// and how often it was accepted. +pub struct MoveCounter { + /// Count the total number of times the move was called. + pub ncalled: u64, + /// Count the number of times the move was accepted since the last update. + pub naccepted: u64, + /// Count the number of times the move was called since the last update. + pub nattempted: u64, + /// The target fraction of accepted over attempted moves. + /// The acceptance can be used to increase the efficiency of a move set. + target_acceptance: Option, +} + +impl MoveCounter { + /// Create a new counter for the move, initializing all counts to zero and + /// setting the `target_acceptance`. + pub fn new(target_acceptance: Option) -> MoveCounter { + MoveCounter{ + ncalled: 0, + naccepted: 0, + nattempted: 0, + target_acceptance: target_acceptance, } } + + /// Set the target acceptance for the move counter. + pub fn set_acceptance(&mut self, target_acceptance: f64) { + self.target_acceptance = Some(target_acceptance); + } + + /// Compute a scaling factor according to the desired acceptance. + pub fn compute_scaling_factor(&self) -> Option { + // Check if there exists an target_acceptance + if let Some(ta) = self.target_acceptance { + // Capture division by zero + if self.nattempted == 0 { return Some(1.0) }; + let quotient = self.naccepted as f64 / self.nattempted as f64 / ta; + // Limit the change + match quotient { + _ if quotient > 1.2 => Some(1.2), + _ if quotient < 0.8 => Some(0.8), + _ => Some(quotient), + } + } else { + None + } + } +} + +impl Default for MoveCounter { + fn default() -> MoveCounter { MoveCounter::new(None) } } #[cfg(test)] diff --git a/src/core/src/sim/mc/moves/mod.rs b/src/core/src/sim/mc/moves/mod.rs index 336372479..15c339eea 100644 --- a/src/core/src/sim/mc/moves/mod.rs +++ b/src/core/src/sim/mc/moves/mod.rs @@ -18,7 +18,7 @@ pub trait MCMove { /// Give a short description of this move fn describe(&self) -> &str; - /// Prepare the move, by selecting the particles to move, and the parameters + /// Prepare the move by selecting the particles to move, and the parameters /// of the move. The `rng` random number generator should be used to /// generate the parameters of the move. /// @@ -27,19 +27,27 @@ pub trait MCMove { fn prepare(&mut self, system: &mut System, rng: &mut Box) -> bool; /// Get the cost of performing this move on `system`. For example in - /// simple NVT simulations, this cost is the energetic difference over - /// `beta`. The cost must be dimensionless, and will be placed in an - /// exponential. The `cache` should be used to compute the cost, or the + /// simple NVT simulations, this cost is the energetic difference between + /// the new and the old state times beta. The cost must be dimmensionless. + /// + /// Note that the cost will be placed in an exponential with a negative sign. + /// For NVT using the Metropolis criterion: + /// cost = beta*(U_new - U_old) -> P_acc = min[1, exp(-cost)]. + /// + /// The `cache` should be used to compute the cost, or the /// `cache.unused` function should be used to ensure that the cache is - /// updated as needed after this move. + /// updated as needed after this move. fn cost(&self, system: &System, beta: f64, cache: &mut EnergyCache) -> f64; /// Apply the move, if it has not already been done in `prepare`. fn apply(&mut self, system: &mut System); - /// Restore the system to it's initial state, if it has been changed in + /// Restore the system to it's initial state if it has been changed in /// `prepare`. fn restore(&mut self, system: &mut System); + + /// Update the sample range for displacements. + fn update_amplitude(&mut self, scaling_factor: Option); } /// Select a random molecule in the system using `rng` as random number diff --git a/src/core/src/sim/mc/moves/rotate.rs b/src/core/src/sim/mc/moves/rotate.rs index 36082a887..432a399f5 100644 --- a/src/core/src/sim/mc/moves/rotate.rs +++ b/src/core/src/sim/mc/moves/rotate.rs @@ -22,8 +22,10 @@ pub struct Rotate { newpos: Vec, /// Normal distribution, for generation of the axis axis_rng: Normal, + /// Maximum values for the range of the range distribution of the angle + theta: f64, /// Range distribution, for generation of the angle - angle_rng: Range, + range: Range, } impl Rotate { @@ -47,7 +49,8 @@ impl Rotate { molid: usize::MAX, newpos: Vec::new(), axis_rng: Normal::new(0.0, 1.0), - angle_rng: Range::new(-theta, theta), + theta: theta, + range: Range::new(-theta, theta), } } } @@ -78,10 +81,11 @@ impl MCMove for Rotate { self.axis_rng.sample(rng), self.axis_rng.sample(rng) ).normalized(); - let theta = self.angle_rng.sample(rng); - + let theta = self.range.sample(rng); self.newpos.clear(); + let molecule = system.molecule(self.molid); + let mut masses = vec![0.0; molecule.size()]; for (i, pi) in molecule.iter().enumerate() { masses[i] = system[pi].mass; @@ -107,6 +111,13 @@ impl MCMove for Rotate { fn restore(&mut self, _: &mut System) { // Nothing to do } + + fn update_amplitude(&mut self, scaling_factor: Option) { + if let Some(s) = scaling_factor { + self.theta *= s; + self.range = Range::new(-self.theta, self.theta); + } + } } /// Rotate the particles at `positions` with the masses in `masses` around the diff --git a/src/core/src/sim/mc/moves/translate.rs b/src/core/src/sim/mc/moves/translate.rs index 92bc7618e..9dceb27b4 100644 --- a/src/core/src/sim/mc/moves/translate.rs +++ b/src/core/src/sim/mc/moves/translate.rs @@ -9,7 +9,7 @@ use std::usize; use super::MCMove; use super::select_molecule; -use types::{Vector3D, Zero}; +use types::Vector3D; use sys::{System, EnergyCache}; /// Monte-Carlo move for translating a molecule @@ -20,21 +20,21 @@ pub struct Translate { molid: usize, /// New positions of the atom in the translated molecule newpos: Vec, - /// Translation vector - delta: Vector3D, + /// Maximum displacement value + dr: f64, /// Translation range for random number generation - delta_range: Range, + range: Range, } impl Translate { - /// Create a new `Translate` move, with maximum displacement of `dr`, - /// translating all the molecules in the system. + /// Create a new `Translate` move, with maximum displacement of `dr`. + /// Translating all the molecules in the system. pub fn new(dr: f64) -> Translate { Translate::create(dr, None) } - /// Create a new `Translate` move, with maximum displacement of `dr`, - /// translating only molecules with `moltype` type. + /// Create a new `Translate` move, with maximum displacement of `dr`. + /// Translating only molecules with `moltype` type. pub fn with_moltype(dr: f64, moltype: u64) -> Translate { Translate::create(dr, Some(moltype)) } @@ -47,8 +47,8 @@ impl Translate { moltype: moltype, molid: usize::MAX, newpos: Vec::new(), - delta: Vector3D::zero(), - delta_range: Range::new(-dr, dr), + dr: dr, + range: Range::new(-dr, dr), } } } @@ -72,15 +72,15 @@ impl MCMove for Translate { return false; } - self.delta = Vector3D::new( - self.delta_range.sample(rng), - self.delta_range.sample(rng), - self.delta_range.sample(rng) + let delta = Vector3D::new( + self.range.sample(rng), + self.range.sample(rng), + self.range.sample(rng) ); self.newpos.clear(); for i in system.molecule(self.molid) { - self.newpos.push(system[i].position + self.delta); + self.newpos.push(system[i].position + delta); } return true; } @@ -98,6 +98,13 @@ impl MCMove for Translate { } fn restore(&mut self, _: &mut System) { - // Nothing to do + // Nothing to do. + } + + fn update_amplitude(&mut self, scaling_factor: Option) { + if let Some(s) = scaling_factor { + self.dr *= s; + self.range = Range::new(-self.dr, self.dr); + } } } From 3cc219c17bdd0eb72d6fa3b2e6fa042829824b9f Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Sat, 26 Nov 2016 13:09:55 +0100 Subject: [PATCH 2/3] Added missing implementation of update_amplitude in testcase of monte_carlo.rs --- src/core/src/sim/mc/monte_carlo.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/src/core/src/sim/mc/monte_carlo.rs b/src/core/src/sim/mc/monte_carlo.rs index ff22d2096..7ffdb6bf3 100644 --- a/src/core/src/sim/mc/monte_carlo.rs +++ b/src/core/src/sim/mc/monte_carlo.rs @@ -271,6 +271,7 @@ mod tests { fn cost(&self, _: &System, _: f64, _: &mut EnergyCache) -> f64 {0.0} fn apply(&mut self, _: &mut System) {} fn restore(&mut self, _: &mut System) {} + fn update_amplitude(&mut self, _:Option) {} } #[test] From a951e2296e537f936e0d8c314877e5a008e629f4 Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Sun, 27 Nov 2016 16:33:08 +0100 Subject: [PATCH 3/3] Moved error handling for invalid acceptances to set_acceptance. Added tests. --- src/core/src/sim/mc/monte_carlo.rs | 61 +++++++++++++++++++++++------- 1 file changed, 47 insertions(+), 14 deletions(-) diff --git a/src/core/src/sim/mc/monte_carlo.rs b/src/core/src/sim/mc/monte_carlo.rs index 7ffdb6bf3..4af39d1bb 100644 --- a/src/core/src/sim/mc/monte_carlo.rs +++ b/src/core/src/sim/mc/monte_carlo.rs @@ -87,12 +87,6 @@ impl MonteCarlo { we can not add new moves." ); } - if target_acceptance >= 1.0 || target_acceptance < 0.0 { - fatal_error!( - "The target acceptance ratio of the move has to be a positive \ - value smaller equal than 1." - ); - } self.moves.push((mcmove, MoveCounter::new(Some(target_acceptance)))); self.frequencies.push(frequency); } @@ -221,25 +215,36 @@ impl MoveCounter { /// Create a new counter for the move, initializing all counts to zero and /// setting the `target_acceptance`. pub fn new(target_acceptance: Option) -> MoveCounter { - MoveCounter{ + let mut counter = MoveCounter{ ncalled: 0, naccepted: 0, nattempted: 0, - target_acceptance: target_acceptance, - } + target_acceptance: None, + }; + counter.set_acceptance(target_acceptance); + counter } /// Set the target acceptance for the move counter. - pub fn set_acceptance(&mut self, target_acceptance: f64) { - self.target_acceptance = Some(target_acceptance); + pub fn set_acceptance(&mut self, target_acceptance: Option) { + // Check if `target_acceptance` has a valid value. + if target_acceptance.is_some() { + let ta = target_acceptance.unwrap(); + if ta >= 1.0 || ta < 0.0 { + fatal_error!( + "The target acceptance ratio of the move has to be a positive value smaller equal than 1." + ) + } + } + self.target_acceptance = target_acceptance; } - + /// Compute a scaling factor according to the desired acceptance. pub fn compute_scaling_factor(&self) -> Option { // Check if there exists an target_acceptance if let Some(ta) = self.target_acceptance { // Capture division by zero - if self.nattempted == 0 { return Some(1.0) }; + if self.nattempted == 0 { return None }; let quotient = self.naccepted as f64 / self.nattempted as f64 / ta; // Limit the change match quotient { @@ -259,7 +264,7 @@ impl Default for MoveCounter { #[cfg(test)] mod tests { - use sim::mc::{MonteCarlo, MCMove}; + use sim::mc::{MonteCarlo, MCMove, MoveCounter}; use sim::Propagator; use sys::{System, EnergyCache}; use rand::Rng; @@ -303,4 +308,32 @@ mod tests { mc.setup(&System::new()); mc.add(Box::new(DummyMove), 1.0); } + + #[test] + #[should_panic] + fn invalid_acceptance() { + let mut mc = MonteCarlo::new(100.0); + mc.add_move_with_acceptance(Box::new(DummyMove), 1.0, 0.5); + mc.moves[0].1.set_acceptance(Some(1.1)); + } + + #[test] + fn valid_acceptance() { + let mut mc = MonteCarlo::new(100.0); + mc.add_move_with_acceptance(Box::new(DummyMove), 1.0, 0.5); + assert_eq!(mc.moves[0].1.target_acceptance, Some(0.5)); + mc.moves[0].1.set_acceptance(None); + assert_eq!(mc.moves[0].1.target_acceptance, None); + } + + #[test] + fn scaling_factor() { + let mut counter = MoveCounter::new(Some(0.5)); + assert_eq!(counter.compute_scaling_factor(), None); + counter.nattempted = 10; + counter.naccepted = 10; + assert_eq!(counter.compute_scaling_factor(), Some(1.2)); + counter.naccepted = 0; + assert_eq!(counter.compute_scaling_factor(), Some(0.8)); + } }