Implement feature: adjustable maximum displacement in Monte-Carlo

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;

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<Box<MCMove>>,
-    /// Cumulative frequencies of the Monte-Carlo moves
+    moves: Vec<(Box<MCMove>, MoveCounter)>,
+    /// Cummulative frequencies of the Monte-Carlo moves
     frequencies: Vec<f64>,
+    /// 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<rand::Rng>,
@@ -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,35 @@ 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<MCMove>, frequency: f64, target_acceptance: f64) {
+        if self.initialized {
+            fatal_error!(
+                "Monte-Carlo simulation has already been initialized, \
+                we can not add new moves."
+            );
+        }
+        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,32 +157,114 @@ 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<f64>,
+}
+
+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<f64>) -> MoveCounter {
+        let mut counter = MoveCounter{
+            ncalled: 0,
+            naccepted: 0,
+            nattempted: 0,
+            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: Option<f64>) {
+        // 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<f64> {
+        // Check if there exists an target_acceptance
+        if let Some(ta) = self.target_acceptance {
+            // Capture division by zero
+            if self.nattempted == 0 { return None };
+            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)]
 mod tests {
-    use sim::mc::{MonteCarlo, MCMove};
+    use sim::mc::{MonteCarlo, MCMove, MoveCounter};
     use sim::Propagator;
     use sys::{System, EnergyCache};
     use rand::Rng;
@@ -165,6 +276,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<f64>) {}
     }
 
     #[test]
@@ -196,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));
+    }
 }

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<Rng>) -> 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<f64>);
 }
 
 /// Select a random molecule in the system using `rng` as random number

src/core/src/sim/mc/moves/rotate.rs

@@ -22,8 +22,10 @@ pub struct Rotate {
     newpos: Vec<Vector3D>,
     /// 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<f64>,
+    range: Range<f64>,
 }
 
 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<f64>) {
+        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