diff --git a/examples/make_input_file_and_run.py b/examples/make_input_file_and_run.py
index 37e3a16..3c900f9 100644
--- a/examples/make_input_file_and_run.py
+++ b/examples/make_input_file_and_run.py
@@ -255,16 +255,22 @@
     'geometry_input': geometry_0D
 }
 
+# Test rustbca_py and 'output_dir' option
 input_data['options']['name'] = 'rustbca_input_file'
+if not os.path.exists(r'./outputs'):
+    os.makedirs('outputs')
+input_data['options']['output_dir'] = 'outputs'
 rustbca_py(input_data, mode)
-s = np.genfromtxt('rustbca_input_filesputtered.output', delimiter=',')
+s = np.genfromtxt('outputs/rustbca_input_filesputtered.output', delimiter=',')
 
 arrays = rustbca_local_py(input_data, mode)
 sputtered = arrays['sputtered']
 
 np.testing.assert_approx_equal(s[0, 2], np.array(arrays['energy'])[sputtered][0])
 
+# reset these options before running from command line
 input_data['options']['name'] = 'input_file'
+input_data['options']['output_dir'] = '.'
 
 # Attempt to cleanup line endings
 input_string = dumps(input_data).replace('\r', '')
diff --git a/examples/test_morse.py b/examples/test_morse.py
index 227203c..0e8e86c 100644
--- a/examples/test_morse.py
+++ b/examples/test_morse.py
@@ -87,7 +87,7 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True):
     mean_free_path_model = "LIQUID"
     interaction_potential = [[{{"KRC_MORSE"={{D=5.4971E-20, r0=2.782E-10, alpha=1.4198E10, k=7E10, x0=0.75E-10}}}}]]
     scattering_integral = [["GAUSS_LEGENDRE"]]
-    root_finder = [[{{"CPR"={{n0=2, nmax=100, epsilon=1E-9, complex_threshold=1E-3, truncation_threshold=1E-9, far_from_zero=1E9, interval_limit=1E-12, derivative_free=true}}}}]]
+    root_finder = [[{{"CPR"={{n0=2, nmax=100, epsilon=1E-3, complex_threshold=1E-3, truncation_threshold=1E-9, far_from_zero=1E9, interval_limit=1E-12, derivative_free=true}}}}]]
     num_threads = 4
     num_chunks = 10
 
@@ -172,7 +172,7 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True):
 #Running and plotting the H-Ni simulations with the Morse potential and updated Es
 num_energies = 15
 energies = np.logspace(-1, 4, num_energies)
-run_sim = True
+run_sim = False
 num_samples = 10000
 R_N = np.zeros(num_energies)
 R_E = np.zeros(num_energies)
@@ -180,8 +180,8 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True):
 R_E_2 = np.zeros(num_energies)
 
 for index, energy in enumerate(energies):
-    R_N[index], R_E[index] = run_krc_morse_potential(energy, index, num_samples=num_samples, run_sim=True)
-    R_N_2[index], R_E_2[index] = run_morse_potential(energy, index, num_samples=num_samples, run_sim=True)
+    R_N[index], R_E[index] = run_krc_morse_potential(energy, index, num_samples=num_samples, run_sim=run_sim)
+    R_N_2[index], R_E_2[index] = run_morse_potential(energy, index, num_samples=num_samples, run_sim=run_sim)
 
 plt.semilogx(energies, R_N, label='R_N Morse-Kr-C H-Ni, Es=1.5eV', color='purple')
 plt.semilogx(energies, R_N_2, label='R_N Morse H-Ni, Es=1.5eV', color='green')
diff --git a/src/bca.rs b/src/bca.rs
index 76b80df..61e747c 100644
--- a/src/bca.rs
+++ b/src/bca.rs
@@ -622,7 +622,7 @@ pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parame
     let g = |r: f64| -> f64 {interactions::distance_of_closest_approach_function_singularity_free(r, a, Za, Zb, relative_energy, impact_parameter, interaction_potential)*
         interactions::scaling_function(r, impact_parameter, interaction_potential)};
 
