bauer/lj_monte_carlo
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

surface tension fix

Browse Source
This commit is contained in:
danijoo
2017-01-23 21:46:45 +01:00
parent 72ca57d357
commit 29768e5d1e
1 changed files with 37 additions and 34 deletions
Download Patch File Download Diff File Expand all files Collapse all files

71
src/surface_tension.rs
View File

@@ -7,11 +7,10 @@ const LJ_EPS : f64 = 1.0;
const LJ_SIG : f64 = 1.0; const LJ_SIG : f64 = 1.0;
fn get_derived_pair_potential(distance: f64) -> f64 { fn get_virial(distance_sqr: f64) -> f64 {
let r = LJ_SIG/distance; let r2 = LJ_SIG.powi(2)/distance_sqr;
let r3 = r*r*r; let r6 = r2 * r2 * r2;
let r4 = r*r*r*r; return 48.0 * LJ_EPS / LJ_SIG * ( r6*r6 - 0.5*r6 )
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 { fn get_distance_with_pbc(x1: f64, x2: f64, length: f64, half_length: f64) -> f64 {
@@ -23,47 +22,50 @@ fn get_distance_with_pbc(x1: f64, x2: f64, length: f64, half_length: f64) -> f64
fn main() { fn main() {
let mut trj_reader = TrjReader::new(&"2phases/10k_1step.xyz".to_string()); let mut trj_reader = TrjReader::new(&"2phases/1kk_100step.xyz".to_string());
let mut frame = trj_reader.next_frame(); let mut frame = trj_reader.next_frame();
let volume = frame.box_x * frame.box_y * frame.box_z; let volume = frame.box_x * frame.box_y * frame.box_z;
let density = frame.num_particles as f64 / volume; let density = frame.num_particles as f64 / volume;
let num_particles = frame.num_particles;
let box_half_x = frame.box_x / 2.0; let box_half_x = frame.box_x / 2.0;
let box_half_y = frame.box_y / 2.0; let box_half_y = frame.box_y / 2.0;
let box_half_z = frame.box_z / 2.0; let box_half_z = frame.box_z / 2.0;
let mut frame_count = 0; let mut frame_count = 0;
let mut trace_xy_sum = 0.0; let mut p_xy_sum = 0.0;
let mut trace_z_sum = 0.0; // let mut p_y_sum = 0.0;
let mut p_normal_sum = 0.0; let mut p_z_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; let variable_without_name = frame.temperature/LJ_EPS * density;
loop { loop {
frame_count += 1; frame_count += 1;
let mut trace_xy = 0.0; let mut trace_xy = 0.0;
// let mut trace_y = 0.0;
let mut trace_z = 0.0; let mut trace_z = 0.0;
for i in 0..frame.num_particles { for i in 0..num_particles {
for j in i+1..frame.num_particles { for j in i+1..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 dist_sqrt = get_particle_distance_squared(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 dist = dist_sqrt.sqrt();
let dx = get_distance_with_pbc(frame.rx[i], frame.rx[j], frame.box_x, box_half_x); 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 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); 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); let virial = get_virial(dist_sqrt);
trace_z += (dz * dz) / dist * get_derived_pair_potential(dist); trace_xy += (dx * dx + dy * dy) / dist * virial;
// trace_y += (dy * dy) / dist * virial;
trace_z += (dz * dz) / dist * virial;
} }
} }
let p_tangial = variable_without_name - 1.0/(2.0*volume)*trace_xy; let p_xy = variable_without_name - 1.0/(2.0*volume)*(trace_xy/num_particles as f64);
let p_normal = variable_without_name - 1.0/volume*trace_z; // let p_yy = variable_without_name - 1.0/volume*(trace_y/num_particles);
let p_diff = p_normal - p_tangial_sum; let p_zz = variable_without_name - 1.0/volume*(trace_z/num_particles as f64);
p_diff_sum += p_diff;
p_tangial_sum += p_tangial; p_xy_sum += p_xy;
p_normal_sum += p_normal; // p_y_sum += p_yy;
trace_xy_sum += trace_xy; p_z_sum += p_zz;
trace_z_sum += trace_z;
/////////////////////////////////// ///////////////////////////////////
@@ -72,18 +74,19 @@ fn main() {
} }
if frame_count % 100 == 0 { if frame_count % 100 == 0 {
// let trace_xy_avg = trace_xy_sum / frame_count as f64; let p_z_avg = p_z_sum / frame_count as f64;
// let trace_z_avg = trace_z_sum / frame_count as f64; let p_xy_avg = p_xy_sum / frame_count as f64;
let p_tangial = p_tangial_sum / frame_count as f64; let p_diff = p_z_avg - p_xy_avg;
let p_normal = p_normal_sum / frame_count as f64; let surface_tension = (frame.box_z/2.0)*p_diff;
let surface_tension = frame.box_z/2.0*(p_diff_sum/frame_count as f64); println!("{} zz: {} xy: {} diff: {} tension: {}", frame_count, p_z_avg, p_xy_avg, p_diff, surface_tension);
println!("{} tangial: {} normal: {} diff: {} tension: {}", frame_count, p_tangial, p_normal, p_diff_sum/frame_count as f64, surface_tension);
frame_count = 0; frame_count = 0;
trace_xy_sum = 0.0; p_z_sum = 0.0;
trace_z_sum = 0.0; p_xy_sum = 0.0;
p_tangial_sum = 0.0; // trace_xy_sum = 0.0;
p_normal_sum = 0.0 // trace_z_sum = 0.0;
// p_tangial_sum = 0.0;
// p_normal_sum = 0.0
} }