some fixes and updates

Cargo.toml

@@ -23,4 +23,7 @@ path = "src/widom.rs"
 name = "surface_tension"
 path = "src/surface_tension.rs"
 
+[[bin]]
+name = "density_z"
+path = "src/density_z.rs"
 

src/density_z.rs

@@ -0,0 +1,97 @@
+mod trajectory;
+use trajectory::*;
+use std::env;
+
+fn main() {
+    // open file
+    let args: Vec<String> = env::args().collect();
+
+    let mut filename : String = "montecarlo.xyz".to_string();
+    let mut skip_frames : usize = 0;
+    let mut slabs : usize = 256;
+    for i in 0..args.len() {
+        if args[i] == "-f" {
+            filename = args[i+1].clone();
+        } else if args[i] == "-s" {
+            skip_frames = args[i+1].parse::<usize>().unwrap();
+        } else if args[i] == "--slabs" {
+            slabs = args[i+1].parse::<usize>().unwrap();
+        }
+    }
+
+    let mut trj_reader = TrjReader::new(&filename);
+
+    // skip some frames
+    println!("# Skipping {} frames.", skip_frames);
+    trj_reader.skip(skip_frames);
+    let mut frame = trj_reader.next_frame();
+
+    let slab_height = frame.box_z / slabs as f64;
+    let slab_volume = slab_height * frame.box_x * frame.box_y;
+    println!("# Density calculation with {} slabs (height={})", slabs, slab_height);
+
+    let mut slab_particles : Vec<f64> = vec![0.0; slabs];
+    let mut frame_count : usize = 0;
+
+    // loop over frames
+    loop  {
+        frame_count += 1;
+        for i in 0..frame.num_particles {
+            let slab_no : usize = get_slab_number_for_position(frame.rz[i], slab_height) - 1;
+            slab_particles[slab_no] += 1.0;
+        }
+
+        if !trj_reader.update_with_next(&mut frame) {
+            break;
+        }
+    }
+
+    println!("# Averaged over {} frames", frame_count);
+    println!("# Position    Density");
+    for i in 0..slab_particles.len() {
+        let p_max = slab_height * (i as f64 + 1.0);
+        let position =(p_max + p_max-slab_height) / 2.0;
+        let particles = slab_particles[i] / frame_count as f64;
+        let density = particles / slab_volume;
+        println!("{}\t{}\t{}", position, density, particles);
+    }
+
+}
+
+pub fn get_slab_number_for_position(z: f64, slab_height: f64) -> usize {
+    let mut slab = 0;
+    let mut z_counter = 0.0_f64;
+    while z_counter  <= z {
+        z_counter += slab_height;
+        slab += 1;
+    }
+    return slab;
+}
+
+#[test]
+fn test_get_slab_no() {
+    let z = 3.5;
+    let slab_height = 1.0;
+    assert_eq!(4, get_slab_number_for_position(z, slab_height));
+
+    let z = 0.0;
+    let slab_height = 1.0;
+    assert_eq!(1, get_slab_number_for_position(z, slab_height));
+
+    let z = 0.1;
+    let slab_height = 0.2;
+    assert_eq!(1, get_slab_number_for_position(z, slab_height));
+
+    let z = 0.0001;
+    let slab_height = 0.2;
+    assert_eq!(1, get_slab_number_for_position(z, slab_height));
+
+
+    let z = 0.19999999_f64;
+    let slab_height = 0.2;
+    assert_eq!(1, get_slab_number_for_position(z, slab_height));
+
+    let z = 1.0;
+    let slab_height = 1.0;
+    assert_eq!(2, get_slab_number_for_position(z, slab_height));
+}

src/energy.rs

