From f95c1da06aacb31bb0adaf5730b69a06049351a5 Mon Sep 17 00:00:00 2001 From: danijoo Date: Sat, 21 Jan 2017 14:03:08 +0100 Subject: [PATCH] initial --- .gitignore | 1 + .idea/compiler.xml | 22 ++ .idea/copyright/profiles_settings.xml | 3 + .idea/libraries/Cargo__mclj_.xml | 10 + .idea/libraries/Rust__mclj_.xml | 13 + .idea/misc.xml | 26 ++ .idea/modules.xml | 8 + .idea/workspace.xml | 525 ++++++++++++++++++++++++++ Cargo.lock | 23 ++ Cargo.toml | 14 + mclj.iml | 16 + src/energy.rs | 88 +++++ src/main.rs | 167 ++++++++ 13 files changed, 916 insertions(+) create mode 100644 .gitignore create mode 100644 .idea/compiler.xml create mode 100644 .idea/copyright/profiles_settings.xml create mode 100644 .idea/libraries/Cargo__mclj_.xml create mode 100644 .idea/libraries/Rust__mclj_.xml create mode 100644 .idea/misc.xml create mode 100644 .idea/modules.xml create mode 100644 .idea/workspace.xml create mode 100644 Cargo.lock create mode 100644 Cargo.toml create mode 100644 mclj.iml create mode 100644 src/energy.rs create mode 100644 src/main.rs diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..eb5a316 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +target diff --git a/.idea/compiler.xml b/.idea/compiler.xml new file mode 100644 index 0000000..96cc43e --- /dev/null +++ b/.idea/compiler.xml @@ -0,0 +1,22 @@ + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/.idea/copyright/profiles_settings.xml b/.idea/copyright/profiles_settings.xml new file mode 100644 index 0000000..e7bedf3 --- /dev/null +++ b/.idea/copyright/profiles_settings.xml @@ -0,0 +1,3 @@ + + + \ No newline at end of file diff --git a/.idea/libraries/Cargo__mclj_.xml b/.idea/libraries/Cargo__mclj_.xml new file mode 100644 index 0000000..e2d24a6 --- /dev/null +++ b/.idea/libraries/Cargo__mclj_.xml @@ -0,0 +1,10 @@ + + + + + + + + + + \ No newline at end of file diff --git a/.idea/libraries/Rust__mclj_.xml b/.idea/libraries/Rust__mclj_.xml new file mode 100644 index 0000000..14758e9 --- /dev/null +++ b/.idea/libraries/Rust__mclj_.xml @@ -0,0 +1,13 @@ + + + + + + + + + + + + + \ No newline at end of file diff --git a/.idea/misc.xml b/.idea/misc.xml new file mode 100644 index 0000000..6959d37 --- /dev/null +++ b/.idea/misc.xml @@ -0,0 +1,26 @@ + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/.idea/modules.xml b/.idea/modules.xml new file mode 100644 index 0000000..e49e34f --- /dev/null +++ b/.idea/modules.xml @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file diff --git a/.idea/workspace.xml b/.idea/workspace.xml new file mode 100644 index 0000000..917aadb --- /dev/null +++ b/.idea/workspace.xml @@ -0,0 +1,525 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + true + DEFINITION_ORDER + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1484992090438 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + No facets are configured + + + + + + + + + + + + + + + 1.8 + + + + + + + + mclj + + + + + + + + Cargo <mclj> + + + + + + + + \ No newline at end of file diff --git a/Cargo.lock b/Cargo.lock new file mode 100644 index 0000000..e8d5b5c --- /dev/null +++ b/Cargo.lock @@ -0,0 +1,23 @@ +[root] +name = "mclj" +version = "0.1.0" +dependencies = [ + "rand 0.3.15 (registry+https://github.com/rust-lang/crates.io-index)", +] + +[[package]] +name = "libc" +version = "0.2.20" +source = "registry+https://github.com/rust-lang/crates.io-index" + +[[package]] +name = "rand" +version = "0.3.15" +source = "registry+https://github.com/rust-lang/crates.io-index" +dependencies = [ + "libc 0.2.20 (registry+https://github.com/rust-lang/crates.io-index)", +] + +[metadata] +"checksum libc 0.2.20 (registry+https://github.com/rust-lang/crates.io-index)" = "684f330624d8c3784fb9558ca46c4ce488073a8d22450415c5eb4f4cfb0d11b5" +"checksum rand 0.3.15 (registry+https://github.com/rust-lang/crates.io-index)" = "022e0636ec2519ddae48154b028864bdce4eaf7d35226ab8e65c611be97b189d" diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..6f478bf --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,14 @@ +[package] +name = "mclj" +version = "0.1.0" +authors = ["Daniel Bauer "] + +[dependencies] +rand = "0.3.15" + +[profile.release] +lto = true +opt-level = 3 + + + diff --git a/mclj.iml b/mclj.iml new file mode 100644 index 0000000..6e67e8a --- /dev/null +++ b/mclj.iml @@ -0,0 +1,16 @@ + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/src/energy.rs b/src/energy.rs new file mode 100644 index 0000000..bd43d29 --- /dev/null +++ b/src/energy.rs @@ -0,0 +1,88 @@ + +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) -> (f64, f64) { + let mut energy = 0.0; + let mut virial = 0.0; + 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 { + 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); + energy += e; + virial += v; + } + } + } + energy += num_particles as f64 * e_corr; + return (energy, virial); +} + +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) -> (f64, f64) { + let mut energy = 0.0; + let mut virial = 0.0; + 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 { + 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); + if dist_squared < cutoff_squared { + let (e,v) = eval_pair_energy(dist_squared); + energy += e; + virial += v; + } + } + 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 { + let mut dx = x1 - x2; + let mut dy = y1 - y2; + 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} + + return dx*dx + dy*dy + dz*dz; +} + +#[test] +fn test_get_particle_distance_squared() { + let (x1, y1, z1) = (0.0, 0.0, 0.0); + 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); + + // 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); +} + +fn eval_pair_energy(dist_squared: 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); + let virial = 48.0 * ::LJ_EPS * ( r62 - 0.5 * r6 ); + return (energy, virial); +} + +#[test] +fn test_eval_pair_energy() { + let (e,v) = eval_pair_energy(1.0); + assert!( (e - 0.0).abs() < 0.00001, "{}", e); + + let (e,v) = eval_pair_energy(2.0); + assert!( (e - -0.4375).abs() < 0.00001, "{}", e); + + let (e,v) = eval_pair_energy(0.5); + assert!( (e - 224.0).abs() < 0.00001, "{}", e); +} \ No newline at end of file diff --git a/src/main.rs b/src/main.rs new file mode 100644 index 0000000..30c03c3 --- /dev/null +++ b/src/main.rs @@ -0,0 +1,167 @@ +#![allow(non_snake_case)] + +extern crate rand; +use rand::Rng; +use rand::distributions::{IndependentSample, Range}; +mod energy; +use energy::*; +use std::io::prelude::*; + +const LJ_EPS : f64 = 1.0; +const LJ_SIG : f64 = 1.0; + +const TAILCORR : bool = true; +const SHIFT: bool = false; + +// easy printing to stderr +macro_rules! println_stderr( + ($($arg:tt)*) => { { + let r = writeln!(&mut ::std::io::stderr(), $($arg)*); + r.expect("failed printing to stderr"); + } } +); + +fn main() { + + // define all the stuff + let sample_steps = 1000000; + let minim_steps = 1000000; + + let num_particles: usize = 512; + let density = 0.7; + let temperature = 0.9; + + let cutoff = 3.0; + let displacement = 0.1; + + + println_stderr!(""); + println_stderr!("################################################################"); + println_stderr!("################## LJ Monte Carlo Simulation #################"); + println_stderr!("################################################################"); + println_stderr!(""); + + + // 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 cutoff_squared = cutoff * cutoff; + + let mut rng = rand::thread_rng(); + let particle_range = Range::new(0, num_particles-1); + + 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; } + } + + 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}, Pressurecprr: {:8.3}", e_corr, 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); + 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!(""); + + for step in 0..minim_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); + + // 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); + 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); + 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; + + + if step_counter % 5000 == 0 && step < minim_steps { + println!("Minim {}\tEnergy: {:.3}\tVirial: {:.3}\tAcceptance:{:.1}\tDisplacement: {:.3}", step_counter, energy, virial, 666, displacement); + } + + // reset sums for sampling + if step == minim_steps-1 { + println!("Starting averaging!"); + step_counter = 0; + energy_sum = 0.0; + virial_sum = 0.0; + } + + } + + + 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 pressure = virial_sum / 3.0 / step_counter as f64 / volume + density * temperature + p_corr; + println!("Steps: {}", step_counter ); + println!("Avg Energy: {:.3}", final_energy); + println!("Energy/Particle: {:.3}", particle_energy); + println!("Virial: {:.3}", final_virial); + println!("Pressure: {:.3}", pressure); + + +} + + +