From cd94c50fb208f356dd8ccf7636d309a1713d188f Mon Sep 17 00:00:00 2001 From: danijoo Date: Mon, 23 Jan 2017 08:59:29 +0100 Subject: [PATCH] surface tension and widom particle insertion (WIP) --- src/energy.rs | 14 +++-- src/main.rs | 2 +- src/surface_tension.rs | 115 ++++++++++++++++++++++++++++++++++ src/trajectory.rs | 9 ++- src/widom.rs | 138 +++++++++++++++++++++++++++++++++++++++++ 5 files changed, 267 insertions(+), 11 deletions(-) create mode 100644 src/surface_tension.rs create mode 100644 src/widom.rs diff --git a/src/energy.rs b/src/energy.rs index dfdf152..d9891c1 100644 --- a/src/energy.rs +++ b/src/energy.rs @@ -5,8 +5,8 @@ pub fn get_total_energy(rx: &[f64], ry: &[f64], rz: &[f64], num_particles: usize let hl_x = l_x / 2.0; let hl_y = l_y / 2.0; let hl_z = l_z / 2.0; - for i in 0..num_particles-1 { - for j in i+1..num_particles-1 { + for i in 0..num_particles { + for j in i+1..num_particles { let dist_squared = get_particle_distance_squared(rx[i], ry[i],rz[i],rx[j],ry[j],rz[j], l_x, l_y, l_z, hl_x, hl_y, hl_z); if dist_squared < cutoff_squared { let (e,v) = eval_pair_energy(dist_squared, e_shift); @@ -25,7 +25,7 @@ pub fn get_particle_energy(rx: &[f64], ry: &[f64], rz: &[f64], p_index: usize, n let hl_x = l_x / 2.0; let hl_y = l_y / 2.0; let hl_z = l_z / 2.0; - for i in 0..num_particles-1 { + for i in 0..num_particles { if i == p_index { continue; } let dist_squared = get_particle_distance_squared(rx[i], ry[i], rz[i], rx[p_index], ry[p_index], rz[p_index], l_x, l_y, l_z, hl_x, hl_y, hl_z); @@ -38,7 +38,7 @@ pub fn get_particle_energy(rx: &[f64], ry: &[f64], rz: &[f64], p_index: usize, n return (energy, virial); } -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 { +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; let mut dy = y1 - y2; let mut dz = z1 - z2; @@ -53,6 +53,10 @@ fn get_particle_distance_squared(x1: f64,y1: f64,z1: f64,x2: f64,y2: f64,z2: f64 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(); +} + #[test] fn test_get_particle_distance_squared() { let (x1, y1, z1) = (0.0, 0.0, 0.0); @@ -65,7 +69,7 @@ fn test_get_particle_distance_squared() { assert!( (get_particle_distance_squared(x1,y1,z1,x2,y2,z2, 9.0,9.0,9.0, 4.5,4.5,4.5) - 16.0) < 0.00001, "{}", get_particle_distance_squared(x1,y1,z1,x2,y2,z2, 9.0,9.0,9.0, 4.5,4.5,4.5)); } -fn eval_pair_energy(dist_squared: f64, e_shift: f64) -> (f64, f64) { +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; let energy = 4.0 * ::LJ_EPS * (r62 - r6) - e_shift; diff --git a/src/main.rs b/src/main.rs index 56abaa9..a363358 100644 --- a/src/main.rs +++ b/src/main.rs @@ -124,7 +124,7 @@ fn main() { let max_displacement = length / 2.0; let mut rng = rand::thread_rng(); - let particle_range = Range::new(0, num_particles-1); + let particle_range = Range::new(0, num_particles); let mut rx : Vec = vec![]; let mut ry : Vec = vec![]; diff --git a/src/surface_tension.rs b/src/surface_tension.rs new file mode 100644 index 0000000..a648bea --- /dev/null +++ b/src/surface_tension.rs @@ -0,0 +1,115 @@ +mod trajectory; +use trajectory::*; +mod energy; +use energy::*; + +const LJ_EPS : f64 = 1.0; +const LJ_SIG : f64 = 1.0; + + +fn get_derived_pair_potential(distance: f64) -> f64 { + let r = LJ_SIG/distance; + let r3 = r*r*r; + let r4 = r*r*r*r; + return 24.0 * LJ_EPS / LJ_SIG * ( r4*r3 - 2.0*(r4*r4*r4*r) ) +} + +fn get_distance_with_pbc(x1: f64, x2: f64, length: f64, half_length: f64) -> f64 { + let mut d = x1-x2; + if d > half_length { d -= length } + else if d < -half_length { d += length } + return d; +} + +fn main() { + + let mut trj_reader = TrjReader::new(&"2phases/10k_1step.xyz".to_string()); + let mut frame = trj_reader.next_frame(); + + let volume = frame.box_x * frame.box_y * frame.box_z; + let density = frame.num_particles as f64 / volume; + + let box_half_x = frame.box_x / 2.0; + let box_half_y = frame.box_y / 2.0; + let box_half_z = frame.box_z / 2.0; + + let mut frame_count = 0; + let mut trace_xy_sum = 0.0; + let mut trace_z_sum = 0.0; + let mut p_normal_sum = 0.0; + let mut p_tangial_sum = 0.0; + let mut p_diff_sum = 0.0; + + let variable_without_name = frame.temperature/LJ_EPS * density; + loop { + frame_count += 1; + + let mut trace_xy = 0.0; + let mut trace_z = 0.0; + for i in 0..frame.num_particles { + for j in i+1..frame.num_particles { + let dist = get_particle_distance(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 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); + trace_xy += (dx * dx + dy * dy) / dist * get_derived_pair_potential(dist); + trace_z += (dz * dz) / dist * get_derived_pair_potential(dist); + } + } + let p_tangial = variable_without_name - 1.0/(2.0*volume)*trace_xy; + let p_normal = variable_without_name - 1.0/volume*trace_z; + let p_diff = p_normal - p_tangial_sum; + p_diff_sum += p_diff; + p_tangial_sum += p_tangial; + p_normal_sum += p_normal; + trace_xy_sum += trace_xy; + trace_z_sum += trace_z; + + /////////////////////////////////// + + if frame_count % 100 == 0 { + print!("."); + } + + if frame_count % 100 == 0 { +// let trace_xy_avg = trace_xy_sum / frame_count as f64; +// let trace_z_avg = trace_z_sum / frame_count as f64; + let p_tangial = p_tangial_sum / frame_count as f64; + let p_normal = p_normal_sum / frame_count as f64; + let surface_tension = frame.box_z/2.0*(p_diff_sum/frame_count as f64); + println!("{} tangial: {} normal: {} diff: {} tension: {}", frame_count, p_tangial, p_normal, p_diff_sum/frame_count as f64, surface_tension); + + frame_count = 0; + trace_xy_sum = 0.0; + trace_z_sum = 0.0; + p_tangial_sum = 0.0; + p_normal_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 diff --git a/src/trajectory.rs b/src/trajectory.rs index 09fe3ef..cb2a4ea 100644 --- a/src/trajectory.rs +++ b/src/trajectory.rs @@ -64,8 +64,7 @@ pub struct TrjReader { pub reader: BufReader, } impl TrjReader { - pub fn new() -> TrjReader { - let filename = "trajectory.xyz"; + pub fn new(filename: &String) -> TrjReader { let file = File::open(filename).expect("Failed to open file."); let reader : BufReader = BufReader::new(file); return TrjReader { reader:reader }; @@ -95,7 +94,7 @@ impl TrjReader { let mut rx = Vec::new(); let mut ry = Vec::new(); let mut rz = Vec::new(); - for i in 0..num_particles-1{ + for i in 0..num_particles{ let atom_line = &mut String::new(); match self.reader.read_line(atom_line) { Ok(size) => { @@ -132,7 +131,7 @@ impl TrjReader { return false; } }, - Err(e) => panic!("no new frame!") + Err(e) => { return false; } } let first_line_vec : Vec<&str> = buffer_string.split_whitespace().collect(); frame.num_particles = first_line_vec[0].parse::().unwrap(); @@ -144,7 +143,7 @@ impl TrjReader { frame.lj_eps = lj[0].parse::().unwrap(); frame.lj_sig = lj[1].parse::().unwrap(); frame.lj_cutoff = lj[2].parse::().unwrap(); - for i in 0..frame.num_particles-1 { + for i in 0..frame.num_particles { let atom_line = &mut String::new(); match self.reader.read_line(atom_line) { Ok(size) => { diff --git a/src/widom.rs b/src/widom.rs new file mode 100644 index 0000000..22e91c2 --- /dev/null +++ b/src/widom.rs @@ -0,0 +1,138 @@ +mod trajectory; +use trajectory::*; +mod energy; +use energy::*; +extern crate rand; +use rand::Rng; + +const LJ_EPS : f64 = 1.0; +const LJ_SIG : f64 = 1.0; + +static MKSA_PLANCKS_CONSTANT_H : f64 = 1.0; +static MASS : f64 = 1.0; + +const INSERTIONS_PER_STEP : usize = 10; +const RUN_AVG_SIZE : usize = 100; + +const SHIFT : bool = true; + +// returns number of particles in liquid and gas phase as tuble +fn count_particles(rz: &[f64], num_particles: usize, liquid_height: f64) -> (usize, usize) { + let mut liquid_count = 0; + let mut gas_count = 0; + for i in 0..num_particles { + if rz[i] < liquid_height { liquid_count += 1; } + else { gas_count += 1; } + } + return (liquid_count, gas_count); +} + +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; + let half_l_y = l_y/2.0; + let half_l_z = l_z/2.0; + for i in 0..num_particles { + + let dist_squared = get_particle_distance_squared(x, y, z, rx[i],ry[i],rz[i], l_x, l_y, l_z, half_l_x, half_l_y, half_l_z); + if dist_squared < cutoff_sqr { + let (e,v) = eval_pair_energy(dist_squared, e_shift); + energy += e; + } + } + return energy; +} + +fn main() { + + let mut trj_reader = TrjReader::new(&"2phases/10k_1step.xyz".to_string()); + let mut frame = trj_reader.next_frame(); + + 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; + + let mut gas_slab = 2.0; + let mut liquid_height = frame.box_z / (gas_slab + 1.0); + let liquid_volume = volume / (gas_slab + 1.0); + 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(); + + + let mut rng = rand::thread_rng(); + + 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 frame_count = 0; + let mut widom_sum_gas = 0.0; + let mut ideal_sum_gas = 0.0; + let mut widom_sum_liquid = 0.0; + let mut ideal_sum_liquid = 0.0; + + let mut avg_count = 0; + loop { + frame_count += 1; + avg_count += 1; + + + // test particle liquid energy + let lx = frame.box_x * rng.gen::(); + let ly = frame.box_y * rng.gen::(); + let lz = liquid_height * rng.gen::(); + 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(); + + // test particle gas energy + let gx = frame.box_x * rng.gen::(); + let gy = frame.box_y * rng.gen::(); + let gz = ((frame.box_z - liquid_height) * rng.gen::()) + liquid_height; + 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 number of particles in each phase + 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; } + else { gas_count += 1.0; } + } + + // calculate ideal gas potential for both phases + let wave = MKSA_PLANCKS_CONSTANT_H / (2.0 * std::f64::consts::PI * MASS * frame.temperature / frame.lj_eps).sqrt(); + ideal_sum_gas += frame.temperature / frame.lj_eps * (gas_volume/(wave.powi(3)* gas_count)).ln(); + ideal_sum_liquid += frame.temperature / frame.lj_eps * (liquid_volume/(wave.powi(3)* liquid_count)).ln(); + + + + if avg_count == 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; + println!("Frame {} gas {} liquid {} particlecount:{}/{}", frame_count, ideal_gas + excess_gas, ideal_liquid + excess_liquid, gas_count, liquid_count); + avg_count = 0; + ideal_sum_gas = 0.0; + ideal_sum_liquid = 0.0; + widom_sum_gas = 0.0; + widom_sum_liquid = 0.0; + + } + + // liquid_gas_sum = 0.0; + // excess_liquid_sum = 0.0; + // gas_gas_sum = 0.0; + // excess_gas_sum = 0.0; + // avg_count = 0; + + + + + if !trj_reader.update_with_next(&mut frame) { break } + } + + + + + +} \ No newline at end of file