surface tension and widom particle insertion (WIP)

src/energy.rs

@@ -5,8 +5,8 @@ pub fn get_total_energy(rx: &[f64], ry: &[f64], rz: &[f64], num_particles: usize
     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 {
+    for i in 0..num_particles {
+        for j in i+1..num_particles {
             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, e_shift);
@@ -25,7 +25,7 @@ pub fn get_particle_energy(rx: &[f64], ry: &[f64], rz: &[f64], p_index: usize, n
     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 i in 0..num_particles {
         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);
@@ -38,7 +38,7 @@ pub fn get_particle_energy(rx: &[f64], ry: &[f64], rz: &[f64], p_index: usize, n
     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 {
+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;
     let mut dy = y1 - y2;
     let mut dz = z1 - z2;
@@ -53,6 +53,10 @@ fn get_particle_distance_squared(x1: f64,y1: f64,z1: f64,x2: f64,y2: f64,z2: f64
     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();
+}
+
 #[test]
 fn test_get_particle_distance_squared() {
     let (x1, y1, z1) = (0.0, 0.0, 0.0);
@@ -65,7 +69,7 @@ fn test_get_particle_distance_squared() {
     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, e_shift: f64) -> (f64, f64) {
+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;
     let energy = 4.0 * ::LJ_EPS * (r62 - r6) - e_shift;

src/main.rs

@@ -124,7 +124,7 @@ fn main() {
     let max_displacement = length / 2.0;
 
     let mut rng = rand::thread_rng();
-    let particle_range = Range::new(0, num_particles-1);
+    let particle_range = Range::new(0, num_particles);
 
     let mut rx : Vec<f64> = vec![];
     let mut ry : Vec<f64> = vec![];

src/surface_tension.rs

@@ -0,0 +1,115 @@
+mod trajectory;
+use trajectory::*;
+mod energy;
+use energy::*;
+
+const LJ_EPS : f64 = 1.0;
+const LJ_SIG : f64 = 1.0;
+
+
+fn get_derived_pair_potential(distance: f64) -> f64 {
+    let r = LJ_SIG/distance;
+    let r3 = r*r*r;
+    let r4 = r*r*r*r;
+    return 24.0 * LJ_EPS / LJ_SIG * ( r4*r3 - 2.0*(r4*r4*r4*r) )
+}
+
+fn get_distance_with_pbc(x1: f64, x2: f64, length: f64, half_length: f64) -> f64 {
+    let mut d = x1-x2;
+    if d > half_length { d -= length }
+    else if d < -half_length { d += length }
+    return d;
+}
+
+fn main() {
+
+    let mut trj_reader = TrjReader::new(&"2phases/10k_1step.xyz".to_string());
+    let mut frame = trj_reader.next_frame();
+
+    let volume = frame.box_x * frame.box_y * frame.box_z;
+    let density = frame.num_particles as f64 / volume;
+
+    let box_half_x = frame.box_x / 2.0;
+    let box_half_y = frame.box_y / 2.0;
+    let box_half_z = frame.box_z / 2.0;
+
+    let mut frame_count = 0;
+    let mut trace_xy_sum = 0.0;
+    let mut trace_z_sum = 0.0;
+    let mut p_normal_sum = 0.0;
+    let mut p_tangial_sum = 0.0;
+    let mut p_diff_sum = 0.0;
+
+    let variable_without_name = frame.temperature/LJ_EPS * density;
+    loop {
+        frame_count += 1;
+
+        let mut trace_xy = 0.0;
+        let mut trace_z = 0.0;
+        for i in 0..frame.num_particles {
+            for j in i+1..frame.num_particles {
+                let dist = get_particle_distance(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 dx = get_distance_with_pbc(frame.rx[i], frame.rx[j], frame.box_x, box_half_x);
+                let dy = get_distance_with_pbc(frame.ry[i], frame.ry[j], frame.box_y, box_half_y);
+                let dz = get_distance_with_pbc(frame.rz[i], frame.rz[j], frame.box_z, box_half_z);
+                trace_xy += (dx * dx  + dy * dy) / dist * get_derived_pair_potential(dist);
+                trace_z += (dz * dz) / dist * get_derived_pair_potential(dist);
+            }
+        }
+        let p_tangial = variable_without_name - 1.0/(2.0*volume)*trace_xy;
+        let p_normal = variable_without_name - 1.0/volume*trace_z;
+        let p_diff = p_normal - p_tangial_sum;
+        p_diff_sum += p_diff;
+        p_tangial_sum += p_tangial;
+        p_normal_sum += p_normal;
+        trace_xy_sum += trace_xy;
+        trace_z_sum += trace_z;
+
+        ///////////////////////////////////
+
+        if frame_count % 100 == 0 {
+            print!(".");
+        }
+
+        if frame_count % 100 == 0 {
+//            let trace_xy_avg = trace_xy_sum / frame_count as f64;
+//            let trace_z_avg = trace_z_sum / frame_count as f64;
+            let p_tangial = p_tangial_sum / frame_count as f64;
+            let p_normal = p_normal_sum / frame_count as f64;
+            let surface_tension = frame.box_z/2.0*(p_diff_sum/frame_count as f64);
+            println!("{}  tangial: {}   normal: {}    diff: {}    tension:   {}", frame_count, p_tangial, p_normal, p_diff_sum/frame_count as f64, surface_tension);
+
+            frame_count = 0;
+            trace_xy_sum = 0.0;
+            trace_z_sum = 0.0;
+            p_tangial_sum = 0.0;
+            p_normal_sum = 0.0
+
+        }
+
+        if !trj_reader.update_with_next(&mut frame) { break }
+    }
+
+}
+
+//        for slab in 0..NUM_SLABS {
+//
+//            // calculate slab density
+//            let slab_end_z = slab_height * slab as f64;
+//            let slab_start_z = slab_end_z - slab_height;
+//            let mut slab_particle_count = 0;
+//            for i in 0..frame.num_particles {
+//                if slab_start_z > frame.rz[i] && frame.rz[i] < slab_end_z {
+//                    slab_particle_count += 1;
+//                }
+//            }
+//            let slab_density = slab_particle_count as f64 / slab_volume;
+//
+//            for i in 0..frame.num_particles {
+//                for j in i+1..frame.num_particles {
+//
+//                }
+//            }
+//
+//            println!("Slab {}\tDensity: {}", slab, slab_density);
+//        }

src/trajectory.rs

@@ -64,8 +64,7 @@ pub struct TrjReader {
     pub reader: BufReader<File>,
 }
 impl TrjReader {
-    pub fn new() -> TrjReader {
-        let filename = "trajectory.xyz";
+    pub fn new(filename: &String) -> TrjReader {
         let file = File::open(filename).expect("Failed to open file.");
         let reader : BufReader<File> = BufReader::new(file);
         return TrjReader { reader:reader };
@@ -95,7 +94,7 @@ impl TrjReader {
         let mut rx = Vec::new();
         let mut ry = Vec::new();
         let mut rz = Vec::new();
-        for i in 0..num_particles-1{
+        for i in 0..num_particles{
             let atom_line = &mut String::new();
             match self.reader.read_line(atom_line) {
                 Ok(size) => {
@@ -132,7 +131,7 @@ impl TrjReader {
                     return false;
                 }
             },
-            Err(e) => panic!("no new frame!")
+            Err(e) => { return false; }
         }
         let first_line_vec : Vec<&str> = buffer_string.split_whitespace().collect();
         frame.num_particles = first_line_vec[0].parse::<usize>().unwrap();
@@ -144,7 +143,7 @@ impl TrjReader {
         frame.lj_eps = lj[0].parse::<f64>().unwrap();
         frame.lj_sig = lj[1].parse::<f64>().unwrap();
         frame.lj_cutoff = lj[2].parse::<f64>().unwrap();
-        for i in 0..frame.num_particles-1 {
+        for i in 0..frame.num_particles {
             let atom_line = &mut String::new();
             match self.reader.read_line(atom_line) {
                 Ok(size) => {

src/widom.rs

@@ -0,0 +1,138 @@
+mod trajectory;
+use trajectory::*;
+mod energy;
+use energy::*;
+extern crate rand;
+use rand::Rng;
+
+const LJ_EPS : f64 = 1.0;
+const LJ_SIG : f64 = 1.0;
+
+static MKSA_PLANCKS_CONSTANT_H : f64 = 1.0;
+static MASS : f64 = 1.0;
+
+const INSERTIONS_PER_STEP : usize = 10;
+const RUN_AVG_SIZE : usize = 100;
+
+const SHIFT : bool = true;
+
+// returns number of particles in liquid and gas phase as tuble
+fn count_particles(rz: &[f64], num_particles: usize, liquid_height: f64) -> (usize, usize) {
+    let mut liquid_count = 0;
+    let mut gas_count = 0;
+    for i in 0..num_particles {
+        if rz[i] < liquid_height { liquid_count += 1; }
+        else { gas_count += 1; }
+    }
+    return (liquid_count, gas_count);
+}
+
+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;
+    let half_l_y = l_y/2.0;
+    let half_l_z = l_z/2.0;
+    for i in 0..num_particles {
+
+        let dist_squared = get_particle_distance_squared(x, y, z, rx[i],ry[i],rz[i], l_x, l_y, l_z, half_l_x, half_l_y, half_l_z);
+        if dist_squared < cutoff_sqr {
+            let (e,v) = eval_pair_energy(dist_squared, e_shift);
+            energy += e;
+        }
+    }
+    return energy;
+}
+
+fn main() {
+
+    let mut trj_reader = TrjReader::new(&"2phases/10k_1step.xyz".to_string());
+    let mut frame = trj_reader.next_frame();
+
+    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;
+
+    let mut gas_slab = 2.0;
+    let mut liquid_height = frame.box_z / (gas_slab + 1.0);
+    let liquid_volume = volume / (gas_slab + 1.0);
+    let gas_volume = volume - liquid_volume;
+
+    let wave = MKSA_PLANCKS_CONSTANT_H / (2.0 * std::f64::consts::PI * MASS * frame.temperature / frame.lj_eps).sqrt();
+
+
+    let mut rng = rand::thread_rng();
+
+    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 mut frame_count = 0;
+    let mut widom_sum_gas = 0.0;
+    let mut ideal_sum_gas = 0.0;
+    let mut widom_sum_liquid = 0.0;
+    let mut ideal_sum_liquid = 0.0;
+
+    let mut avg_count = 0;
+    loop  {
+        frame_count += 1;
+        avg_count += 1;
+
+
+        // test particle liquid energy
+        let lx = frame.box_x * rng.gen::<f64>();
+        let ly = frame.box_y * rng.gen::<f64>();
+        let lz = liquid_height * rng.gen::<f64>();
+        let widom_e_liquid = get_particle_insertion_energy(&frame.rx, &frame.ry, &frame.rz, frame.num_particles, lx, ly, lz, frame.box_x, frame.box_y, frame.box_z, cutoff_sqr, e_shift);
+        widom_sum_liquid += (-beta*widom_e_liquid).exp();
+
+        // test particle gas energy
+        let gx = frame.box_x * rng.gen::<f64>();
+        let gy = frame.box_y * rng.gen::<f64>();
+        let gz = ((frame.box_z - liquid_height) * rng.gen::<f64>()) + liquid_height;
+        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 number of particles in each phase
+        let mut liquid_count = 0.0;
+        let mut gas_count = 0.0;
+        for i in 0..frame.num_particles {
+            if frame.rz[i] < liquid_height { liquid_count += 1.0; }
+            else { gas_count += 1.0; }
+        }
+
+        // calculate ideal gas potential for both phases
+        let wave = MKSA_PLANCKS_CONSTANT_H / (2.0 * std::f64::consts::PI * MASS * frame.temperature / frame.lj_eps).sqrt();
+        ideal_sum_gas += frame.temperature / frame.lj_eps * (gas_volume/(wave.powi(3)* gas_count)).ln();
+        ideal_sum_liquid += frame.temperature / frame.lj_eps * (liquid_volume/(wave.powi(3)* liquid_count)).ln();
+
+
+
+        if avg_count == RUN_AVG_SIZE {
+            let ideal_gas = ideal_sum_gas / avg_count as f64;
+            let ideal_liquid = ideal_sum_liquid / 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;
+            println!("Frame {}   gas {}   liquid {}  particlecount:{}/{}", frame_count, ideal_gas + excess_gas, ideal_liquid + excess_liquid, gas_count, liquid_count);
+            avg_count = 0;
+            ideal_sum_gas = 0.0;
+            ideal_sum_liquid = 0.0;
+            widom_sum_gas = 0.0;
+            widom_sum_liquid = 0.0;
+
+        }
+
+        //            liquid_gas_sum = 0.0;
+        //            excess_liquid_sum = 0.0;
+        //            gas_gas_sum = 0.0;
+        //            excess_gas_sum = 0.0;
+        //            avg_count = 0;
+
+
+
+
+        if !trj_reader.update_with_next(&mut frame) { break }
+    }
+
+
+
+
+
+}