-    let upper_bound = impact_parameter + interactions::crossing_point_doca(interaction_potential);
+    let upper_bound = 10.0*impact_parameter + interactions::crossing_point_doca(interaction_potential);
 
     let roots = match derivative_free {
         true => find_roots_with_secant_polishing(&g, &f, 1e-15, upper_bound,
@@ -636,8 +636,7 @@ pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parame
             truncation_threshold, interval_limit, far_from_zero)
         }
     }.with_context(|| format!("Numerical error: CPR Rootfinder failed to converge when calculating distance of closest approach for Er = {} eV p = {} A using {}.",
-        relative_energy/EV, impact_parameter/ANGSTROM, interaction_potential))
-    .unwrap();
+        relative_energy/EV, impact_parameter/ANGSTROM, interaction_potential))?;
 
     let max_root = roots.iter().cloned().fold(f64::NAN, f64::max)/a;
 
diff --git a/src/consts.rs b/src/consts.rs
index 829fb7a..111ae14 100644
--- a/src/consts.rs
+++ b/src/consts.rs
@@ -44,4 +44,6 @@ pub const TABLE_SIZE: usize = Z_MAX*(Z_MAX + 1)/2;
 /// Gauss-Legendre Points
 pub const GL_X: [f64; 5] = [0./2. + 1./2., -0.538469/2. + 1./2., 0.538469/2. + 1./2., -0.90618/2. + 1./2., 0.90618/2. + 1./2.];
 /// Gauss-Legendre Weights
-pub const GL_W: [f64; 5] = [0.568889/2., 0.478629/2., 0.478629/2., 0.236927/2., 0.236927/2.];
\ No newline at end of file
+pub const GL_W: [f64; 5] = [0.568889/2., 0.478629/2., 0.478629/2., 0.236927/2., 0.236927/2.];
+/// Enforced maximum number density
+pub const MAX_DENSITY: f64 = 1e32;
\ No newline at end of file
diff --git a/src/geometry.rs b/src/geometry.rs
index fb121c6..f7f31f9 100644
--- a/src/geometry.rs
+++ b/src/geometry.rs
@@ -65,8 +65,9 @@ impl Geometry for Mesh0D {
 
         let densities: Vec<f64> = input.densities.iter().map(|element| element/(length_unit).powi(3)).collect();
         assert!(!densities.is_empty(), "Input Error: density list empty.");
-
+        
         let total_density: f64 = densities.iter().sum();
+        assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density);
 
         let energy_barrier_thickness = 1./total_density.cbrt()/SQRTPI*2.;
 
@@ -154,6 +155,7 @@ impl Geometry for Mesh1D {
         let densities: Vec<Vec<f64>> = geometry_input.densities
             .iter()
             .map( |row| row.iter().map(|element| element/(length_unit).powi(3)).collect() ).collect();
+    
 
         //Assert all layer density lists are equal length
         assert!(
@@ -170,6 +172,7 @@ impl Geometry for Mesh1D {
             layer_bottom += layer_thickness*length_unit;
 
             let total_density: f64 = densities.iter().sum();
+            assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density);
             let concentrations: Vec<f64> = densities.iter().map(|&density| density/total_density).collect::<Vec<f64>>();
 
             layers.push(Layer1D::new(layer_top, layer_bottom, densities, concentrations, ck));
@@ -308,6 +311,7 @@ impl Geometry for HomogeneousMesh2D {
         let densities: Vec<f64> = input.densities.iter().map(|element| element/(length_unit).powi(3)).collect();
 
         let total_density: f64 = densities.iter().sum();
+        assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density);
 
         let energy_barrier_thickness = 1./total_density.cbrt()/SQRTPI*2.;
 
@@ -476,6 +480,7 @@ impl Geometry for Mesh2D {
                 y3*length_unit,
             );
             let total_density: f64 = densities.iter().sum();
+            assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density);
             let concentrations: Vec<f64> = densities.iter().map(|&density| density/total_density).collect::<Vec<f64>>();
 
             cells.push(Cell2D::new(coordinate_set_converted, densities, concentrations, ck));
diff --git a/src/input.rs b/src/input.rs
index 0a95d6e..b79b87b 100644
--- a/src/input.rs
+++ b/src/input.rs
@@ -212,6 +212,10 @@ fn default_rootfinder() -> Vec<Vec<Rootfinder>> {
     vec![vec![Rootfinder::DEFAULTNEWTON]]
 }
 
+fn default_output_dir() -> String {
+    "./".to_string()
+}
+
 /// Rustbca's internal representation of the simulation-level options.
 #[cfg(not(feature = "distributions"))]
 #[derive(Deserialize, Clone)]
@@ -250,7 +254,9 @@ pub struct Options {
     #[serde(default = "default_false")]
     pub track_energy_losses: bool,
     #[serde(default = "default_seed")]
-    pub seed: i32
+    pub seed: i32,
+    #[serde(default = "default_output_dir")]
+    pub output_dir: String,
 }
 
 #[cfg(not(feature = "distributions"))]
@@ -275,6 +281,7 @@ impl Options {
             track_displacements: false,
             track_energy_losses: false,
             seed: default_seed(),
+            output_dir: default_output_dir(),
         }
     }
 }
