From 60733f3531e372c92cffb599da163ec848960364 Mon Sep 17 00:00:00 2001 From: Daniel Bauer Date: Sat, 4 Feb 2017 00:48:22 +0100 Subject: [PATCH] comments --- src/energy.rs | 28 ++- src/main.rs | 537 +++++++++++++++++++++-------------------- src/surface_tension.rs | 50 ++-- src/trajectory.rs | 4 +- src/widom.rs | 75 +++--- 5 files changed, 365 insertions(+), 329 deletions(-) diff --git a/src/energy.rs b/src/energy.rs index 360f77c..44a897d 100644 --- a/src/energy.rs +++ b/src/energy.rs @@ -1,4 +1,7 @@ + #![allow(dead_code)] +/// Calculates the total energy and virial of a system containing num_particles with coords rx,ry,rz +/// of size l_x, l_y, l_z and given cutoff + corrections pub fn get_total_energy(rx: &[f64], ry: &[f64], rz: &[f64], num_particles: usize, l_x: f64, l_y: f64, l_z: f64, cutoff_squared: f64, e_corr: f64, e_shift: f64) -> (f64, f64) { let mut energy = 0.0; let mut virial = 0.0; @@ -19,6 +22,8 @@ pub fn get_total_energy(rx: &[f64], ry: &[f64], rz: &[f64], num_particles: usize return (energy, virial); } +/// Calculates the particle energy and virial for particle at p_index in system containing num_particles with coords rx,ry,rz +/// of size l_x, l_y, l_z and given cutoff + corrections pub fn get_particle_energy(rx: &[f64], ry: &[f64], rz: &[f64], p_index: usize, num_particles: usize, l_x: f64, l_y: f64, l_z: f64, cutoff_squared: f64, e_shift: f64) -> (f64, f64) { let mut energy = 0.0; let mut virial = 0.0; @@ -38,6 +43,7 @@ pub fn get_particle_energy(rx: &[f64], ry: &[f64], rz: &[f64], p_index: usize, n return (energy, virial); } +// squared distance between 2 particles regarding the minimum image convention 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).abs(); let mut dy = (y1 - y2).abs(); @@ -50,13 +56,17 @@ pub fn get_particle_distance_squared(x1: f64,y1: f64,z1: f64,x2: f64,y2: f64,z2: if dz > hl_z { dz -= l_y } else if dz < -hl_z { dz += l_z} 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(); +// one dimensional distance with applied minimum image convention +pub fn get_distance_with_pbc(x1: f64, x2: f64, length: f64, half_length: f64) -> f64 { + let mut d = (x1-x2).abs(); + if d > half_length { d -= length } + else if d < -half_length { d += length } + return d; } + #[test] fn test_get_particle_distance_squared() { let (x1, y1, z1) = (0.0, 0.0, 0.0); @@ -81,6 +91,7 @@ fn test_get_particle_distance_squared() { assert!(dist - 12.0 < 0.00001); } +/// calculate the lj energy and virial between two particles from given square distance 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; @@ -101,11 +112,12 @@ fn test_eval_pair_energy() { assert!( (e - 224.0).abs() < 0.00001, "{}", e); } - -pub fn eval_virial(distance: f64, LJ_EPS: f64, LJ_SIG: 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 ); +/// calculate the virial between two particles from given square distance. If energy is required too, +/// see eval_pair_energy which does energy and virial +pub fn eval_virial(distance: f64, lj_eps: f64, lj_sig: 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] diff --git a/src/main.rs b/src/main.rs index 2d22452..415912e 100644 --- a/src/main.rs +++ b/src/main.rs @@ -11,13 +11,19 @@ use argparse::{ArgumentParser, Store, StoreFalse, StoreTrue}; mod trajectory; use trajectory::*; - +// LJ params const LJ_EPS : f64 = 1.0; const LJ_SIG : f64 = 1.0; - + // intended acceptance rate = 33% const TRIES_INTENDED : f64 = 3.0; + +// factor for auto displacement scaling const DISP_SCALE_FACTOR : f64 = 0.1; +const SCALE_INTERVAL : usize = 5000; + +const EQUILIBRATION_OUTPUT_INTERVAL : usize = 5000; +const SAMPLING_OUTPUT_INTERVAL : usize = 5000; // easy printing to stderr macro_rules! println_stderr( @@ -27,7 +33,278 @@ macro_rules! println_stderr( } } ); -fn parse_cmd_args(NUM_STEPS: &mut usize, NUM_MINIM_STEPS: &mut usize, +fn main() { + println_stderr!(""); + println_stderr!("################################################################"); + println_stderr!("################## LJ Monte Carlo Simulation #################"); + println_stderr!("################################################################"); + println_stderr!(""); + + /** Definition of default run parameters **/ + let mut eq_steps = 1000000; + let mut sample_steps = 100000; + + let mut num_particles: usize = 512; + let mut density = 0.7; + let mut temperature = 0.9; + + let mut cutoff = 3.0; + + let mut TAILCORR : bool = true; + let mut SHIFT: bool = true; + + let mut displacement = 0.1; // max particle displacement in one dimension + let mut SCALE: bool = true; // switch for displacement scaling + + // scale factor in z for vaccuum space + let mut vacuum_slab = 0.0; + + // output config + let mut output_prefix = "montecarlo".to_string(); // .xyz will be append + let mut output_interval : i64 = 100; + let mut output_minim : bool = false; + + // parse cmd line arguments and override defaults + parse_cmd_args(&mut sample_steps, &mut eq_steps, &mut num_particles, + &mut density, &mut temperature, + &mut cutoff, &mut displacement, &mut SCALE, &mut TAILCORR, &mut SHIFT, + &mut output_prefix, &mut output_interval, &mut output_minim, + &mut vacuum_slab); + + + /** Initialize the system **/ + let beta = 1.0/temperature; + + let mut volume = (num_particles as f64)/ density; + let length = volume.cbrt(); + let (l_x, l_y, mut l_z) = (length, length, length); + + let cutoff_squared = cutoff * cutoff; + let max_displacement = length / 2.0; // displacement wont be scaled over that + + // initialize randomness - TODO seed? + let mut rng = rand::thread_rng(); + let particle_range = Range::new(0, num_particles); + + // randomly place particles in the box + let mut rx : Vec = vec![]; + let mut ry : Vec = vec![]; + let mut rz : Vec = vec![]; + loop { + rx.push(l_x * rng.gen::()); + ry.push(l_y * rng.gen::()); + rz.push(l_z * rng.gen::()); + if rx.len() == num_particles { break; } + } + + // 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; + } + } + + // calculation of shift and tailcorrections + let e_shift = if SHIFT { 4.0 * LJ_EPS * ( (LJ_SIG/cutoff).powi(12) - (LJ_SIG/cutoff).powi(6) ) } else { 0.0 }; + let e_corr = if TAILCORR { 8.0/3.0*std::f64::consts::PI*density*LJ_EPS*LJ_SIG.powi(3)*((1.0/3.0*(LJ_SIG/cutoff).powi(9)) - (LJ_SIG/cutoff).powi(3)) } else { 0.0 }; + let p_corr = if TAILCORR { 16.0/3.0*std::f64::consts::PI*density.powi(2)*LJ_EPS*LJ_SIG.powi(3)*((2.0/3.0*(LJ_SIG/cutoff).powi(9)) - (LJ_SIG/cutoff).powi(3)) } else { 0.0 }; + + println_stderr!("Particles: {}, Density: {}, Temperature: {}", num_particles, density, temperature); + println_stderr!("System volume: {:8.3}, Dimensions {:.3}/{:.3}/{:.3}", volume, l_x, l_y, l_z); + println_stderr!("Minimization steps: {}, Sampling steps: {}", eq_steps, sample_steps); + println_stderr!("LJ params eps: {}, sigma: {}, cutoff: {}", LJ_EPS, LJ_SIG, cutoff); + println_stderr!("Tailcorr: {:8.3}, Shift: {:8.3}, Pressurecorr: {:8.3}", e_corr, e_shift, p_corr); + + // energy and average sums + let (mut energy, mut virial) = get_total_energy(&rx, &ry, &rz, num_particles, l_x, l_y, l_z, cutoff_squared, e_corr, e_shift); + let mut energy_sum = 0.0; + let mut virial_sum = 0.0; + let mut step_counter = 0; + let mut accept_counter = 0; + + + // prepare and write first trajectory frame + let mut trajectory : XYZTrajectory = XYZTrajectory::new(&format!("{}.xyz", output_prefix)); + if output_minim { trajectory.write(&rx, &ry, &rz, num_particles, l_x, l_y, l_z, temperature, LJ_EPS, LJ_SIG, cutoff, true); } + + + println_stderr!(""); + println_stderr!("################################################################"); + println_stderr!("######################## Equilibration #######################"); + println_stderr!("################################################################"); + println_stderr!(""); + + + // START OF METROPOLIS + /*****************************************************************************************/ + + for step in 0..eq_steps+sample_steps { + + // select rnd particle + let rnd_index = particle_range.ind_sample(&mut rng); + + // store old position + let oldX = rx[rnd_index]; + let oldY = ry[rnd_index]; + let oldZ = rz[rnd_index]; + + // old particle energy + let (old_particle_energy, old_particle_virial) = get_particle_energy(&rx, &ry, &rz, rnd_index, num_particles, l_x, l_y, l_z, cutoff_squared, e_shift); + + // rnd displacement and PBC + rx[rnd_index] += ( rng.gen::() - 0.5 ) * displacement; + 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 ry[rnd_index] < 0.0 { 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 } + + // 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); + + let dE = new_particle_energy - old_particle_energy; + + // acceptance rule + if dE < 0.0 || rng.gen::() < (-beta * dE).exp() { + accept_counter += 1; + energy += dE; + virial += new_particle_virial - old_particle_virial; + + // recalculate total energy every 1000 steps to account for rounding errors in particle energy function + if step % 10000 == 0 { + let (e, v) = get_total_energy(&rx, &ry, &rz, num_particles, l_x, l_y, l_z, cutoff_squared, e_corr, e_shift); + energy = e; + virial = v; + } + } else { + // restore old positions if move is rejected + rx[rnd_index] = oldX; + ry[rnd_index] = oldY; + rz[rnd_index] = oldZ; + } + + // update average sums + step_counter += 1; + energy_sum += energy; + virial_sum += virial; + + // reset average sums for sampling + if step == eq_steps-1 { + println_stderr!(""); + println_stderr!("################################################################"); + println_stderr!("########################## Sampling ##########################"); + println_stderr!("################################################################"); + println_stderr!(""); + step_counter = 0; + accept_counter = 0; + energy_sum = 0.0; + virial_sum = 0.0; + } + + // Everything below here is not part of the metropolis sampling (extras) + + // displacement scaling during equilibration for good acceptance ratios + if SCALE && step < eq_steps && step % SCALE_INTERVAL == 0 { + let tries_per_step : f64 = step_counter as f64 /accept_counter as f64; + + // will increase the max displacement if the acceptance rate is too high and vice versa + let scale_factor = (TRIES_INTENDED/tries_per_step * DISP_SCALE_FACTOR).abs(); + if tries_per_step < TRIES_INTENDED - 0.2 && displacement < max_displacement { + displacement += displacement * scale_factor; + } else if tries_per_step > TRIES_INTENDED + 0.2 && displacement > 0.0 { + displacement -= displacement * scale_factor; + } + step_counter = 0; + accept_counter = 0; + energy_sum = 0.0; + } + + // print some output during equilibration + if step < eq_steps && step_counter % EQUILIBRATION_OUTPUT_INTERVAL == 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; + let avg_energy = energy_sum / step_counter as f64; + let avg_virial = virial_sum / step_counter as f64; + println_stderr!("Eq {:<10} Energy: {:<30.3} Virial: {:<30.3} Accept.: {:<4.1}% dr: {:.3}", step+1, avg_energy, avg_virial, acceptance_rate, displacement); + } + + // print some output during sampling + if step > eq_steps && step_counter % SAMPLING_OUTPUT_INTERVAL == 0 { + println_stderr!("Step {:<10} Energy: {:<30.3}", step_counter, energy); + } + + // write trajectory + if step as i64 % output_interval == 0 { + if step > eq_steps || output_minim { + trajectory.write(&rx, &ry, &rz, num_particles, l_x, l_y, l_z, temperature, LJ_EPS, LJ_SIG, cutoff, true); + } + } + } + + // END OF METROPOLIS + /*****************************************************************************************/ + println_stderr!("Done sampling!"); + + let final_energy = energy_sum/step_counter as f64; + let particle_energy = final_energy / num_particles as f64; + let final_virial = virial_sum / 3.0 / step_counter as f64 / volume; + let pressure = virial_sum / 3.0 / step_counter as f64 / volume + density * temperature + p_corr; + let final_acceptance_rate = 1.0/((accept_counter as f64)/(step_counter as f64)) * 100.0; + + println_stderr!(""); + println_stderr!("################################################################"); + println_stderr!("########################## Results ###########################"); + println_stderr!("################################################################"); + println_stderr!(""); + println!( +"Minimization: {} +Steps: {} + +# Lennard Jones Params +epsilon: {} +sigma: {} +cutoff: {} + +# System +Particles: {} +Density: {} +Temperature: {} +Volume: {} +Box dimension: {:.3}/{:.3}/{:.3} +Max Displacement: {} + +# Correction +Energy correction: {} +Shift: {} +P-Correction: {} + +# Averages +Tries: {} +Accepted: {} +Acceptance: {:.2}% +Energy: {} +Energy per particle: {} +Virial: {} +Pressure: {}", + eq_steps, sample_steps, + LJ_EPS, LJ_SIG, cutoff, + num_particles, density, temperature, volume, l_x, l_y, l_z, displacement, + e_corr, e_shift, p_corr, + step_counter, accept_counter, final_acceptance_rate, final_energy, particle_energy, final_virial, pressure); + + trajectory.write(&rx, &ry, &rz, num_particles, l_x, l_y, l_z, temperature, LJ_EPS, LJ_SIG, cutoff, true); +} + +// Parse command line arguments +fn parse_cmd_args(NUM_STEPS: &mut usize, NUM_eq_steps: &mut usize, NUM_PARTICLES: &mut usize, DENSITY: &mut f64, TEMPERATURE: &mut f64, CUTOFF: &mut f64, MAX_DISP_START: &mut f64, SCALE: &mut bool, TAILCORR: &mut bool, SHIFT: &mut bool, OUTPUT_PREFIX: &mut String, OUTPUT_INTERVAL: &mut i64, OUTPUT_MINIM: &mut bool, @@ -37,7 +314,7 @@ fn parse_cmd_args(NUM_STEPS: &mut usize, NUM_MINIM_STEPS: &mut usize, ap.refer(NUM_STEPS) .add_option(&["-n", "--nsteps"], Store, "Simulation steps: Number of steps for averaging" ); - ap.refer(NUM_MINIM_STEPS) + ap.refer(NUM_eq_steps) .add_option(&["-m", "--nminimsteps"], Store, "Minimization steps: Number of steps before averaging starts"); ap.refer(NUM_PARTICLES) @@ -78,255 +355,3 @@ fn parse_cmd_args(NUM_STEPS: &mut usize, NUM_MINIM_STEPS: &mut usize, "Disable lj shifting"); ap.parse_args_or_exit(); } - -fn main() { - - // define all the stuff - let mut minim_steps = 1000000; - let mut sample_steps = 100000; - - let mut num_particles: usize = 512; - let mut density = 0.7; - let mut temperature = 0.9; - - let mut cutoff = 3.0; - let mut displacement = 0.1; - - let mut TAILCORR : bool = true; - let mut SHIFT: bool = true; - let mut SCALE: bool = true; - - let mut vacuum_slab = 0.0; - - let mut output_prefix = "montecarlo".to_string(); - let mut output_interval : i64 = 100; - let mut output_minim : bool = false; - parse_cmd_args(&mut sample_steps, &mut minim_steps, &mut num_particles, - &mut density, &mut temperature, - &mut cutoff, &mut displacement, &mut SCALE, &mut TAILCORR, &mut SHIFT, - &mut output_prefix, &mut output_interval, &mut output_minim, - &mut vacuum_slab); - - println_stderr!(""); - println_stderr!("################################################################"); - println_stderr!("################## LJ Monte Carlo Simulation #################"); - println_stderr!("################################################################"); - println_stderr!(""); - - - // initialize stuff - let beta = 1.0/temperature; - let mut volume = (num_particles as f64)/ density; - let length = volume.cbrt(); - let (l_x, l_y, mut l_z) = (length, length, length); - - let cutoff_squared = cutoff * cutoff; - let max_displacement = length / 2.0; - - let mut rng = rand::thread_rng(); - let particle_range = Range::new(0, num_particles); - - let mut rx : Vec = vec![]; - let mut ry : Vec = vec![]; - let mut rz : Vec = vec![]; - loop { - rx.push(l_x * rng.gen::()); - ry.push(l_y * rng.gen::()); - rz.push(l_z * rng.gen::()); - if rx.len() == num_particles { break; } - } - - // 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 }; - let e_corr = if TAILCORR { 8.0/3.0*std::f64::consts::PI*density*LJ_EPS*LJ_SIG.powi(3)*((1.0/3.0*(LJ_SIG/cutoff).powi(9)) - (LJ_SIG/cutoff).powi(3)) } else { 0.0 }; - let p_corr = if TAILCORR { 16.0/3.0*std::f64::consts::PI*density.powi(2)*LJ_EPS*LJ_SIG.powi(3)*((2.0/3.0*(LJ_SIG/cutoff).powi(9)) - (LJ_SIG/cutoff).powi(3)) } else { 0.0 }; - - println_stderr!("Particles: {}, Density: {}, Temperature: {}", num_particles, density, temperature); - println_stderr!("System volume: {:8.3}, Dimensions {:.3}/{:.3}/{:.3}", volume, l_x, l_y, l_z); - println_stderr!("Minimization steps: {}, Sampling steps: {}", minim_steps, sample_steps); - println_stderr!("LJ params eps: {}, sigma: {}, cutoff: {}", LJ_EPS, LJ_SIG, cutoff); - println_stderr!("Tailcorr: {:8.3}, Shift: {:8.3}, Pressurecorr: {:8.3}", e_corr, e_shift, p_corr); - - let (mut energy, mut virial) = get_total_energy(&rx, &ry, &rz, num_particles, l_x, l_y, l_z, cutoff_squared, e_corr, e_shift); - let mut energy_sum = 0.0; - let mut virial_sum = 0.0; - let mut step_counter = 0; - let mut accept_counter = 0; - - println_stderr!(""); - println_stderr!("################################################################"); - println_stderr!("##################### Energy Minimization ####################"); - println_stderr!("################################################################"); - println_stderr!(""); - - // prepare and write first trajectory frame - let mut trajectory : XYZTrajectory = XYZTrajectory::new(&format!("{}.xyz", output_prefix)); - if output_minim { trajectory.write(&rx, &ry, &rz, num_particles, l_x, l_y, l_z, temperature, LJ_EPS, LJ_SIG, cutoff, true); } - -// let vacuum_scale_step = (minim_steps as f64 * 0.1) as usize; - for step in 0..minim_steps+sample_steps { - - // first minimizate solvent phase, then add vacuum slab -// if vacuum_slab > 0.0 && step == vacuum_scale_step { // increase space in z -// let scale = vacuum_slab + 1.0; -// l_z *= scale; -// volume *= scale; -// density /= scale; -// } - // select rnd particle - let rnd_index = particle_range.ind_sample(&mut rng); - - // store old position - let oldX = rx[rnd_index]; - let oldY = ry[rnd_index]; - let oldZ = rz[rnd_index]; - - // old particle energy - let (old_particle_energy, old_particle_virial) = get_particle_energy(&rx, &ry, &rz, rnd_index, num_particles, l_x, l_y, l_z, cutoff_squared, e_shift); - - // rnd displacement and PBC - rx[rnd_index] += ( rng.gen::() - 0.5 ) * displacement; - 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 ry[rnd_index] < 0.0 { 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 } - - // 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); - let dE = new_particle_energy - old_particle_energy; - - //accept move - if rng.gen::() < (-beta * dE).exp() { - accept_counter += 1; - if step % 1000 == 0 { // calculate total energy every 1000 steps to account for rounding errors - let (e, v) = get_total_energy(&rx, &ry, &rz, num_particles, l_x, l_y, l_z, cutoff_squared, e_corr, e_shift); - energy = e; - virial = v; - } else { - energy += dE; - virial += new_particle_virial - old_particle_virial; - } - } else { // or restore old position - rx[rnd_index] = oldX; - ry[rnd_index] = oldY; - rz[rnd_index] = oldZ; - } - - // update sums for averaging - step_counter += 1; - energy_sum += energy; - virial_sum += virial; - - // print some output during minimization - 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; - let avg_energy = energy_sum / step_counter as f64; - let avg_virial = virial_sum / step_counter as f64; - println_stderr!("Minim {:<10} Energy: {:<30.3} Virial: {:<30.3} Accept.: {:<4.1}% dr: {:.3}", step+1, avg_energy, avg_virial, acceptance_rate, displacement); - - if SCALE { - let scale_factor = (TRIES_INTENDED/tries_per_step * DISP_SCALE_FACTOR).abs(); - if tries_per_step < TRIES_INTENDED - 0.2 && displacement < max_displacement { - displacement += displacement * scale_factor; - } else if tries_per_step > TRIES_INTENDED + 0.2 && displacement > 0.0 { - displacement -= displacement * scale_factor; - } - step_counter = 0; - accept_counter = 0; - energy_sum = 0.0; - } - } - - // reset sums for sampling - if step == minim_steps-1 { - println_stderr!(""); - println_stderr!("################################################################"); - println_stderr!("########################## Sampling ##########################"); - println_stderr!("################################################################"); - println_stderr!(""); - step_counter = 0; - accept_counter = 0; - energy_sum = 0.0; - virial_sum = 0.0; - } - - if step > minim_steps && step_counter % 5000 == 0 { - println_stderr!("Step {:<10} Energy: {:<30.3}", step_counter, energy); - } - - // write trajectory maybe - if step as i64 % output_interval == 0 { - if step > minim_steps || output_minim { - trajectory.write(&rx, &ry, &rz, num_particles, l_x, l_y, l_z, temperature, LJ_EPS, LJ_SIG, cutoff, true); - } - } - - } - - - let final_energy = energy_sum/step_counter as f64; - let particle_energy = final_energy / num_particles as f64; - let final_virial = virial_sum / 3.0 / step_counter as f64 / volume; - let pressure = virial_sum / 3.0 / step_counter as f64 / volume + density * temperature + p_corr; - let final_acceptance_rate = 1.0/((accept_counter as f64)/(step_counter as f64)) * 100.0; - - println_stderr!("Done sampling!"); - println_stderr!(""); - println_stderr!("################################################################"); - println_stderr!("########################## Results ###########################"); - println_stderr!("################################################################"); - println_stderr!(""); - println!( -"Minimization: {} -Steps: {} - -# Lennard Jones Params -epsilon: {} -sigma: {} -cutoff: {} - -# System -Particles: {} -Density: {} -Temperature: {} -Volume: {} -Box dimension: {:.3}/{:.3}/{:.3} -Max Displacement: {} - -# Correction -Energy correction: {} -Shift: {} -P-Correction: {} - -# Averages -Tries: {} -Accepted: {} -Acceptance: {:.2}% -Energy: {} -Energy per particle: {} -Virial: {} -Pressure: {}", - minim_steps, sample_steps, - LJ_EPS, LJ_SIG, cutoff, - num_particles, density, temperature, volume, l_x, l_y, l_z, displacement, - e_corr, e_shift, p_corr, - step_counter, accept_counter, final_acceptance_rate, final_energy, particle_energy, final_virial, pressure); - - trajectory.write(&rx, &ry, &rz, num_particles, l_x, l_y, l_z, temperature, LJ_EPS, LJ_SIG, cutoff, true); -} diff --git a/src/surface_tension.rs b/src/surface_tension.rs index ccb7890..fbdba5b 100644 --- a/src/surface_tension.rs +++ b/src/surface_tension.rs @@ -7,32 +7,14 @@ use std::env; const LJ_EPS : f64 = 1.0; const LJ_SIG : f64 = 1.0; - -fn get_distance_with_pbc(x1: f64, x2: f64, length: f64, half_length: f64) -> f64 { - 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 ); -} - +const AVG_OUTPUT_INTERVAL : usize = 10; fn main() { - - // parse args - let args: Vec = env::args().collect(); let mut filename = "montecarlo.xyz".to_string(); let mut skip: usize = 0; + + // parse cmd line args + let args: Vec = env::args().collect(); for i in 0..args.len() { if args[i] == "-f" { filename = args[i + 1].clone(); @@ -41,17 +23,16 @@ fn main() { } } - // open file and skip to requiested position + // open file and skip to requested position let mut trj_reader = TrjReader::new(&filename); if skip > 0 { trj_reader.skip(skip) }; - // trajectory information + // read first trajectory and system params 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; let num_particles = frame.num_particles; + println!("{:?}", frame); let box_half_x = frame.box_x / 2.0; let box_half_y = frame.box_y / 2.0; @@ -63,6 +44,7 @@ fn main() { let variable_without_name = frame.temperature/LJ_EPS * density; + println!("Calculating surface tension"); println!("~~~ THIS IS A RUNNING AVERAGE! ~~~"); loop { frame_count += 1; @@ -71,6 +53,7 @@ fn main() { let mut trace_z = 0.0; for i in 0..num_particles { for j in i+1..num_particles { + // this needs some optimization for speed let dist_sqrt = get_particle_distance_squared(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 dist = dist_sqrt.sqrt(); let dx = get_distance_with_pbc(frame.rx[i], frame.rx[j], frame.box_x, box_half_x); @@ -88,7 +71,7 @@ fn main() { p_z_sum += p_zz; /////////////////////////////////// - if frame_count % 10 == 0 { + if frame_count % AVG_OUTPUT_INTERVAL == 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; @@ -100,7 +83,20 @@ fn main() { // p_xy_sum = 0.0; } + // read next frame if !trj_reader.update_with_next(&mut frame) { break } } } + +/// calc surface tension from box z size and pressure tensor +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 ); +} diff --git a/src/trajectory.rs b/src/trajectory.rs index 27a0307..1ff80a8 100644 --- a/src/trajectory.rs +++ b/src/trajectory.rs @@ -1,5 +1,5 @@ #![allow(dead_code)] -#![allow(unused_must_use)] // hate. +#![allow(unused_must_use)] #![allow(unused_variables)] use std::error::Error; @@ -9,7 +9,6 @@ use std::path::Path; use std::fmt; use std::io::BufReader; - pub struct XYZTrajectory { file: File, @@ -31,7 +30,6 @@ impl XYZTrajectory { } - pub fn write(&mut self, rx: &[f64], ry: &[f64], rz: &[f64], num_particles: usize, box_x : f64, box_y : f64, box_z : f64, temp: f64 ,lj_eps : f64, lj_sig : f64, lj_cutoff : f64, flush: bool) { self.file.write(format!("{} ## Box: {} {} {} Temp: {} LJ: {}/{}/{}\n", num_particles, box_x,box_y,box_z,temp, lj_eps, lj_sig, lj_cutoff).as_bytes()); for i in 0..num_particles { diff --git a/src/widom.rs b/src/widom.rs index ac60e86..3abc3d8 100644 --- a/src/widom.rs +++ b/src/widom.rs @@ -1,3 +1,5 @@ +#![allow(unused_variables)] + mod trajectory; use trajectory::*; mod energy; @@ -14,30 +16,19 @@ static MASS : f64 = 1.0; const RUN_AVG_SIZE : usize = 1; -const SHIFT : bool = true; - -fn eval_ideal_potential(temperature: f64, lj_eps: f64, volume: f64, particles: f64, thermal_wavelength3: f64) -> f64 { - let density = volume/particles; - return -temperature/lj_eps * ( density*thermal_wavelength3 ).ln(); -} - -#[test] -fn test_eval_ideal_potential() { - let expected = 12.476649; - let result = eval_ideal_potential(2.0,1.0,10.0,512.0,0.1); - assert!((expected-result).abs() < 0.00001, "{}", result); -} - fn main() { - // parse args - let args: Vec = env::args().collect(); + // default values let mut filename = "montecarlo.xyz".to_string(); - let mut skip: usize = 0; - let mut insertions: usize = 100; + let mut skip: usize = 0; // skip some frames + let mut insertions: usize = 100; // number of particle insertions per step per phase // liquid phase boundaries let mut liquid_start = 0.0; let mut liquid_end = 0.0; + let mut shift = true; + + // parse command line arguments + let args: Vec = env::args().collect(); for i in 0..args.len() { if args[i] == "-f" { filename = args[i + 1].clone(); @@ -49,18 +40,19 @@ fn main() { liquid_end = args[i + 1].parse::().unwrap(); } else if args[i] == "-n" { insertions = args[i+1].parse::().unwrap(); + } else if args[i] == "--noshift" { + shift = false } } - // open file and skip to requiested position + // open file and skip to requested position let mut trj_reader = TrjReader::new(&filename); if skip > 0 { trj_reader.skip(skip) }; - // trajectory information + // Get first frame and read system configuration 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; @@ -77,7 +69,7 @@ fn main() { let tw3 = ((2.0 * std::f64::consts::PI * MASS * frame.temperature / frame.lj_eps)/MKSA_PLANCKS_CONSTANT_H.powi(2)).powf(3.0/2.0); // LJ shift - 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 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 }; // average counters @@ -88,8 +80,6 @@ fn main() { let mut ideal_pot_liquid_sum = 0.0; let mut avg_count = 0; - // loop over all frames - let mut counter = 0; loop { frame_count += 1; @@ -116,14 +106,13 @@ fn main() { let gx = frame.box_x * rng.gen::(); let gy = frame.box_y * rng.gen::(); let mut gz = frame.box_z * rng.gen::(); - while gz > liquid_start && gz < liquid_end { // retry until we have a particle in gas - gz= frame.box_z * rng.gen::(); + while gz > liquid_start && gz < liquid_end { // retry until we hit the gas phase + gz= frame.box_z * rng.gen::(); } - 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 potentials + // ideal gas potentials ideal_pot_gas_sum += eval_ideal_potential(frame.temperature, frame.lj_eps, gas_volume, gas_count, tw3); ideal_pot_liquid_sum += eval_ideal_potential(frame.temperature, frame.lj_eps, liquid_volume, liquid_count, tw3); @@ -133,12 +122,14 @@ fn main() { if avg_count / insertions > RUN_AVG_SIZE { let ideal_gas_potential = ideal_pot_gas_sum / avg_count as f64; let ideal_liquid_potential = ideal_pot_liquid_sum / 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 mut gas_total = ideal_gas_potential + excess_gas; - let liquid_total = ideal_liquid_potential + excess_liquid; - println!("Frame {}\tg_ex: {:5}\tl_ex: {:5}\tg_tot: {:5}\tl_tot: {:5}\t\tnparticles: {}/{}", frame_count, excess_gas, excess_liquid, if gas_total.is_infinite() { 0.0 } else { 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); + let excess_gas_potential = -(widom_sum_gas/avg_count as f64).ln()/beta; + let excess_liquid_potential = -(widom_sum_liquid/avg_count as f64).ln()/beta; + let gas_total = ideal_gas_potential + excess_gas_potential; + let liquid_total = ideal_liquid_potential + excess_liquid_potential; + println!("Frame {}\tg_ex: {:5}\tl_ex: {:5}\tg_tot: {:5}\tl_tot: {:5}\t\tnparticles: {}/{}", + frame_count, excess_gas_potential, excess_liquid_potential, + if gas_total.is_infinite() { 0.0 } else { gas_total }, + liquid_total, gas_count, liquid_count); // reset averages for next round avg_count = 0; @@ -148,6 +139,7 @@ fn main() { widom_sum_liquid = 0.0; } + // jump to next frame if !trj_reader.update_with_next(&mut frame) { break; } @@ -155,7 +147,7 @@ fn main() { } - +/// calculates the energy for a hypothetic particle inserted at x,y,z 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; @@ -170,3 +162,16 @@ fn get_particle_insertion_energy(rx: &[f64], ry: &[f64], rz: &[f64], num_particl } return energy; } + +/// calculates the ideal gas chemical potential +fn eval_ideal_potential(temperature: f64, lj_eps: f64, volume: f64, particles: f64, thermal_wavelength3: f64) -> f64 { + let density = volume/particles; + return -temperature/lj_eps * ( density*thermal_wavelength3 ).ln(); +} + +#[test] +fn test_eval_ideal_potential() { + let expected = 12.476649; + let result = eval_ideal_potential(2.0,1.0,10.0,512.0,0.1); + assert!((expected-result).abs() < 0.00001, "{}", result); +}