Monte carlo simulation for hard spheres written in Rust language

montecarlo.rs

//
// Monte carlo simulation for hard spheres written in Rust language.
//

//To compile and execute this (on Mac OSX):
// rustc -O test.rs -L/opt/local/lib
// ./test
//

extern mod core;

extern mod gsl {
  fn gsl_rng_uniform(rng: *gsl_rng) -> c_double;
  const gsl_rng_default: *gsl_rng_type;
  fn gsl_rng_env_setup() -> *gsl_rng_type;
  fn gsl_rng_alloc(t: *gsl_rng_type) -> *gsl_rng;
}

use libc::*;

struct gsl_rng_type {
    name: @str,
    max: uint,
    min: uint,
    size: size_t,
    set: @fn(state: @c_void, seed: uint),
    get_double: @fn(state: @c_void) -> c_double
  }

struct gsl_rng{
  _type: @gsl_rng_type,
  state: @c_void
}


//extern mod c_float_utils;
use core::uint::*;
use core::rand::*;
use core::io::*;
use core::vec::*;
use core::str::*;
//use c_float_utils::*;

struct Coord {
  x: float,
  y: float,
  z: float
}

impl Coord {
  fn norm_sq(&self) -> float {
    self.x * self.x + self.y*self.y + self.z * self.z
  }
}

//If I use ~Coord, this slows down significantly for some reason.
// #[inline(always)] // <- This does not have effect on speed with -O compile option.
fn diff_coord(a: &Coord, b: &Coord) -> Coord {
  Coord{x:a.x-b.x,y:a.y-b.y,z:a.z-b.z}
}


fn init_particles(L: float, Nside: uint) -> ~[Coord] {
  let mut count: uint =0;
  let mut c: ~[Coord] = vec::with_capacity(Nside*Nside*Nside+2);
  c.push(Coord{x:0f,y:0f,z:0f});
  for range(0,Nside) |i| {
    for range(0,Nside) |j| {
      for range(0,Nside) |k| {
        c.push(Coord{
          x:((i+1) as float)*L/(Nside as float),
          y:((j+1) as float)*L/(Nside as float),
          z:((k+1) as float)*L/(Nside as float)});
      }
    }
  }
  c
}

fn round(v: float) -> float { v.to_int() as float }

fn distance_sq(c1: &Coord, c2: &Coord, L: float) -> float{
  let mut delta = diff_coord(c1,c2);
  
  //Use of the following closure slows down signigicantly.
  //let normalize = |v: float| -> float {v - L*round(v/L)};
  //delta.x = normalize(delta.x);

  delta.x -= L*round(delta.x/L);
  delta.y -= L*round(delta.y/L);
  delta.z -= L*round(delta.z/L);

  delta.norm_sq()
}

fn judge_overlap(c: &[Coord], c2: &Coord, L: float, idx: uint) -> bool {
  distance_sq(&c[idx],c2,L) < 1f
}

fn print_trajectory(traj: &[Coord]) {
  let mut i = 0;
  for vec::each(traj) |t| {
    i += 1;
    println(fmt!("%d\t%f\t%f\t%f",i+1 as int,t.x,t.y,t.z));
  }
}

fn vec_repeats<T: Copy>(n: uint, obj: T) -> ~[T] {
  let mut v: ~[T] = vec::with_capacity(n);
  for range(0,n) |_| {
    v.push(copy obj);
  }
  v
}

fn mk_range(s: uint,e: uint) -> ~[uint] {
  let mut v: ~[uint] = vec::with_capacity(e-s);
  for range(s,e) |i| {
    v.push(i);
  }
  v
}

fn mk_ranr() -> *gsl_rng {
  let mut ranr: *gsl_rng;
 unsafe{
  let ranT: *gsl_rng_type = gsl::gsl_rng_default;
  gsl::gsl_rng_env_setup();
  ranr = gsl::gsl_rng_alloc(ranT);
  gsl::gsl_rng_uniform(ranr);
  }
  ranr
}

fn main() {
  let rng = Rng();

  const Nside: uint = 10;
  const N: uint = Nside * Nside * Nside;

  const rho: float = 0.5;
  const L: float = 12.6;
  //const L: float = std::cmath::pow(1.0/rho*N,0.33333);
  print(fmt!("Nside = %d, d = %f, N = %d, L = %f\n",Nside as int,1.0,N as int,L));
  let mut c: ~[Coord] = init_particles(L,Nside);

  const disp_size: float=1f;
  const nsweeps:  uint = 1000;

  let mut traj: ~[Coord] = vec::with_capacity(1000);
  let ns = mk_range(1,N+1);
  
  let ranr: *gsl_rng = mk_ranr();
  
  for range(1,nsweeps+1) |sweep| {
    if sweep % 50 == 0 {
      println(sweep.to_str());
    }
    for range(1,N+1) |step| {
      //let imove: uint = (rng.gen_f64()*(N as f64)).to_int() as uint;
      let imove: uint = f64::ceil(gsl::gsl_rng_uniform(ranr)*(N as f64)) as uint;
      
      //let get_disp = |d: float| -> float {d*(rng.gen_f64() as float-0.5)};
      let ctrial = unsafe{
      let get_disp = |d: float| -> float {d*(gsl::gsl_rng_uniform(ranr) as float-0.5)};
      
      Coord{
            x: c[imove].x + get_disp(disp_size),
            y: c[imove].y + get_disp(disp_size),
            z: c[imove].z + get_disp(disp_size)}
            };

      let mut overlap = false;
      //Check overlap for all pairs (imove,1),(imove,2),...,(imove,N) except (imove,imove)
      for vec::each(ns) |&i| {
        if imove == i || judge_overlap(c,&ctrial,L,i) {
          overlap = true;
          break;
        }
      }
      if overlap {
        c[imove] = ctrial;
      }
//This is slower somehow.
/*      if ns.all(|&i| {
        imove == i || !judge_overlap(c,&ctrial,L,i)
      }) {
        c[imove] = ctrial;
      }*/
    }
    traj.push(c[445]);
  }

  print_trajectory(traj);
}