@@ -74,6 +74,11 @@ fn test_get_particle_distance_squared() {
     let a = get_particle_distance_squared(x3,y3,z3,x4,y4,z4, 1.418983411970384,1.418983411970384,2.83796682394076,0.709491705985192,0.709491705985192,1.418983411970384);
     let b = get_particle_distance_squared(x4,y4,z4,x3,y3,z3, 1.418983411970384,1.418983411970384,2.83796682394076,0.709491705985192,0.709491705985192,1.418983411970384);
     assert!( (a-b).abs() < 0.0000000001);
+
+    let (x1, y1, z1) = (1.0, 1.0, 1.0);
+    let (x2, y2, z2) = (99.0, 99.0, 99.0);
+    let dist = get_particle_distance_squared(x1,y1,z1,x2,y2,z2, 100.0, 100.0, 100.0, 50.0, 50.0, 50.0);
+    assert!(dist - 12.0 < 0.00001);
 }
 
 pub fn eval_pair_energy(dist_squared: f64, e_shift: f64) -> (f64, f64) {

src/main.rs

@@ -136,12 +136,16 @@ fn main() {
         if rx.len() == num_particles { break; }
     }
 
-    // scale box in z for vacuum space above
+    // 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 };
@@ -196,11 +200,11 @@ fn main() {
         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 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 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 }
+        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);

src/trajectory.rs

@@ -178,4 +178,30 @@ impl TrjReader {
         return true;
     }
 
+    // skip x frames
+    pub fn skip(&mut self, skip: usize) {
+        if skip < 1 { return };
+
+        // find number of particles
+        let buffer_string = &mut String::new();
+        match self.reader.read_line(buffer_string) {
+            Ok(size) => {
+                if size == 0 {
+                    panic!("no new frame!")
+                }
+            },
+            Err(e) => panic!("no new frame!")
+        }
+        let first_line_vec : Vec<&str> = buffer_string.split_whitespace().collect();
+        let num_particles = first_line_vec[0].parse::<usize>().unwrap();
+
+        let lines_to_skip = (num_particles + 1) * skip - 1;
+        let mut skipped = 0;
+        loop {
+            self.reader.read_line(&mut String::new());
+            skipped += 1;
+            if skipped == lines_to_skip { break; }
+        }
+    }
+
 }

src/widom.rs

@@ -13,33 +13,59 @@ static MKSA_PLANCKS_CONSTANT_H : f64 = 1.0;
 static MASS : f64 = 1.0;
 
 const INSERTIONS_PER_STEP : usize = 1000;
-const RUN_AVG_SIZE : usize = 10;
+const RUN_AVG_SIZE : usize = 100;
 
 const SHIFT : bool = true;
 
