diff --git a/src/main.rs b/src/main.rs index a363358..4401341 100644 --- a/src/main.rs +++ b/src/main.rs @@ -223,7 +223,8 @@ fn main() { 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; - println_stderr!("Minim {:<10} Energy: {:<30.3} Accept.: {:<4.1}% dr: {:.3}", step+1, energy, acceptance_rate, displacement); + let avg_energy = energy_sum / step_counter as f64; + println_stderr!("Minim {:<10} Energy: {:<30.3} Accept.: {:<4.1}% dr: {:.3}", step+1, avg_energy, acceptance_rate, displacement); if SCALE { let scale_factor = (TRIES_INTENDED/tries_per_step * DISP_SCALE_FACTOR).abs(); @@ -234,6 +235,7 @@ fn main() { } step_counter = 0; accept_counter = 0; + energy_sum = 0.0; } } diff --git a/src/trajectory.rs b/src/trajectory.rs index cb2a4ea..709a322 100644 --- a/src/trajectory.rs +++ b/src/trajectory.rs @@ -6,7 +6,7 @@ use std::error::Error; use std::io::prelude::*; use std::fs::File; use std::path::Path; - +use std::fmt; use std::io::BufReader; @@ -60,6 +60,21 @@ pub struct Frame { pub lj_cutoff: f64, } +impl fmt::Debug for Frame { + fn fmt(&self, fmt: &mut fmt::Formatter) -> fmt::Result { + fmt.debug_struct("Frame") + .field("num_particles", &self.num_particles) + .field("box_x", &self.box_x) + .field("box_y", &self.box_y) + .field("box_z", &self.box_z) + .field("temperature", &self.temperature) + .field("lj_eps", &self.lj_eps) + .field("lj_sig", &self.lj_sig) + .field("lj_cutoff", &self.lj_cutoff) + .finish() + } +} + pub struct TrjReader { pub reader: BufReader, } diff --git a/src/widom.rs b/src/widom.rs index 22e91c2..85fb25d 100644 --- a/src/widom.rs +++ b/src/widom.rs @@ -4,6 +4,7 @@ mod energy; use energy::*; extern crate rand; use rand::Rng; +use std::env; const LJ_EPS : f64 = 1.0; const LJ_SIG : f64 = 1.0; @@ -11,22 +12,104 @@ 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 INSERTIONS_PER_STEP : usize = 1000; +const RUN_AVG_SIZE : usize = 10; 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; } +fn main() { + // open file + let args: Vec = env::args().collect(); + let mut trj_reader = TrjReader::new(&args[1]); + + 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 + let gas_slab = args[2].parse::().unwrap(); + let 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 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; + 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; + + // loop over frames + let mut avg_count = 0; + loop { + frame_count += 1; + + // 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; } + } + + // test particle insertion multiple times + for i in 0..INSERTIONS_PER_STEP { + avg_count += 1; + + // liquid test partciles + 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 ideal gas potential for both phases + 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(); + + } + + // print running averages + 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; + 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; + ideal_sum_gas = 0.0; + ideal_sum_liquid = 0.0; + widom_sum_gas = 0.0; + widom_sum_liquid = 0.0; + } + + if !trj_reader.update_with_next(&mut frame) { + break; + } } - 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; @@ -36,103 +119,8 @@ fn get_particle_insertion_energy(rx: &[f64], ry: &[f64], rz: &[f64], num_particl 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; + energy += eval_pair_energy(dist_squared, e_shift).0; } } 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