@@ -331,7 +338,9 @@ pub struct Options {
     pub y_num: usize,
     pub z_num: usize,
     #[serde(default = "default_seed")]
-    pub seed: i32
+    pub seed: i32,
+    #[serde(default = "default_output_dir")]
+    pub output_dir: String,
 }
 
 #[cfg(feature = "distributions")]
@@ -370,7 +379,8 @@ impl Options {
             x_num: 0,
             y_num: 0,
             z_num: 0,
-            seed: default_seed()
+            seed: default_seed(),
+            output_dir: default_output_dir(),
         }
     }
 }
@@ -422,8 +432,8 @@ pub fn process_input_file<T: Geometry>(input: <T as Geometry>::InputFileFormat)
     assert!(material.m.len() == material.Eb.len(), "Input error: material input arrays of unequal length.");
     assert!(material.m.len() == material.Es.len(), "Input error: material input arrays of unequal length.");
 
-    if material.interaction_index.len() <= 1 {
-        material.interaction_index = vec![0; material.m.len()];
+    if material.Ed.len() <= 1 {
+        material.interaction_index = vec![material.interaction_index[0]; material.m.len()];
     }
 
     if material.Ed.len() <= 1 {
diff --git a/src/lib.rs b/src/lib.rs
index e300238..4ebb80f 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -1196,19 +1196,19 @@ pub fn compound_bca_list_tracked_py<'py>(energies: Vec<f64>, ux: Vec<f64>, uy: V
 ///Returns:
 ///    vx, vy, vz (float): final x, y, and z velocity in m/s. When ion implants in material, vx, vy, and vz will all be zero.
 #[pyfunction]
-pub fn reflect_single_ion_py<'py>(ion: &Bound<'py, PyDict>, target: &Bound<'py, PyDict>, vx: f64, vy: f64, vz: f64) -> (f64, f64, f64){
+pub fn reflect_single_ion_py<'py>(ion: &Bound<'py, PyDict>, target: &Bound<'py, PyDict>, vx: f64, vy: f64, vz: f64) -> PyResult<(f64, f64, f64)> {
     
-    let Z1: f64 = ion.get_item("Z").unwrap().expect("Error: Cannot get key 'Z' from ion dict.").extract().unwrap();
-    let m1: f64 = ion.get_item("m").unwrap().expect("Error: Cannot get key 'm' from ion dict.").extract().unwrap();
-    let Es1: f64 = ion.get_item("Es").unwrap().expect("Error: Cannot get key 'Es' from ion dict.").extract().unwrap();
-    let Ec1: f64 = ion.get_item("Ec").unwrap().expect("Error: Cannot get key 'Ec' from ion dict.").extract().unwrap();
-
-    let Z2: f64 = target.get_item("Z").unwrap().expect("Error: Cannot get key 'Z' from target dict.").extract().unwrap();
-    let m2: f64 = target.get_item("m").unwrap().expect("Error: Cannot get key 'm' from target dict.").extract().unwrap();
-    let Es2: f64 = target.get_item("Es").unwrap().expect("Error: Cannot get key 'Es' from target dict.").extract().unwrap();
-    let Ec2: f64 = target.get_item("Ec").unwrap().expect("Error: Cannot get key 'Ec' from target dict.").extract().unwrap();
-    let Eb2: f64 = target.get_item("Eb").unwrap().expect("Error: Cannot get key 'Eb' from target dict.").extract().unwrap();
-    let n2: f64 = target.get_item("n").unwrap().expect("Error: Cannot get key 'n' from target dict.").extract().unwrap();
+    let Z1: f64 = get_value_from_dict!(ion, "Z")?;
+    let m1: f64 = get_value_from_dict!(ion, "m")?;
+    let Es1: f64 = get_value_from_dict!(ion, "Es")?;
+    let Ec1: f64 = get_value_from_dict!(ion, "Ec")?;
+
+    let Z2: f64 = get_value_from_dict!(target, "Z")?;
+    let m2: f64 = get_value_from_dict!(target, "m")?;
+    let Es2: f64 = get_value_from_dict!(target, "Es")?;
+    let Ec2: f64 = get_value_from_dict!(target, "Ec")?;
+    let Eb2: f64 = get_value_from_dict!(target, "Eb")?;
+    let n2: f64 = get_value_from_dict!(target, "n")?;
 
     assert!(vx > 0.0, "Input error: vx must be greater than zero for incident particles to hit surface at x=0.");
 
@@ -1271,9 +1271,9 @@ pub fn reflect_single_ion_py<'py>(ion: &Bound<'py, PyDict>, target: &Bound<'py,
     let vz2 = output[0].dir.z*reflected_velocity;
 
     if output[0].E > 0.0 && output[0].dir.x < 0.0 && output[0].left && output[0].incident {
-        (vx2, vy2, vz2)
+        Ok((vx2, vy2, vz2))
     } else {
-        (0.0, 0.0, 0.0)
+        Ok((0.0, 0.0, 0.0))
     }
 }
 
diff --git a/src/material.rs b/src/material.rs
index 6be7b0e..26610ac 100644
--- a/src/material.rs
+++ b/src/material.rs
@@ -229,17 +229,6 @@ impl <T: Geometry> Material<T> {
         }
     }
 
