diff --git a/Cargo.toml b/Cargo.toml index 19fb844..dfff524 100644 --- a/Cargo.toml +++ b/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" diff --git a/src/density_z.rs b/src/density_z.rs new file mode 100644 index 0000000..8ca4ecd --- /dev/null +++ b/src/density_z.rs @@ -0,0 +1,97 @@ +mod trajectory; +use trajectory::*; +use std::env; + +fn main() { + // open file + let args: Vec = 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::().unwrap(); + } else if args[i] == "--slabs" { + slabs = args[i+1].parse::().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 = 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)); +} \ No newline at end of file diff --git a/src/energy.rs b/src/energy.rs index bcb3491..c27394c 100644 --- a/src/energy.rs +++ b/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) { diff --git a/src/main.rs b/src/main.rs index 80f1534..15b0397 100644 --- a/src/main.rs +++ b/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::() - 0.5 ) * displacement; rz[rnd_index] += ( rng.gen::() - 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); diff --git a/src/trajectory.rs b/src/trajectory.rs index 709a322..27a0307 100644 --- a/src/trajectory.rs +++ b/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::().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; } + } + } + } diff --git a/src/widom.rs b/src/widom.rs index 85fb25d..60fc12d 100644 --- a/src/widom.rs +++ b/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 = 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::().unwrap(); + } else if args[i] == "-ls" { + liquid_start = args[i + 1].parse::().unwrap(); + } else if args[i] == "-le" { + liquid_end = 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) }; 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); + // rnd number generator for particle insertion + let mut rng = rand::thread_rng(); + + + // 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::(); let ly = frame.box_y * rng.gen::(); - let lz = liquid_height * rng.gen::(); - + let lz = liquid_start + (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(); +// widom_sum_liquid += widom_e_liquid; - // test particle gas energy + // gas test particle + let upper_or_lower = rng.gen::(); 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 gz; + if upper_or_lower { + gz = rng.gen::() * liquid_start; + } else { + gz = rng.gen::() * (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 - 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(); + // calculate ideal gas potential + ideal_sum_gas += -frame.temperature / frame.lj_eps * ((gas_volume/gas_count) * wave3).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;