[WIP] Additional input validation checks/warnings

examples/make_input_file_and_run.py

@@ -144,6 +144,7 @@
             }
         ]
     ],
+    'seed': 0 # if <0, will generate a seed from thread-local PRNG; if >0, will be used as seed to PRNG
 }
 
 # material parameters are per-species

src/input.rs

@@ -426,6 +426,14 @@ where <T as Geometry>::InputFileFormat: Deserialize<'static> + 'static {
         material.Ed = vec![material.Ed[0]; material.m.len()];
     }
 
+    for Ed_ in &material.Ed {
+        assert!(*Ed_ >= 0.0, "Input Error: Ed = {}; Ed cannot be less than zero.", Ed_);
+    }
+
+    for Eb_ in &material.Eb {
+        assert!(*Eb_ >= 0.0, "Input Error: Eb = {}; Eb cannot be less than zero.", Eb_);
+    }
+
     //Check that incompatible options are not on simultaneously
     if options.high_energy_free_flight_paths {
         assert!(options.electronic_stopping_mode == ElectronicStoppingMode::INTERPOLATED,
@@ -657,27 +665,27 @@ where <T as Geometry>::InputFileFormat: Deserialize<'static> + 'static {
                             y: match y {
                                 Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*length_unit},
                                 Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap();  uniform.sample(&mut rng)*length_unit},
-                                Distributions::POINT(x) => x*length_unit,
+                                Distributions::POINT(y) => y*length_unit,
                             },
                             z: match z {
                                 Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*length_unit},
                                 Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap();  uniform.sample(&mut rng)*length_unit},
-                                Distributions::POINT(x) => x*length_unit,
+                                Distributions::POINT(z) => z*length_unit,
                             },
                             ux: match cosx {
                                 Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*length_unit},
                                 Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap();  uniform.sample(&mut rng)*length_unit},
-                                Distributions::POINT(x) => x*length_unit
+                                Distributions::POINT(ux) => ux
                             },
                             uy: match cosy {
                                 Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*length_unit},
                                 Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap();  uniform.sample(&mut rng)*length_unit},
-                                Distributions::POINT(x) => x*length_unit,
+                                Distributions::POINT(uy) => uy,
                             },
                             uz: match cosz {
                                 Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*length_unit},
                                 Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap();  uniform.sample(&mut rng)*length_unit},
-                                Distributions::POINT(x) => x*length_unit,
+                                Distributions::POINT(uz) => uz,
                             },
                             interaction_index,
                             tag: 0,

src/material.rs

@@ -298,7 +298,7 @@ impl <T: Geometry> Material<T> {
 
             let stopping_power = match electronic_stopping_mode {
                 //Biersack-Varelas Interpolation
-                ElectronicStoppingMode::INTERPOLATED => 1./(1./S_high + 1./S_low)*ck,
+                ElectronicStoppingMode::INTERPOLATED => 1./(1./S_high + 1./(S_low*ck)),
                 //Oen-Robinson
                 ElectronicStoppingMode::LOW_ENERGY_LOCAL => S_low*ck,
                 //Lindhard-Scharff

src/tests.rs

@@ -975,15 +975,15 @@ fn test_momentum_conservation() {
                         println!();
 
                         //These values are in  [angstrom amu / second], so very large.
-                        assert!(approx_eq!(f64, initial_momentum.x/ANGSTROM/AMU, final_momentum.x/ANGSTROM/AMU, epsilon = 1000.));
-                        assert!(approx_eq!(f64, initial_momentum.y/ANGSTROM/AMU, final_momentum.y/ANGSTROM/AMU, epsilon = 1000.));
-                        assert!(approx_eq!(f64, initial_momentum.z/ANGSTROM/AMU, final_momentum.z/ANGSTROM/AMU, epsilon = 1000.));
+                        assert!(approx_eq!(f64, initial_momentum.x/ANGSTROM/AMU, final_momentum.x/ANGSTROM/AMU, epsilon = 10.));
+                        assert!(approx_eq!(f64, initial_momentum.y/ANGSTROM/AMU, final_momentum.y/ANGSTROM/AMU, epsilon = 10.));
+                        assert!(approx_eq!(f64, initial_momentum.z/ANGSTROM/AMU, final_momentum.z/ANGSTROM/AMU, epsilon = 10.));
 
                         assert!(!particle_1.E.is_nan());
                         assert!(!particle_2.E.is_nan());
                         assert!(!initial_momentum.x.is_nan());
-                        assert!(!initial_momentum.x.is_nan());
-                        assert!(!initial_momentum.x.is_nan());
+                        assert!(!initial_momentum.y.is_nan());
+                        assert!(!initial_momentum.z.is_nan());
                     }
                 }
             }