diff --git a/src/meas.rs b/src/meas.rs
index 3e51acc..b87ac13 100644
--- a/src/meas.rs
+++ b/src/meas.rs
@@ -50,7 +50,11 @@ impl AbstractMeasurement for Label {
     }
 }
 
-fn sum_over_region<T: Default + Sum<T>, F: Fn(Meters, &Cell) -> T>(state: &dyn GenericSim, r: &dyn Region, f: F) -> T {
+fn sum_over_region<T, F>(state: &dyn GenericSim, r: &dyn Region, f: F) -> T
+where 
+    T: Sum<T> + Default + Send,
+    F: Fn(Meters, &Cell) -> T + Sync,
+{
     state.map_sum_enumerated(|coord, cell| {
         if r.contains(coord) {
             f(coord, cell)
diff --git a/src/sim.rs b/src/sim.rs
index d887ec6..2fc44bb 100644
--- a/src/sim.rs
+++ b/src/sim.rs
@@ -6,6 +6,7 @@ use dyn_clone::{self, DynClone};
 use log::trace;
 
 use ndarray::{Array3, Zip};
+use rayon::prelude::*;
 use serde::{Serialize, Deserialize};
 use std::convert::From;
 use std::iter::Sum;
@@ -56,22 +57,26 @@ impl<'a> dyn GenericSim + 'a {
     /// Apply `F` to each Cell, and sum the results.
     pub fn map_sum<F, Ret>(&self, f: F) -> Ret
     where
-          F: Fn(&Cell<mat::Static>) -> Ret,
-          Ret: Sum<Ret>,
+          F: Fn(&Cell<mat::Static>) -> Ret + Sync,
+          Ret: Sum<Ret> + Send,
     {
         self.map_sum_enumerated(|_at: Index, cell| f(cell))
     }
 
     pub fn map_sum_enumerated<C, F, Ret>(&self, f: F) -> Ret
     where C: Coord,
-          F: Fn(C, &Cell<mat::Static>) -> Ret,
-          Ret: Sum<Ret>
+          F: Fn(C, &Cell<mat::Static>) -> Ret + Sync,
+          Ret: Sum<Ret> + Send,
     {
         let (w, h, d) = (self.width(), self.height(), self.depth());
-        (0..d).map(|z| (0..h).map(|y| (0..w).map(|x| {
+        (0..d).into_par_iter().map(
+            |z| (0..h).into_par_iter().map_with(z,
+            |&mut z, y| (0..w).into_par_iter().map_with((z, y),
+            |&mut (z, y), x|
+        {
             let at = Index(Vec3u::new(x, y, z));
             f(C::from_index(at, self.feature_size()), &self.get(at))
-        }).sum()).sum()).sum()
+        }))).flatten().flatten().sum()
     }
 
     pub fn volume_of_region<R: Region>(&self, region: &R) -> u32 {
@@ -81,8 +86,8 @@ impl<'a> dyn GenericSim + 'a {
     /// Apply `F` to each Cell, and sum the results.
     pub fn map_sum_over<F, Ret, Reg>(&self, region: &Reg, f: F) -> Ret
     where
-          F: Fn(&Cell<mat::Static>) -> Ret,
-          Ret: Sum<Ret> + Default,
+          F: Fn(&Cell<mat::Static>) -> Ret + Sync,
+          Ret: Sum<Ret> + Default + Send,
           Reg: Region
     {
         self.map_sum_over_enumerated(region, |_at: Index, cell| f(cell))
@@ -90,12 +95,16 @@ impl<'a> dyn GenericSim + 'a {
 
     pub fn map_sum_over_enumerated<C, F, Ret, Reg>(&self, region: &Reg, f: F) -> Ret
     where C: Coord,
-          F: Fn(C, &Cell<mat::Static>) -> Ret,
-          Ret: Sum<Ret> + Default,
+          F: Fn(C, &Cell<mat::Static>) -> Ret + Sync,
+          Ret: Sum<Ret> + Default + Send,
           Reg: Region,
     {
         let (w, h, d) = (self.width(), self.height(), self.depth());
-        (0..d).map(|z| (0..h).map(|y| (0..w).map(|x| {
+        (0..d).into_par_iter().map(
+            |z| (0..h).into_par_iter().map_with(z,
+            |&mut z, y| (0..w).into_par_iter().map_with((z, y),
+            |&mut (z, y), x|
+        {
             let at = Index(Vec3u::new(x, y, z));
             let meters = at.to_meters(self.feature_size());
             if region.contains(meters) {
@@ -103,7 +112,7 @@ impl<'a> dyn GenericSim + 'a {
             } else {
                 Default::default()
             }
-        }).sum()).sum()).sum()
+        }))).flatten().flatten().sum()
     }
 
     pub fn current<C: Coord>(&self, c: C) -> Vec3 {