+
+
 fn main() {
-    // open file
+    // parse args
     let args: Vec<String> = env::args().collect();
-    let mut trj_reader = TrjReader::new(&args[1]);
+    let mut filename = "montecarlo.xyz".to_string();
+    let mut skip: usize = 0;
+
+    // liquid phase boundaries
+    let mut liquid_start = 0.0;
+    let mut liquid_end = 0.0;
+    for i in 0..args.len() {
+        if args[i] == "-f" {
+            filename = args[i + 1].clone();
+        } else if args[i] == "-s" {
+            skip = args[i + 1].parse::<usize>().unwrap();
+        } else if args[i] == "-ls" {
+            liquid_start = args[i + 1].parse::<f64>().unwrap();
+        } else if args[i] == "-le" {
+            liquid_end = args[i + 1].parse::<f64>().unwrap();
+        }
+    }
+
+    // open file and skip to requiested position
+    let mut trj_reader = TrjReader::new(&filename);
+    if skip > 0 { trj_reader.skip(skip) };
 
     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;
 
-    // account for gas phase slab
+    // rnd number generator for particle insertion
-    let gas_slab = args[2].parse::<f64>().unwrap();
+    let mut rng = rand::thread_rng();
-    let liquid_height = frame.box_z / (gas_slab + 1.0);
+
-    let liquid_volume = volume / (gas_slab + 1.0);
+
+    // calculate liquid and gas volume
+    let liquid_height : f64 = liquid_end - liquid_start;
+    let liquid_volume = liquid_height * frame.box_x * frame.box_y;
     let gas_volume = volume - liquid_volume;
 
-    let wave = MKSA_PLANCKS_CONSTANT_H / (2.0 * std::f64::consts::PI * MASS * frame.temperature / frame.lj_eps).sqrt();
+    // lambda = h/sqrt(2*pi*m*kb*T) = (2*pi*m*kb*T/h^2)^3/2
+    // this is lambda^3
+    let wave3 =  ((2.0 * std::f64::consts::PI * MASS * frame.temperature / frame.lj_eps)/MKSA_PLANCKS_CONSTANT_H.powi(2)).powf(3.0/2.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 };
 
-    let mut rng = rand::thread_rng();
-
 
     // average counters
     let mut frame_count = 0;
@@ -47,8 +73,6 @@ fn main() {
     let mut ideal_sum_gas = 0.0;
     let mut widom_sum_liquid = 0.0;
     let mut ideal_sum_liquid = 0.0;
-
-    // loop over frames
     let mut avg_count = 0;
     loop  {
         frame_count += 1;
@@ -57,32 +81,40 @@ fn main() {
         let mut liquid_count = 0.0;
         let mut gas_count = 0.0;
         for i in 0..frame.num_particles {
-            if frame.rz[i] < liquid_height { liquid_count += 1.0; }
+            if frame.rz[i] > liquid_start && frame.rz[i] < liquid_end { liquid_count += 1.0; }
             else { gas_count += 1.0; }
         }
 
+
         // test particle insertion multiple times
         for i in 0..INSERTIONS_PER_STEP {
             avg_count += 1;
 
-            // liquid test partciles
+            // liquid test partcile
             let lx = frame.box_x * rng.gen::<f64>();
             let ly = frame.box_y * rng.gen::<f64>();
-            let lz = liquid_height * rng.gen::<f64>();
+            let lz = liquid_start + (liquid_height * rng.gen::<f64>());
-
             let widom_e_liquid = get_particle_insertion_energy(&frame.rx, &frame.ry, &frame.rz, frame.num_particles, lx, ly, lz, frame.box_x, frame.box_y, frame.box_z, cutoff_sqr, e_shift);
             widom_sum_liquid += (-beta*widom_e_liquid).exp();
+//            widom_sum_liquid += widom_e_liquid;
 
-            // test particle gas energy
+            // gas test particle
+            let upper_or_lower = rng.gen::<bool>();
             let gx = frame.box_x * rng.gen::<f64>();
             let gy = frame.box_y * rng.gen::<f64>();
-            let gz = ((frame.box_z - liquid_height) * rng.gen::<f64>()) + liquid_height;
+            let gz;
+            if upper_or_lower {
+                gz = rng.gen::<f64>() * liquid_start;
+            } else {
+                gz = rng.gen::<f64>() * (frame.box_z - liquid_end) + liquid_end;
+            }
             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();
+//            widom_sum_gas += widom_e_gas;
 
-            // calculate ideal gas potential for both phases
+            // calculate ideal gas potential
-            ideal_sum_gas += frame.temperature / frame.lj_eps * (gas_volume/(wave.powi(3)* gas_count)).ln();
+            ideal_sum_gas += -frame.temperature / frame.lj_eps * ((gas_volume/gas_count) * wave3).ln();
-            ideal_sum_liquid += frame.temperature / frame.lj_eps * (liquid_volume/(wave.powi(3)* liquid_count)).ln();
+            ideal_sum_liquid += frame.temperature / frame.lj_eps * ((liquid_volume/liquid_count) * wave3).ln();
 
         }
 
@@ -90,9 +122,14 @@ fn main() {
         if avg_count / INSERTIONS_PER_STEP > RUN_AVG_SIZE {
             let ideal_gas = ideal_sum_gas / avg_count as f64;
             let ideal_liquid = ideal_sum_liquid / 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 excess_gas = -(widom_sum_gas/avg_count as f64).ln()/beta;
             let excess_liquid = -(widom_sum_liquid/avg_count as f64).ln()/beta;
-            println!("Frame {}\tgas {}\tliquid {}\t particles gas/liquid:{}/{}", frame_count, ideal_gas + excess_gas, ideal_liquid + excess_liquid, gas_count, liquid_count);
+            let gas_total = ideal_gas + excess_gas;
+            let liquid_total = ideal_liquid + excess_liquid;
+            println!("Frame {}\tg_ex: {:5}\tl_ex: {:5}\tg_tot: {:5}\tl_tot: {:5}\t\tnparticles: {}/{}", frame_count, excess_gas, excess_liquid, 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);
 
             // reset averages for next round
             avg_count = 0;