-    ///The minimum cutoff energy of all species that make up the material.
-    pub fn minimum_cutoff_energy(&self) -> f64 {
-        let mut min_Ec = self.Ec.iter().sum::<f64>();
-        for Ec in self.Ec.iter() {
-            if min_Ec > *Ec {
-                min_Ec = *Ec;
-            }
-        }
-        min_Ec
-    }
-
     ///Choose the parameters of a target atom as a concentration-weighted random draw from the species in the triangle that contains or is nearest to (x, y).
     pub fn choose(&self, x: f64, y: f64, z: f64, rng: &mut ChaCha8Rng) -> (usize, f64, f64, f64, f64, f64, usize) {
         let random_number: f64 = rng.random::<f64>();
diff --git a/src/output.rs b/src/output.rs
index 463c0c5..b94d25e 100644
--- a/src/output.rs
+++ b/src/output.rs
@@ -1,5 +1,25 @@
 use super::*;
 use std::fs::File;
+use std::path::Path;
+#[cfg(feature = "distributions")]
+extern crate ndarray;
+#[cfg(feature = "distributions")]
+use ndarray::prelude::*;
+
+macro_rules! open_output_file {
+    ($output_path:expr, $name:expr, $options:expr) => {{
+    //Open output files for streaming output
+    let file_path = $output_path.join(format!("{}{}{}", $options.name, $name, ".output"));
+    let file = OpenOptions::new()
+        .write(true)
+        .create(true)
+        .truncate(true)
+        .open(file_path)
+        .context(format!("Could not open {} output file.", $options.name,))
+        .unwrap();
+    BufWriter::with_capacity($options.write_buffer_size, file)
+    }}
+}
 
 #[derive(Clone, Debug)]
 pub struct OutputUnits {
@@ -21,12 +41,6 @@ pub fn energy_angle_from_particle(particle: &particle::Particle, units: &OutputU
     (energy, angle)
 }
 
-#[cfg(feature = "distributions")]
-extern crate ndarray;
-
-#[cfg(feature = "distributions")]
-use ndarray::prelude::*;
-
 /// Distribution tracker for tracking EADs and implantation distributions
 #[derive(Serialize)]
 #[cfg(feature = "distributions")]
@@ -67,17 +81,22 @@ impl Distributions {
     /// Write distributions to toml
     pub fn print(&self, options: &Options) {
 
+        let output_path = Path::new(&options.output_dir);
+        assert!(output_path.try_exists().unwrap());
+        assert!(output_path.is_dir());
+
+        let file_path = output_path.join(format!("{}{}{}", options.name, "distributions", ".toml"));
         let distribution_output_file = OpenOptions::new()
             .write(true)
             .create(true)
             .truncate(true)
-            .open(format!("{}{}", options.name, "distributions.toml"))
+            .open(file_path)
             .context("Could not open distributions output file.")
             .unwrap();
-        let mut distribution_file_stream = BufWriter::with_capacity(8000, distribution_output_file);
+        let mut distribution_file_stream = BufWriter::with_capacity(options.write_buffer_size, distribution_output_file);
         let toml = toml::to_string(&self).unwrap();
+        
         writeln!(distribution_file_stream, "{}", toml).unwrap();
-
     }
 
     /// Updates distributions with a single particle
@@ -203,14 +222,15 @@ pub struct SummaryPerSpecies {
 impl SummaryPerSpecies {
     pub fn new(options: &Options) -> SummaryPerSpecies {
 
-        let summary_output_file = OpenOptions::new()
-            .write(true)
-            .create(true)
-            .truncate(true)
-            .open(format!("{}{}", options.name, "summary.output"))
-            .context("Could not open output file.")
-            .unwrap();
-        let writer = BufWriter::with_capacity(8000, summary_output_file);
+        let output_path = Path::new(&options.output_dir);
+        assert!(output_path.try_exists().unwrap());
+        assert!(output_path.is_dir());
+
+        let writer = open_output_file!(
+            output_path,
+            "summary",
+            options
+        );
 
         SummaryPerSpecies {
             m: vec![],
@@ -259,69 +279,53 @@ impl SummaryPerSpecies {
 
 /// Open list output files for streaming write
 pub fn open_output_lists(options: &Options) -> OutputListStreams {
-    //Open output files for streaming output
-    let reflected_file = OpenOptions::new()
-        .write(true)
-        .create(true)
-        .truncate(true)
-        .open(format!("{}{}", options.name, "reflected.output"))
-        .context("Could not open output file.")
-        .unwrap();
-    let reflected_file_stream = BufWriter::with_capacity(options.write_buffer_size, reflected_file);
 
-    let sputtered_file = OpenOptions::new()
-        .write(true)
-        .create(true)
-        .truncate(true)
-        .open(format!("{}{}", options.name, "sputtered.output"))
-        .context("Could not open output file.")
-        .unwrap();
-    let sputtered_file_stream = BufWriter::with_capacity(options.write_buffer_size, sputtered_file);
+    let output_path = Path::new(&options.output_dir);
+    assert!(output_path.try_exists().unwrap());
+    assert!(output_path.is_dir());
 
-    let deposited_file = OpenOptions::new()
-        .write(true)
-        .create(true)
-        .truncate(true)
-        .open(format!("{}{}", options.name, "deposited.output"))
-        .context("Could not open output file.")
-        .unwrap();
-    let deposited_file_stream = BufWriter::with_capacity(options.write_buffer_size, deposited_file);
-
-    let trajectory_file = OpenOptions::new()
-        .write(true)
-        .create(true)
-        .truncate(true)
-        .open(format!("{}{}", options.name, "trajectories.output"))
-        .context("Could not open output file.")
-        .unwrap();
-    let trajectory_file_stream = BufWriter::with_capacity(options.write_buffer_size, trajectory_file);
-
-    let trajectory_data = OpenOptions::new()
-        .write(true)
-        .create(true)
-        .truncate(true)
-        .open(format!("{}{}", options.name, "trajectory_data.output"))
-        .context("Could not open output file.")
-        .unwrap();
-    let trajectory_data_stream = BufWriter::with_capacity(options.write_buffer_size, trajectory_data);
-
-    let displacements_file = OpenOptions::new()
-        .write(true)
-        .create(true)
-        .truncate(true)
-        .open(format!("{}{}", options.name, "displacements.output"))
-        .context("Could not open output file.")
-        .unwrap();
-    let displacements_file_stream = BufWriter::with_capacity(options.write_buffer_size, displacements_file);
-
-    let energy_loss_file = OpenOptions::new()
-        .write(true)
-        .create(true)
-        .truncate(true)
-        .open(format!("{}{}", options.name, "energy_loss.output"))
-        .context("Could not open output file.")
-        .unwrap();
-    let energy_loss_file_stream = BufWriter::with_capacity(options.write_buffer_size, energy_loss_file);
+    //Open output files for streaming output
+    let reflected_file_stream = open_output_file!(
+        output_path,
+        "reflected",
+        options
+    );
+
+    let sputtered_file_stream = open_output_file!(
+        output_path,
+        "sputtered",
+        options
+    );
+
+    let deposited_file_stream = open_output_file!(
+        output_path,
+        "deposited",
+        options
+    );
+
+    let trajectory_file_stream = open_output_file!(
+        output_path,
+        "trajectories",
+        options
+    );
+
+    let trajectory_data_stream = open_output_file!(
+        output_path,
+        "trajectory_data",
+        options
+    );
+
+    let displacements_file_stream = open_output_file!(
+        output_path,
+        "displacements",
+        options
+    );
+
+    let energy_loss_file_stream = open_output_file!(
+        output_path,
+        "energy_loss",
+        options
+    );
 
     OutputListStreams {
         reflected_file_stream,
diff --git a/src/parry.rs b/src/parry.rs
index 8686683..ae1e886 100644
--- a/src/parry.rs
+++ b/src/parry.rs
@@ -77,6 +77,7 @@ impl Geometry for ParryBall {
         let electronic_stopping_correction_factor = input.electronic_stopping_correction_factor;
         let densities: Vec<f64> = input.densities.iter().map(|element| element/(length_unit).powi(3)).collect();
         let total_density: f64 = densities.iter().sum();
+        assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density);
         let energy_barrier_thickness = total_density.powf(-1./3.)/SQRTPI*2.;
         let concentrations: Vec<f64> = densities.iter().map(|&density| density/total_density).collect::<Vec<f64>>();
         let radius = input.radius*length_unit;
@@ -201,6 +202,7 @@ impl Geometry for ParryTriMesh {
         let electronic_stopping_correction_factor = input.electronic_stopping_correction_factor;
         let densities: Vec<f64> = input.densities.iter().map(|element| element/(length_unit).powi(3)).collect();
         let total_density: f64 = densities.iter().sum();
+        assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density);
         let energy_barrier_thickness = total_density.powf(-1./3.)/SQRTPI*2.;
         let concentrations: Vec<f64> = densities.iter().map(|&density| density/total_density).collect::<Vec<f64>>();
         let points = input.vertices.iter().map(|p| Point::new(p[0]*length_unit , p[1]*length_unit , p[2]*length_unit)).collect();
diff --git a/src/sphere.rs b/src/sphere.rs
index a9326fc..bb51d05 100644
--- a/src/sphere.rs
+++ b/src/sphere.rs
@@ -63,13 +63,14 @@ impl Geometry for Sphere {
             "NM" => NM,
             "M" => 1.,
             _ => input.length_unit.parse()
-                .unwrap_or_else(|_| panic!("Input errror: could nor parse length unit {}. Use a valid float or one of ANGSTROM, NM, MICRON, CM, MM, M",
+                .unwrap_or_else(|_| panic!("Input errror: could nor parse length unit {}. Use a valid float or one of ANGSTROM, NM, MICRON, CM, MM, M",
                         &input.length_unit.as_str())),
         };
 
         let electronic_stopping_correction_factor = input.electronic_stopping_correction_factor;
         let densities: Vec<f64> = input.densities.iter().map(|element| element/(length_unit).powi(3)).collect();
         let total_density: f64 = densities.iter().sum();
+        assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density);
         let energy_barrier_thickness = total_density.powf(-1./3.)/SQRTPI*2.;
         let concentrations: Vec<f64> = densities.iter().map(|&density| density/total_density).collect::<Vec<f64>>();
         let radius = input.radius*length_unit;
diff --git a/src/tests.rs b/src/tests.rs
index d346226..afba5de 100644
--- a/src/tests.rs
+++ b/src/tests.rs
@@ -230,6 +230,7 @@ fn test_distributions() {
         track_displacements: false,
         track_energy_losses: true,
         seed: 0,
+        output_dir: ".".to_string(),
         energy_min: 0.0,
         energy_max: 10.0,
         energy_num: 11,
@@ -885,6 +886,7 @@ fn test_momentum_conservation() {
                             track_displacements: false,
                             track_energy_losses: false,
                             seed: 0,
+                            output_dir: ".".to_string(),
                         };
 
                         #[cfg(feature = "distributions")]
@@ -907,6 +909,7 @@ fn test_momentum_conservation() {
                             track_displacements: false,
                             track_energy_losses: false,
                             seed: 0,
+                            output_dir: ".".to_string(),
                             energy_min: 0.0,
                             energy_max: 10.0,
                             energy_num: 11,
@@ -1119,6 +1122,7 @@ fn test_quadrature() {
         track_displacements: false,
         track_energy_losses: false,
         seed: 0,
+        output_dir: ".".to_string(),
     };
 
     #[cfg(feature = "distributions")]
@@ -1141,6 +1145,7 @@ fn test_quadrature() {
         track_displacements: false,
         track_energy_losses: false,
         seed: 0,
+        output_dir: ".".to_string(),
         energy_min: 0.0,
         energy_max: 10.0,
         energy_num: 11,