From bc455738f47dcb2e4754eaea9a1daf89b0720391 Mon Sep 17 00:00:00 2001 From: danijoo Date: Sat, 21 Jan 2017 14:36:54 +0100 Subject: [PATCH] support of differnt box dimensions --- src/energy.rs | 16 +++++++--------- src/main.rs | 15 +++++++-------- 2 files changed, 14 insertions(+), 17 deletions(-) diff --git a/src/energy.rs b/src/energy.rs index bd43d29..68f796e 100644 --- a/src/energy.rs +++ b/src/energy.rs @@ -44,11 +44,11 @@ fn get_particle_distance_squared(x1: f64,y1: f64,z1: f64,x2: f64,y2: f64,z2: f64 let mut dz = z1 - z2; if dx > hl_x { dx -= l_x } - else if dx < -hl_x{ dx += -l_x} - if dy > hl_y { dy -= l_y } - else if dy < -hl_y{ dy += -l_y} - if dz > hl_z { dz -= l_z } - else if dz < -hl_z{ dz += -l_z} + else if dx < -hl_x { dx += l_x } + if dy > hl_y { dy -= l_y} + else if dy < -hl_y{ dy += l_y } + if dz > hl_z { dz -= l_y } + else if dz < -hl_z { dz += l_z} return dx*dx + dy*dy + dz*dz; } @@ -59,12 +59,10 @@ fn test_get_particle_distance_squared() { let (x2, y2, z2) = (5.0, 0.0, 0.0); // no pbc - let dist = get_particle_distance_squared(x1,y1,y1,y2,z2,z2, 20.0, 20.0, 20.0, 10.0, 10.0, 10.0); - assert!( (dist - 25.0) < 0.00001, "{}", dist); + assert!( (get_particle_distance_squared(x1,y1,z1,x2,y2,z2, 20.0, 20.0, 20.0, 10.0, 10.0, 10.0) - 25.0) < 0.00001); // with pbc - let dist = get_particle_distance_squared(x1,y1,y1,y2,z1,z2, 9.0, 9.0, 9.0, 4.5, 4.5, 4.5); - assert!( (dist - 16.0) < 0.00001, "{}", dist); + 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) -> (f64, f64) { diff --git a/src/main.rs b/src/main.rs index 30c03c3..003b031 100644 --- a/src/main.rs +++ b/src/main.rs @@ -24,8 +24,8 @@ macro_rules! println_stderr( fn main() { // define all the stuff - let sample_steps = 1000000; - let minim_steps = 1000000; + let minim_steps = 1000000; + let sample_steps = 100000; let num_particles: usize = 512; let density = 0.7; @@ -45,9 +45,8 @@ fn main() { // initialize stuff let beta = 1.0/temperature; let volume = (num_particles as f64)/ density; - let l_x = volume.cbrt(); - let l_y = l_x; - let l_z = l_x; + let length = volume.cbrt(); + let (l_x, l_y, l_z) = (length, length, length); let cutoff_squared = cutoff * cutoff; let mut rng = rand::thread_rng(); @@ -67,7 +66,7 @@ fn main() { 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!("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}, Pressurecprr: {:8.3}", e_corr, SHIFT, p_corr); @@ -135,7 +134,7 @@ fn main() { virial_sum += virial; - if step_counter % 5000 == 0 && step < minim_steps { + if step_counter % 5000 == 0 { println!("Minim {}\tEnergy: {:.3}\tVirial: {:.3}\tAcceptance:{:.1}\tDisplacement: {:.3}", step_counter, energy, virial, 666, displacement); } @@ -152,7 +151,7 @@ fn main() { 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 / num_particles as f64 / volume; + 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; println!("Steps: {}", step_counter ); println!("Avg Energy: {:.3}", final_energy);