comments

src/energy.rs

@@ -1,4 +1,7 @@
+ #![allow(dead_code)]
 
+/// Calculates the total energy and virial of a system containing num_particles with coords rx,ry,rz
+/// of size l_x, l_y, l_z and given cutoff + corrections
 pub fn get_total_energy(rx: &[f64], ry: &[f64], rz: &[f64], num_particles: usize, l_x: f64, l_y: f64, l_z: f64, cutoff_squared: f64, e_corr: f64, e_shift: f64) -> (f64, f64) {
     let mut energy = 0.0;
     let mut virial = 0.0;
@@ -19,6 +22,8 @@ pub fn get_total_energy(rx: &[f64], ry: &[f64], rz: &[f64], num_particles: usize
     return (energy, virial);
 }
 
+/// Calculates the particle energy and virial for particle at p_index in system containing num_particles with coords rx,ry,rz
+/// of size l_x, l_y, l_z and given cutoff + corrections
 pub fn get_particle_energy(rx: &[f64], ry: &[f64], rz: &[f64], p_index: usize, num_particles: usize, l_x: f64, l_y: f64, l_z: f64, cutoff_squared: f64, e_shift: f64) -> (f64, f64) {
     let mut energy = 0.0;
     let mut virial = 0.0;
@@ -38,6 +43,7 @@ pub fn get_particle_energy(rx: &[f64], ry: &[f64], rz: &[f64], p_index: usize, n
     return (energy, virial);
 }
 
+// squared distance between 2 particles regarding the minimum image convention
 pub fn get_particle_distance_squared(x1: f64,y1: f64,z1: f64,x2: f64,y2: f64,z2: f64, l_x: f64, l_y: f64, l_z: f64, hl_x: f64, hl_y: f64, hl_z: f64) -> f64 {
     let mut dx = (x1 - x2).abs();
     let mut dy = (y1 - y2).abs();
@@ -50,13 +56,17 @@ pub fn get_particle_distance_squared(x1: f64,y1: f64,z1: f64,x2: f64,y2: f64,z2:
     if dz > hl_z { dz -= l_y }
         else if dz < -hl_z { dz += l_z}
     return dx*dx + dy*dy + dz*dz;
-
 }
 
-pub fn get_particle_distance(x1: f64,y1: f64,z1: f64,x2: f64,y2: f64,z2: f64, l_x: f64, l_y: f64, l_z: f64, hl_x: f64, hl_y: f64, hl_z: f64) -> f64{
-    return get_particle_distance_squared(x1,y1,z1,x2,y2,z2, l_x, l_y, l_z, hl_x, hl_y, hl_z).sqrt();
+// one dimensional distance with applied minimum image convention
+pub fn get_distance_with_pbc(x1: f64, x2: f64, length: f64, half_length: f64) -> f64 {
+    let mut d = (x1-x2).abs();
+    if d > half_length { d -= length }
+    else if d < -half_length { d += length }
+    return d;
 }
 
+
 #[test]
 fn test_get_particle_distance_squared() {
     let (x1, y1, z1) = (0.0, 0.0, 0.0);
@@ -81,6 +91,7 @@ fn test_get_particle_distance_squared() {
     assert!(dist - 12.0 < 0.00001);
 }
 
+/// calculate the lj energy and virial between two particles from given square distance
 pub fn eval_pair_energy(dist_squared: f64, e_shift: f64) -> (f64, f64) {
     let r6 = ::LJ_SIG/(dist_squared * dist_squared * dist_squared);
     let r62 = r6*r6;
@@ -101,11 +112,12 @@ fn test_eval_pair_energy() {
     assert!( (e - 224.0).abs() < 0.00001, "{}",  e);
 }
 
-
-pub fn eval_virial(distance: f64, LJ_EPS: f64, LJ_SIG: f64) -> f64 {
-    let r7 = (LJ_SIG/distance).powi(7);
-    let r13 = (LJ_SIG/distance).powi(13);
-    return 24.0 * LJ_EPS / LJ_SIG * ( r7-2.0*r13 );
+/// calculate the virial between two particles from given square distance. If energy is required too,
+/// see eval_pair_energy which does energy and virial
+pub fn eval_virial(distance: f64, lj_eps: f64, lj_sig: f64) -> f64 {
+    let r7 = (lj_sig/distance).powi(7);
+    let r13 = (lj_sig/distance).powi(13);
+    return 24.0 * lj_eps / lj_sig * ( r7-2.0*r13 );
 }
 
 #[test]

src/main.rs

@@ -11,13 +11,19 @@ use argparse::{ArgumentParser, Store, StoreFalse, StoreTrue};
 mod trajectory;
 use trajectory::*;
 
-
+// LJ params
 const LJ_EPS : f64 = 1.0;
 const LJ_SIG : f64 = 1.0;
 
-
+ // intended acceptance rate = 33%
 const TRIES_INTENDED : f64 = 3.0;
+
+// factor for auto displacement scaling
 const DISP_SCALE_FACTOR : f64 = 0.1;
+const SCALE_INTERVAL : usize = 5000;
+
+const EQUILIBRATION_OUTPUT_INTERVAL : usize = 5000;
+const SAMPLING_OUTPUT_INTERVAL : usize = 5000;
 
 // easy printing to stderr
 macro_rules! println_stderr(
@@ -27,7 +33,278 @@ macro_rules! println_stderr(
     } }
 );
 
-fn parse_cmd_args(NUM_STEPS: &mut usize, NUM_MINIM_STEPS: &mut usize,
+fn main() {
+    println_stderr!("");
+    println_stderr!("################################################################");
+    println_stderr!("##################  LJ Monte Carlo Simulation  #################");
+    println_stderr!("################################################################");
+    println_stderr!("");
+
+    /** Definition of default run parameters **/
+    let mut eq_steps  = 1000000;
+    let mut sample_steps = 100000;
+
+    let mut num_particles: usize = 512;
+    let mut density = 0.7;
+    let mut temperature = 0.9;
+
+    let mut cutoff = 3.0;
+
+    let mut TAILCORR : bool = true;
+    let mut SHIFT: bool = true;
+
+    let mut displacement = 0.1; // max particle displacement in one dimension
+    let mut SCALE: bool = true; // switch for displacement scaling
+
+    // scale factor in z for vaccuum space
+    let mut vacuum_slab = 0.0;
+
+    // output config
+    let mut output_prefix = "montecarlo".to_string(); // .xyz will be append
+    let mut output_interval : i64 = 100;
+    let mut output_minim : bool = false;
+
+    // parse cmd line arguments and override defaults
+    parse_cmd_args(&mut sample_steps, &mut eq_steps, &mut num_particles,
+                   &mut density, &mut temperature,
+                   &mut cutoff, &mut displacement, &mut SCALE, &mut TAILCORR, &mut SHIFT,
+                   &mut output_prefix, &mut output_interval, &mut output_minim,
+                   &mut vacuum_slab);
+
+
+    /** Initialize the system **/
+    let beta = 1.0/temperature;
+
+    let mut volume = (num_particles as f64)/ density;
+    let length  = volume.cbrt();
+    let (l_x, l_y, mut l_z) = (length, length, length);
+
+    let cutoff_squared = cutoff * cutoff;
+    let max_displacement = length / 2.0; // displacement wont be scaled over that
+
+    // initialize randomness - TODO seed?
+    let mut rng = rand::thread_rng();
+    let particle_range = Range::new(0, num_particles);
+
+    // randomly place particles in the box
+    let mut rx : Vec<f64> = vec![];
+    let mut ry : Vec<f64> = vec![];
+    let mut rz : Vec<f64> = vec![];
+    loop {
+        rx.push(l_x * rng.gen::<f64>());
+        ry.push(l_y * rng.gen::<f64>());
+        rz.push(l_z * rng.gen::<f64>());
+        if rx.len() == num_particles { break; }
+    }
+
+    // scale box in z for vacuum space and move particles in the middle of the box
+    if vacuum_slab > 0.0 {
+        let scale = vacuum_slab + 1.0;
+        l_z *= scale;
+        volume *= scale;
+        density /= scale;
+        let move_z = l_z/scale*vacuum_slab/2.0;
+        for i in 0..num_particles {
+            rz[i] += move_z;
+        }
+    }
+
+    // calculation of shift and tailcorrections
+    let e_shift = if SHIFT { 4.0 * LJ_EPS * ( (LJ_SIG/cutoff).powi(12) - (LJ_SIG/cutoff).powi(6) ) } else { 0.0 };
+    let e_corr = if TAILCORR { 8.0/3.0*std::f64::consts::PI*density*LJ_EPS*LJ_SIG.powi(3)*((1.0/3.0*(LJ_SIG/cutoff).powi(9)) - (LJ_SIG/cutoff).powi(3)) } else { 0.0 };
+    let p_corr = if TAILCORR { 16.0/3.0*std::f64::consts::PI*density.powi(2)*LJ_EPS*LJ_SIG.powi(3)*((2.0/3.0*(LJ_SIG/cutoff).powi(9)) - (LJ_SIG/cutoff).powi(3)) } else { 0.0 };
+
+    println_stderr!("Particles: {}, Density: {}, Temperature: {}", num_particles, density, temperature);
+    println_stderr!("System volume: {:8.3}, Dimensions {:.3}/{:.3}/{:.3}", volume, l_x, l_y, l_z);
+    println_stderr!("Minimization steps: {}, Sampling steps: {}", eq_steps, sample_steps);
+    println_stderr!("LJ params eps: {}, sigma: {}, cutoff: {}", LJ_EPS, LJ_SIG, cutoff);
+    println_stderr!("Tailcorr: {:8.3}, Shift: {:8.3}, Pressurecorr: {:8.3}", e_corr, e_shift, p_corr);
+
+    // energy and average sums
+    let (mut energy, mut virial) = get_total_energy(&rx, &ry, &rz, num_particles, l_x, l_y, l_z, cutoff_squared, e_corr, e_shift);
+    let mut energy_sum = 0.0;
+    let mut virial_sum = 0.0;
+    let mut step_counter = 0;
+    let mut accept_counter = 0;
+
+
+    // prepare and write first trajectory frame
+    let mut trajectory : XYZTrajectory = XYZTrajectory::new(&format!("{}.xyz", output_prefix));
+    if output_minim { trajectory.write(&rx, &ry, &rz, num_particles, l_x, l_y, l_z, temperature, LJ_EPS, LJ_SIG, cutoff, true); }
+
+
+    println_stderr!("");
+    println_stderr!("################################################################");
+    println_stderr!("########################  Equilibration  #######################");
+    println_stderr!("################################################################");
+    println_stderr!("");
+
+
+    // START OF METROPOLIS
+    /*****************************************************************************************/
+
+    for step in 0..eq_steps+sample_steps {
+
+        // select rnd particle
+        let rnd_index = particle_range.ind_sample(&mut rng);
+
+        // store old position
+        let oldX = rx[rnd_index];
+        let oldY = ry[rnd_index];
+        let oldZ = rz[rnd_index];
+
+        // old particle energy
+        let (old_particle_energy, old_particle_virial) = get_particle_energy(&rx, &ry, &rz, rnd_index, num_particles, l_x, l_y, l_z, cutoff_squared, e_shift);
+
+        // rnd displacement and PBC
+        rx[rnd_index] += ( rng.gen::<f64>() - 0.5 ) * displacement;
+        ry[rnd_index] += ( rng.gen::<f64>() - 0.5 ) * displacement;
+        rz[rnd_index] += ( rng.gen::<f64>() - 0.5 ) * displacement;
+        if rx[rnd_index] < 0.0 { rx[rnd_index] += l_x }
+        if rx[rnd_index] >= l_x { rx[rnd_index] -= l_x }
+        if ry[rnd_index] < 0.0 { ry[rnd_index] += l_y }
+        if ry[rnd_index] >= l_y { ry[rnd_index] -= l_y }
+        if rz[rnd_index] < 0.0 { rz[rnd_index] += l_z }
+        if rz[rnd_index] >= l_z { rz[rnd_index] -= l_z }
+
+        // calculate energy difference
+        let (new_particle_energy, new_particle_virial) = get_particle_energy(&rx, &ry, &rz, rnd_index, num_particles, l_x, l_y, l_z, cutoff_squared, e_shift);
+
+        let dE = new_particle_energy - old_particle_energy;
+
+        // acceptance rule
+        if dE < 0.0 || rng.gen::<f64>() < (-beta * dE).exp() {
+            accept_counter += 1;
+            energy += dE;
+            virial += new_particle_virial - old_particle_virial;
+
+            // recalculate total energy every 1000 steps to account for rounding errors in particle energy function
+            if step % 10000 == 0 {
+                let (e, v) = get_total_energy(&rx, &ry, &rz, num_particles, l_x, l_y, l_z, cutoff_squared, e_corr, e_shift);
+                energy = e;
+                virial = v;
+            }
+        } else {
+            // restore old positions if move is rejected
+            rx[rnd_index] = oldX;
+            ry[rnd_index] = oldY;
+            rz[rnd_index] = oldZ;
+        }
+
+        // update average sums
+        step_counter += 1;
+        energy_sum += energy;
+        virial_sum += virial;
+
+        // reset average sums for sampling
+        if step == eq_steps-1 {
+            println_stderr!("");
+            println_stderr!("################################################################");
+            println_stderr!("##########################  Sampling  ##########################");
+            println_stderr!("################################################################");
+            println_stderr!("");
+            step_counter = 0;
+            accept_counter = 0;
+            energy_sum = 0.0;
+            virial_sum = 0.0;
+        }
+
+        // Everything below here is not part of the metropolis sampling (extras)
+
+        // displacement scaling during equilibration for good acceptance ratios
+        if SCALE && step < eq_steps && step % SCALE_INTERVAL == 0 {
+            let tries_per_step : f64 = step_counter as f64 /accept_counter as f64;
+
+            // will increase the max displacement if the acceptance rate is too high and vice versa
+            let scale_factor = (TRIES_INTENDED/tries_per_step * DISP_SCALE_FACTOR).abs();
+            if tries_per_step < TRIES_INTENDED - 0.2 && displacement < max_displacement {
+                displacement += displacement * scale_factor;
+            } else if tries_per_step > TRIES_INTENDED + 0.2 && displacement > 0.0 {
+                displacement -= displacement * scale_factor;
+            }
+            step_counter = 0;
+            accept_counter = 0;
+            energy_sum = 0.0;
+        }
+
+        // print some output during equilibration
+        if step < eq_steps && step_counter % EQUILIBRATION_OUTPUT_INTERVAL == 0 && step != 0 {
+            let tries_per_step : f64 = step_counter as f64 /accept_counter as f64;
+            let acceptance_rate = 1.0/tries_per_step * 100.0;
+            let avg_energy = energy_sum / step_counter as f64;
+            let avg_virial = virial_sum / step_counter as f64;
+            println_stderr!("Eq {:<10} Energy: {:<30.3} Virial: {:<30.3} Accept.: {:<4.1}%   dr: {:.3}", step+1, avg_energy, avg_virial, acceptance_rate, displacement);
+        }
+
+        // print some output during sampling
+        if step > eq_steps && step_counter % SAMPLING_OUTPUT_INTERVAL == 0 {
+            println_stderr!("Step  {:<10} Energy: {:<30.3}", step_counter, energy);
+        }
+
+        // write trajectory
+        if step as i64 % output_interval == 0 {
+            if step > eq_steps || output_minim {
+                trajectory.write(&rx, &ry, &rz, num_particles, l_x, l_y, l_z, temperature, LJ_EPS, LJ_SIG, cutoff, true);
+            }
+        }
+    }
+
+    // END OF METROPOLIS
+    /*****************************************************************************************/
+    println_stderr!("Done sampling!");
+
+    let final_energy = energy_sum/step_counter as f64;
+    let particle_energy = final_energy / num_particles as f64;
+    let final_virial = virial_sum / 3.0 / step_counter as f64 / volume;
+    let pressure = virial_sum / 3.0 / step_counter as f64 / volume + density * temperature + p_corr;
+    let final_acceptance_rate = 1.0/((accept_counter as f64)/(step_counter as f64)) * 100.0;
+
+    println_stderr!("");
+    println_stderr!("################################################################");
+    println_stderr!("##########################  Results  ###########################");
+    println_stderr!("################################################################");
+    println_stderr!("");
+    println!(
+"Minimization: {}
+Steps: {}
+
+# Lennard Jones Params
+epsilon: {}
+sigma: {}
+cutoff: {}
+
+# System
+Particles: {}
+Density: {}
+Temperature: {}
+Volume: {}
+Box dimension: {:.3}/{:.3}/{:.3}
+Max Displacement: {}
+
+# Correction
+Energy correction: {}
+Shift: {}
+P-Correction: {}
+
+# Averages
+Tries: {}
+Accepted: {}
+Acceptance: {:.2}%
+Energy: {}
+Energy per particle: {}
+Virial: {}
+Pressure: {}",
+         eq_steps, sample_steps,
+        LJ_EPS, LJ_SIG, cutoff,
+        num_particles, density, temperature, volume, l_x, l_y, l_z, displacement,
+        e_corr, e_shift, p_corr,
+        step_counter, accept_counter, final_acceptance_rate, final_energy, particle_energy, final_virial, pressure);
+
+    trajectory.write(&rx, &ry, &rz, num_particles, l_x, l_y, l_z, temperature, LJ_EPS, LJ_SIG, cutoff, true);
+}
+
+// Parse command line arguments
+fn parse_cmd_args(NUM_STEPS: &mut usize, NUM_eq_steps: &mut usize,
                   NUM_PARTICLES: &mut usize, DENSITY: &mut f64, TEMPERATURE: &mut f64,
                   CUTOFF: &mut f64, MAX_DISP_START: &mut f64, SCALE: &mut bool, TAILCORR: &mut bool, SHIFT: &mut bool,
                   OUTPUT_PREFIX: &mut String, OUTPUT_INTERVAL: &mut i64, OUTPUT_MINIM: &mut bool,
@@ -37,7 +314,7 @@ fn parse_cmd_args(NUM_STEPS: &mut usize, NUM_MINIM_STEPS: &mut usize,
     ap.refer(NUM_STEPS)
         .add_option(&["-n", "--nsteps"], Store,
                     "Simulation steps: Number of steps for averaging" );
-    ap.refer(NUM_MINIM_STEPS)
+    ap.refer(NUM_eq_steps)
         .add_option(&["-m", "--nminimsteps"], Store,
                     "Minimization steps: Number of steps before averaging starts");
     ap.refer(NUM_PARTICLES)
@@ -78,255 +355,3 @@ fn parse_cmd_args(NUM_STEPS: &mut usize, NUM_MINIM_STEPS: &mut usize,
                     "Disable lj shifting");
     ap.parse_args_or_exit();
 }
-
-fn main() {
-
-    // define all the stuff
-    let mut minim_steps  = 1000000;
-    let mut sample_steps = 100000;
-
-    let mut num_particles: usize = 512;
-    let mut density = 0.7;
-    let mut temperature = 0.9;
-
-    let mut cutoff = 3.0;
-    let mut displacement = 0.1;
-
-    let mut  TAILCORR : bool = true;
-    let mut SHIFT: bool = true;
-    let mut SCALE: bool = true;
-
-    let mut vacuum_slab = 0.0;
-
-    let mut output_prefix = "montecarlo".to_string();
-    let mut output_interval : i64 = 100;
-    let mut output_minim : bool = false;
-    parse_cmd_args(&mut sample_steps, &mut minim_steps, &mut num_particles,
-                   &mut density, &mut temperature,
-                   &mut cutoff, &mut displacement, &mut SCALE, &mut TAILCORR, &mut SHIFT,
-                   &mut output_prefix, &mut output_interval, &mut output_minim,
-                   &mut vacuum_slab);
-
-    println_stderr!("");
-    println_stderr!("################################################################");
-    println_stderr!("##################  LJ Monte Carlo Simulation  #################");
-    println_stderr!("################################################################");
-    println_stderr!("");
-
-
-    // initialize stuff
-    let beta = 1.0/temperature;
-    let mut volume = (num_particles as f64)/ density;
-    let length  = volume.cbrt();
-    let (l_x, l_y, mut l_z) = (length, length, length);
-
-    let cutoff_squared = cutoff * cutoff;
-    let max_displacement = length / 2.0;
-
-    let mut rng = rand::thread_rng();
-    let particle_range = Range::new(0, num_particles);
-
-    let mut rx : Vec<f64> = vec![];
-    let mut ry : Vec<f64> = vec![];
-    let mut rz : Vec<f64> = vec![];
-    loop {
-        rx.push(l_x * rng.gen::<f64>());
-        ry.push(l_y * rng.gen::<f64>());
-        rz.push(l_z * rng.gen::<f64>());
-        if rx.len() == num_particles { break; }
-    }
-
-    // scale box in z for vacuum space and move particles in the middle of the box
-    if vacuum_slab > 0.0 {
-        let scale = vacuum_slab + 1.0;
-        l_z *= scale;
-        volume *= scale;
-        density /= scale;
-        let move_z = l_z/scale*vacuum_slab/2.0;
-        for i in 0..num_particles {
-            rz[i] += move_z;
-        }
-    }
-
-    let e_shift = if SHIFT { 4.0 * LJ_EPS * ( (LJ_SIG/cutoff).powi(12) - (LJ_SIG/cutoff).powi(6) ) } else { 0.0 };
-    let e_corr = if TAILCORR { 8.0/3.0*std::f64::consts::PI*density*LJ_EPS*LJ_SIG.powi(3)*((1.0/3.0*(LJ_SIG/cutoff).powi(9)) - (LJ_SIG/cutoff).powi(3)) } else { 0.0 };
-    let p_corr = if TAILCORR { 16.0/3.0*std::f64::consts::PI*density.powi(2)*LJ_EPS*LJ_SIG.powi(3)*((2.0/3.0*(LJ_SIG/cutoff).powi(9)) - (LJ_SIG/cutoff).powi(3)) } else { 0.0 };
-
-    println_stderr!("Particles: {}, Density: {}, Temperature: {}", num_particles, density, temperature);
-    println_stderr!("System volume: {:8.3}, Dimensions {:.3}/{:.3}/{:.3}", volume, l_x, l_y, l_z);
-    println_stderr!("Minimization steps: {}, Sampling steps: {}", minim_steps, sample_steps);
-    println_stderr!("LJ params eps: {}, sigma: {}, cutoff: {}", LJ_EPS, LJ_SIG, cutoff);
-    println_stderr!("Tailcorr: {:8.3}, Shift: {:8.3}, Pressurecorr: {:8.3}", e_corr, e_shift, p_corr);
-
-    let (mut energy, mut virial) = get_total_energy(&rx, &ry, &rz, num_particles, l_x, l_y, l_z, cutoff_squared, e_corr, e_shift);
-    let mut energy_sum = 0.0;
-    let mut virial_sum = 0.0;
-    let mut step_counter = 0;
-    let mut accept_counter = 0;
-
-    println_stderr!("");
-    println_stderr!("################################################################");
-    println_stderr!("#####################  Energy Minimization  ####################");
-    println_stderr!("################################################################");
-    println_stderr!("");
-
-    // prepare and write first trajectory frame
-    let mut trajectory : XYZTrajectory = XYZTrajectory::new(&format!("{}.xyz", output_prefix));
-    if output_minim { trajectory.write(&rx, &ry, &rz, num_particles, l_x, l_y, l_z, temperature, LJ_EPS, LJ_SIG, cutoff, true); }
-
-//    let vacuum_scale_step = (minim_steps as f64 * 0.1) as usize;
-    for step in 0..minim_steps+sample_steps {
-
-        // first minimizate solvent phase, then add vacuum slab
-//        if vacuum_slab > 0.0 && step == vacuum_scale_step { // increase space in z
-//            let scale = vacuum_slab + 1.0;
-//            l_z *= scale;
-//            volume *= scale;
-//            density /= scale;
-//        }
-        // select rnd particle
-        let rnd_index = particle_range.ind_sample(&mut rng);
-
-        // store old position
-        let oldX = rx[rnd_index];
-        let oldY = ry[rnd_index];
-        let oldZ = rz[rnd_index];
-
-        // old particle energy
-        let (old_particle_energy, old_particle_virial) = get_particle_energy(&rx, &ry, &rz, rnd_index, num_particles, l_x, l_y, l_z, cutoff_squared, e_shift);
-
-        // rnd displacement and PBC
-        rx[rnd_index] += ( rng.gen::<f64>() - 0.5 ) * displacement;
-        ry[rnd_index] += ( rng.gen::<f64>() - 0.5 ) * displacement;
-        rz[rnd_index] += ( rng.gen::<f64>() - 0.5 ) * displacement;
-        if rx[rnd_index] < 0.0 { rx[rnd_index] += l_x }
-        if rx[rnd_index] >= l_x { rx[rnd_index] -= l_x }
-        if ry[rnd_index] < 0.0 { ry[rnd_index] += l_y }
-        if ry[rnd_index] >= l_y { ry[rnd_index] -= l_y }
-        if rz[rnd_index] < 0.0 { rz[rnd_index] += l_z }
-        if rz[rnd_index] >= l_z { rz[rnd_index] -= l_z }
-
-        // calculate energy difference
-        let (new_particle_energy, new_particle_virial) = get_particle_energy(&rx, &ry, &rz, rnd_index, num_particles, l_x, l_y, l_z, cutoff_squared, e_shift);
-        let dE = new_particle_energy - old_particle_energy;
-
-        //accept move
-        if rng.gen::<f64>() < (-beta * dE).exp() {
-            accept_counter += 1;
-            if step % 1000 == 0 { // calculate total energy every 1000 steps to account for rounding errors
-                let (e, v) = get_total_energy(&rx, &ry, &rz, num_particles, l_x, l_y, l_z, cutoff_squared, e_corr, e_shift);
-                energy = e;
-                virial = v;
-            } else {
-                energy += dE;
-                virial += new_particle_virial - old_particle_virial;
-            }
-        } else { // or restore old position
-            rx[rnd_index] = oldX;
-            ry[rnd_index] = oldY;
-            rz[rnd_index] = oldZ;
-        }
-
-        // update sums for averaging
-        step_counter += 1;
-        energy_sum += energy;
-        virial_sum += virial;
-
-        // print some output during minimization
-        if step < minim_steps && step_counter % 5000 == 0 && step != 0 {
-            let tries_per_step : f64 = step_counter as f64 /accept_counter as f64;
-            let acceptance_rate = 1.0/tries_per_step * 100.0;
-            let avg_energy = energy_sum / step_counter as f64;
-            let avg_virial = virial_sum / step_counter as f64;
-            println_stderr!("Minim {:<10} Energy: {:<30.3} Virial: {:<30.3} Accept.: {:<4.1}%   dr: {:.3}", step+1, avg_energy, avg_virial, acceptance_rate, displacement);
-
-            if SCALE {
-                let scale_factor = (TRIES_INTENDED/tries_per_step * DISP_SCALE_FACTOR).abs();
-                if tries_per_step < TRIES_INTENDED - 0.2 && displacement < max_displacement {
-                    displacement += displacement * scale_factor;
-                } else if tries_per_step > TRIES_INTENDED + 0.2 && displacement > 0.0 {
-                    displacement -= displacement * scale_factor;
-                }
-                step_counter = 0;
-                accept_counter = 0;
-                energy_sum = 0.0;
-            }
-        }
-
-        // reset sums for sampling
-        if step == minim_steps-1 {
-            println_stderr!("");
-            println_stderr!("################################################################");
-            println_stderr!("##########################  Sampling  ##########################");
-            println_stderr!("################################################################");
-            println_stderr!("");
-            step_counter = 0;
-            accept_counter = 0;
-            energy_sum = 0.0;
-            virial_sum = 0.0;
-        }
-
-        if step > minim_steps && step_counter % 5000 == 0 {
-            println_stderr!("Step  {:<10} Energy: {:<30.3}", step_counter, energy);
-        }
-
-        // write trajectory maybe
-        if step as i64 % output_interval == 0 {
-            if step > minim_steps || output_minim {
-                trajectory.write(&rx, &ry, &rz, num_particles, l_x, l_y, l_z, temperature, LJ_EPS, LJ_SIG, cutoff, true);
-            }
-        }
-
-    }
-
-
-    let final_energy = energy_sum/step_counter as f64;
-    let particle_energy = final_energy / num_particles as f64;
-    let final_virial = virial_sum / 3.0 / step_counter as f64 / volume;
-    let pressure = virial_sum / 3.0 / step_counter as f64 / volume + density * temperature + p_corr;
-    let final_acceptance_rate = 1.0/((accept_counter as f64)/(step_counter as f64)) * 100.0;
-
-    println_stderr!("Done sampling!");
-    println_stderr!("");
-    println_stderr!("################################################################");
-    println_stderr!("##########################  Results  ###########################");
-    println_stderr!("################################################################");
-    println_stderr!("");
-    println!(
-"Minimization: {}
-Steps: {}
-
-# Lennard Jones Params
-epsilon: {}
-sigma: {}
-cutoff: {}
-
-# System
-Particles: {}
-Density: {}
-Temperature: {}
-Volume: {}
-Box dimension: {:.3}/{:.3}/{:.3}
-Max Displacement: {}
-
-# Correction
-Energy correction: {}
-Shift: {}
-P-Correction: {}
-
-# Averages
-Tries: {}
-Accepted: {}
-Acceptance: {:.2}%
-Energy: {}
-Energy per particle: {}
-Virial: {}
-Pressure: {}",
-         minim_steps, sample_steps,
-        LJ_EPS, LJ_SIG, cutoff,
-        num_particles, density, temperature, volume, l_x, l_y, l_z, displacement,
-        e_corr, e_shift, p_corr,
-        step_counter, accept_counter, final_acceptance_rate, final_energy, particle_energy, final_virial, pressure);
-
-    trajectory.write(&rx, &ry, &rz, num_particles, l_x, l_y, l_z, temperature, LJ_EPS, LJ_SIG, cutoff, true);
-}

src/surface_tension.rs

@@ -7,32 +7,14 @@ use std::env;
 const LJ_EPS : f64 = 1.0;
 const LJ_SIG : f64 = 1.0;
 
-
-fn get_distance_with_pbc(x1: f64, x2: f64, length: f64, half_length: f64) -> f64 {
-    let mut d = (x1-x2).abs();
-    if d > half_length { d -= length }
-    else if d < -half_length { d += length }
-    return d;
-}
-
-fn eval_surface_tension(box_z: f64, p_zz: f64, p_xy: f64) -> f64 {
-    return box_z / 2.0 * (p_zz - p_xy);
-}
-
-#[test]
-fn test_eval_surface_tension() {
-    let expected = 2.0;
-    let result = eval_surface_tension(2.0,5.0,3.0);
-    assert!( (result-expected).abs() < 0.0001, "{}", result );
-}
-
+const AVG_OUTPUT_INTERVAL : usize = 10;
 
 fn main() {
-
-    // parse args
-    let args: Vec<String> = env::args().collect();
     let mut filename = "montecarlo.xyz".to_string();
     let mut skip: usize = 0;
+
+    // parse cmd line args
+    let args: Vec<String> = env::args().collect();
     for i in 0..args.len() {
         if args[i] == "-f" {
             filename = args[i + 1].clone();
@@ -41,17 +23,16 @@ fn main() {
         }
     }
 
-    // open file and skip to requiested position
+    // open file and skip to requested position
     let mut trj_reader = TrjReader::new(&filename);
     if skip > 0 { trj_reader.skip(skip) };
 
-    // trajectory information
+    // read first trajectory and system params
     let mut frame = trj_reader.next_frame();
-    println!("{:?}", frame);
-
     let volume = frame.box_x * frame.box_y * frame.box_z;
     let density = frame.num_particles as f64 / volume;
     let num_particles = frame.num_particles;
+    println!("{:?}", frame);
 
     let box_half_x = frame.box_x / 2.0;
     let box_half_y = frame.box_y / 2.0;
@@ -63,6 +44,7 @@ fn main() {
 
     let variable_without_name = frame.temperature/LJ_EPS * density;
 
+    println!("Calculating surface tension");
     println!("~~~ THIS IS A RUNNING AVERAGE! ~~~");
     loop {
         frame_count += 1;
@@ -71,6 +53,7 @@ fn main() {
         let mut trace_z = 0.0;
         for i in 0..num_particles {
             for j in i+1..num_particles {
+                // this needs some optimization for speed
                 let dist_sqrt = get_particle_distance_squared(frame.rx[i], frame.ry[i], frame.rz[i], frame.rx[j], frame.ry[j], frame.rz[j], frame.box_x, frame.box_y, frame.box_z, box_half_x, box_half_y, box_half_z);
                 let dist = dist_sqrt.sqrt();
                 let dx = get_distance_with_pbc(frame.rx[i], frame.rx[j], frame.box_x, box_half_x);
@@ -88,7 +71,7 @@ fn main() {
         p_z_sum += p_zz;
 
         ///////////////////////////////////
-        if frame_count % 10 == 0 {
+        if frame_count % AVG_OUTPUT_INTERVAL == 0 {
             let p_z_avg = p_z_sum / frame_count as f64;
             let p_xy_avg = p_xy_sum / frame_count as f64;
             let p_diff = p_z_avg - p_xy_avg;
@@ -100,7 +83,20 @@ fn main() {
             // p_xy_sum = 0.0;
         }
 
+        // read next frame
         if !trj_reader.update_with_next(&mut frame) { break }
     }
 
 }
+
+/// calc surface tension from box z size and pressure tensor
+fn eval_surface_tension(box_z: f64, p_zz: f64, p_xy: f64) -> f64 {
+    return box_z / 2.0 * (p_zz - p_xy);
+}
+
+#[test]
+fn test_eval_surface_tension() {
+    let expected = 2.0;
+    let result = eval_surface_tension(2.0,5.0,3.0);
+    assert!( (result-expected).abs() < 0.0001, "{}", result );
+}

src/trajectory.rs

@@ -1,5 +1,5 @@
 #![allow(dead_code)]
-#![allow(unused_must_use)] // hate.
+#![allow(unused_must_use)]
 #![allow(unused_variables)]
 
 use std::error::Error;
@@ -9,7 +9,6 @@ use std::path::Path;
 use std::fmt;
 use std::io::BufReader;
 
-
 pub struct XYZTrajectory {
     file: File,
 
@@ -31,7 +30,6 @@ impl XYZTrajectory {
 
     }
 
-
     pub fn write(&mut self, rx: &[f64], ry: &[f64], rz: &[f64], num_particles: usize, box_x : f64, box_y : f64, box_z : f64, temp: f64 ,lj_eps : f64, lj_sig : f64, lj_cutoff : f64, flush: bool) {
         self.file.write(format!("{} ## Box: {} {} {} Temp: {} LJ: {}/{}/{}\n", num_particles, box_x,box_y,box_z,temp, lj_eps, lj_sig, lj_cutoff).as_bytes());
         for i  in 0..num_particles {

src/widom.rs

@@ -1,3 +1,5 @@
+#![allow(unused_variables)]
+
 mod trajectory;
 use trajectory::*;
 mod energy;
@@ -14,30 +16,19 @@ static MASS : f64 = 1.0;
 
 const RUN_AVG_SIZE : usize = 1;
 
-const SHIFT : bool = true;
-
-fn eval_ideal_potential(temperature: f64, lj_eps: f64, volume: f64, particles: f64, thermal_wavelength3: f64) -> f64 {
-    let density = volume/particles;
-    return -temperature/lj_eps * ( density*thermal_wavelength3 ).ln();
-}
-
-#[test]
-fn test_eval_ideal_potential() {
-    let expected = 12.476649;
-    let result = eval_ideal_potential(2.0,1.0,10.0,512.0,0.1);
-    assert!((expected-result).abs() < 0.00001, "{}", result);
-}
-
 fn main() {
-    // parse args
-    let args: Vec<String> = env::args().collect();
+    // default values
     let mut filename = "montecarlo.xyz".to_string();
-    let mut skip: usize = 0;
-    let mut insertions: usize = 100;
+    let mut skip: usize = 0; // skip some frames
+    let mut insertions: usize = 100; // number of particle insertions per step per phase
 
     // liquid phase boundaries
     let mut liquid_start = 0.0;
     let mut liquid_end = 0.0;
+    let mut shift = true;
+
+    // parse command line arguments
+    let args: Vec<String> = env::args().collect();
     for i in 0..args.len() {
         if args[i] == "-f" {
             filename = args[i + 1].clone();
@@ -49,18 +40,19 @@ fn main() {
             liquid_end = args[i + 1].parse::<f64>().unwrap();
         } else if args[i] == "-n" {
             insertions = args[i+1].parse::<usize>().unwrap();
+        } else if args[i] == "--noshift" {
+            shift = false
         }
     }
 
-    // open file and skip to requiested position
+    // open file and skip to requested position
     let mut trj_reader = TrjReader::new(&filename);
     if skip > 0 { trj_reader.skip(skip) };
 
-    // trajectory information
+    // Get first frame and read system configuration
     let mut frame = trj_reader.next_frame();
     println!("{:?}", frame);
 
-    // get some non changing values
     let volume = frame.box_x * frame.box_y * frame.box_z;
     let beta = 1.0/frame.temperature;
     let cutoff_sqr = frame.lj_cutoff * frame.lj_cutoff;
@@ -77,7 +69,7 @@ fn main() {
     let tw3 =  ((2.0 * std::f64::consts::PI * MASS * frame.temperature / frame.lj_eps)/MKSA_PLANCKS_CONSTANT_H.powi(2)).powf(3.0/2.0);
 
     // LJ shift
-    let e_shift = if SHIFT { 4.0 * LJ_EPS * ( (LJ_SIG/frame.lj_cutoff).powi(12) - (LJ_SIG/frame.lj_cutoff).powi(6) ) } else { 0.0 };
+    let e_shift = if shift { 4.0 * LJ_EPS * ( (LJ_SIG/frame.lj_cutoff).powi(12) - (LJ_SIG/frame.lj_cutoff).powi(6) ) } else { 0.0 };
 
 
     // average counters
@@ -88,8 +80,6 @@ fn main() {
     let mut ideal_pot_liquid_sum = 0.0;
     let mut avg_count = 0;
 
-    // loop over all frames
-    let mut counter = 0;
     loop  {
         frame_count += 1;
 
@@ -116,14 +106,13 @@ fn main() {
             let gx = frame.box_x * rng.gen::<f64>();
             let gy = frame.box_y * rng.gen::<f64>();
             let mut gz = frame.box_z * rng.gen::<f64>();
-            while gz > liquid_start && gz < liquid_end { // retry until we have a particle in gas
-                 gz= frame.box_z * rng.gen::<f64>();
+            while gz > liquid_start && gz < liquid_end { // retry until we hit the gas phase
+                gz= frame.box_z * rng.gen::<f64>();
             }
-
             let widom_e_gas = get_particle_insertion_energy(&frame.rx, &frame.ry, &frame.rz, frame.num_particles, gx, gy, gz, frame.box_x, frame.box_y, frame.box_z, cutoff_sqr, e_shift);
             widom_sum_gas += (-beta*widom_e_gas).exp();
 
-            // calculate ideal potentials
+            // ideal gas potentials
             ideal_pot_gas_sum += eval_ideal_potential(frame.temperature, frame.lj_eps, gas_volume, gas_count, tw3);
             ideal_pot_liquid_sum += eval_ideal_potential(frame.temperature, frame.lj_eps, liquid_volume, liquid_count, tw3);
 
@@ -133,12 +122,14 @@ fn main() {
         if avg_count / insertions > RUN_AVG_SIZE {
             let ideal_gas_potential = ideal_pot_gas_sum / avg_count as f64;
             let ideal_liquid_potential = ideal_pot_liquid_sum / avg_count as f64;
-            let excess_gas = -(widom_sum_gas/avg_count as f64).ln()/beta;
-            let excess_liquid = -(widom_sum_liquid/avg_count as f64).ln()/beta;
-            let mut gas_total = ideal_gas_potential + excess_gas;
-            let liquid_total = ideal_liquid_potential + excess_liquid;
-            println!("Frame {}\tg_ex: {:5}\tl_ex: {:5}\tg_tot: {:5}\tl_tot: {:5}\t\tnparticles: {}/{}", frame_count, excess_gas, excess_liquid, if gas_total.is_infinite() { 0.0 } else { gas_total } , liquid_total, gas_count, liquid_count);
-//            println!("Frame {}\tgas {}\tliquid {}\t particles gas/liquid:{}/{}", frame_count, ideal_gas + excess_gas, ideal_liquid + excess_liquid, gas_count, liquid_count);
+            let excess_gas_potential = -(widom_sum_gas/avg_count as f64).ln()/beta;
+            let excess_liquid_potential = -(widom_sum_liquid/avg_count as f64).ln()/beta;
+            let gas_total = ideal_gas_potential + excess_gas_potential;
+            let liquid_total = ideal_liquid_potential + excess_liquid_potential;
+            println!("Frame {}\tg_ex: {:5}\tl_ex: {:5}\tg_tot: {:5}\tl_tot: {:5}\t\tnparticles: {}/{}",
+                frame_count, excess_gas_potential, excess_liquid_potential,
+                if gas_total.is_infinite() { 0.0 } else { gas_total },
+                liquid_total, gas_count, liquid_count);
 
             // reset averages for next round
             avg_count = 0;
@@ -148,6 +139,7 @@ fn main() {
             widom_sum_liquid = 0.0;
         }
 
+        // jump to next frame
         if !trj_reader.update_with_next(&mut frame) {
             break;
         }
@@ -155,7 +147,7 @@ fn main() {
 
 }
 
-
+/// calculates the energy for a hypothetic particle inserted at x,y,z
 fn get_particle_insertion_energy(rx: &[f64], ry: &[f64], rz: &[f64], num_particles: usize, x: f64, y: f64, z: f64, l_x: f64,l_y: f64, l_z: f64, cutoff_sqr: f64, e_shift: f64) -> f64 {
     let mut energy = 0.0;
     let half_l_x = l_x/2.0;
@@ -170,3 +162,16 @@ fn get_particle_insertion_energy(rx: &[f64], ry: &[f64], rz: &[f64], num_particl
     }
     return energy;
 }
+
+/// calculates the ideal gas chemical potential
+fn eval_ideal_potential(temperature: f64, lj_eps: f64, volume: f64, particles: f64, thermal_wavelength3: f64) -> f64 {
+    let density = volume/particles;
+    return -temperature/lj_eps * ( density*thermal_wavelength3 ).ln();
+}
+
+#[test]
+fn test_eval_ideal_potential() {
+    let expected = 12.476649;
+    let result = eval_ideal_potential(2.0,1.0,10.0,512.0,0.1);
+    assert!((expected-result).abs() < 0.00001, "{}", result);
+}