surface tension reimplemented as running average (and maybe correct values now)

src/surface_tension.rs

@@ -2,28 +2,65 @@ mod trajectory;
 use trajectory::*;
 mod energy;
 use energy::*;
+use std::env;
 
 const LJ_EPS : f64 = 1.0;
 const LJ_SIG : f64 = 1.0;
 
 
-fn get_virial(distance_sqr: f64) -> f64 {
-    let r2 = LJ_SIG.powi(2)/distance_sqr;
-    let r6 = r2 * r2 * r2;
-    return 48.0 * LJ_EPS / LJ_SIG * ( r6*r6 - 0.5*r6 )
+fn get_virial(distance: 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]
+fn test_get_viral() {
+    let result = get_virial(2.5);
+    let expected = 0.038999477;
+    assert!( (result - expected) < 0.0001, "{}", result );
 }
 
 fn get_distance_with_pbc(x1: f64, x2: f64, length: f64, half_length: f64) -> f64 {
-    let mut d = x1-x2;
+    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 );
+}
+
+
 fn main() {
 
-    let mut trj_reader = TrjReader::new(&"2phases/1kk_100step.xyz".to_string());
+    // parse args
+    let args: Vec<String> = env::args().collect();
+    let mut filename = "montecarlo.xyz".to_string();
+    let mut skip: usize = 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();
+        }
+    }
+
+    // open file and skip to requiested position
+    let mut trj_reader = TrjReader::new(&filename);
+    if skip > 0 { trj_reader.skip(skip) };
+
+    // trajectory information
     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;
@@ -35,16 +72,15 @@ fn main() {
 
     let mut frame_count = 0;
     let mut p_xy_sum = 0.0;
-//    let mut p_y_sum = 0.0;
     let mut p_z_sum = 0.0;
 
     let variable_without_name = frame.temperature/LJ_EPS * density;
 
+    println!("~~~ THIS IS A RUNNING AVERAGE! ~~~");
     loop {
         frame_count += 1;
 
         let mut trace_xy = 0.0;
-//        let mut trace_y = 0.0;
         let mut trace_z = 0.0;
         for i in 0..num_particles {
             for j in i+1..num_particles {
@@ -53,66 +89,31 @@ fn main() {
                 let dx = get_distance_with_pbc(frame.rx[i], frame.rx[j], frame.box_x, box_half_x);
                 let dy = get_distance_with_pbc(frame.ry[i], frame.ry[j], frame.box_y, box_half_y);
                 let dz = get_distance_with_pbc(frame.rz[i], frame.rz[j], frame.box_z, box_half_z);
-                let virial = get_virial(dist_sqrt);
+                let virial = get_virial(dist);
                 trace_xy += (dx * dx + dy * dy) / dist * virial;
-//                trace_y += (dy * dy) / dist * virial;
                 trace_z += (dz * dz) / dist * virial;
             }
         }
-        let p_xy = variable_without_name - 1.0/(2.0*volume)*(trace_xy/num_particles as f64);
-//        let p_yy = variable_without_name - 1.0/volume*(trace_y/num_particles);
-        let p_zz = variable_without_name - 1.0/volume*(trace_z/num_particles as f64);
+        let p_xy = variable_without_name - 1.0/(2.0*volume)*(trace_xy);
+        let p_zz = variable_without_name - 1.0/volume*(trace_z);
 
         p_xy_sum += p_xy;
-//        p_y_sum += p_yy;
         p_z_sum += p_zz;
 
         ///////////////////////////////////
-
-        if frame_count % 100 == 0 {
-            print!(".");
-        }
-
-        if frame_count % 100 == 0 {
+        if frame_count % 10 == 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;
-            let surface_tension = (frame.box_z/2.0)*p_diff;
-            println!("{}  zz: {}   xy: {}    diff: {}    tension:   {}", frame_count, p_z_avg, p_xy_avg, p_diff, surface_tension);
-
-            frame_count = 0;
-            p_z_sum = 0.0;
-            p_xy_sum = 0.0;
-//            trace_xy_sum = 0.0;
-//            trace_z_sum = 0.0;
-//            p_tangial_sum = 0.0;
-//            p_normal_sum = 0.0
+            let surface_tension = eval_surface_tension(frame.box_z, p_z_avg, p_xy_avg);
+            println!("Frame {}\t\tzz: {:.5}\txy: {:.5}\tdifference: {:.5}\t\ttension: {:.5}", frame_count, p_z_avg, p_xy_avg, p_diff, surface_tension);
 
+            // frame_count = 0;
+            // p_z_sum = 0.0;
+            // p_xy_sum = 0.0;
         }
 
         if !trj_reader.update_with_next(&mut frame) { break }
     }
 
 }
-
-//        for slab in 0..NUM_SLABS {
-//
-//            // calculate slab density
-//            let slab_end_z = slab_height * slab as f64;
-//            let slab_start_z = slab_end_z - slab_height;
-//            let mut slab_particle_count = 0;
-//            for i in 0..frame.num_particles {
-//                if slab_start_z > frame.rz[i] && frame.rz[i] < slab_end_z {
-//                    slab_particle_count += 1;
-//                }
-//            }
-//            let slab_density = slab_particle_count as f64 / slab_volume;
-//
-//            for i in 0..frame.num_particles {
-//                for j in i+1..frame.num_particles {
-//
-//                }
-//            }
-//
-//            println!("Slab {}\tDensity: {}", slab, slab_density);
-//        }