initial

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);
+}

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<f64> = vec![];
+    let mut ry : Vec<f64> = vec![];
+    let mut rz : Vec<f64> = vec![];
+    loop {
+        rx.push(l_x * rng.gen::<f64>());
+        ry.push(l_y * rng.gen::<f64>());
+        rz.push(l_z * rng.gen::<f64>());
+        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::<f64>() - 0.5 ) * displacement;
+        ry[rnd_index] += ( rng.gen::<f64>() - 0.5 ) * displacement;
+        rz[rnd_index] += ( rng.gen::<f64>() - 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::<f64>() < (-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);
+
+
+}
+
+
+