From e6eaa62a623c5cc9e955060bd595e54819167e59 Mon Sep 17 00:00:00 2001 From: danijoo Date: Thu, 2 Feb 2017 18:56:00 +0100 Subject: [PATCH] surface tension reimplemented as running average (and maybe correct values now) --- src/surface_tension.rs | 105 +++++++++++++++++++++-------------------- 1 file changed, 53 insertions(+), 52 deletions(-) diff --git a/src/surface_tension.rs b/src/surface_tension.rs index e355b31..8d80a37 100644 --- a/src/surface_tension.rs +++ b/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 = 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::().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); -// } \ No newline at end of file