diff --git a/src/energy.rs b/src/energy.rs index c27394c..360f77c 100644 --- a/src/energy.rs +++ b/src/energy.rs @@ -99,4 +99,19 @@ fn test_eval_pair_energy() { let (e,v) = eval_pair_energy(0.5, 0.0); assert!( (e - 224.0).abs() < 0.00001, "{}", e); -} \ No newline at end of file +} + + +pub fn eval_virial(distance: f64, LJ_EPS: f64, LJ_SIG: f64) -> f64 { + let r7 = (LJ_SIG/distance).powi(7); + let r13 = (LJ_SIG/distance).powi(13); + return 24.0 * LJ_EPS / LJ_SIG * ( r7-2.0*r13 ); +} + +#[test] +fn test_eval_virial() { + let dist = 1.5; + let result = eval_virial(dist, 1.0, 1.0, 0.0); + let expected = 1.1580288; + assert!( (result - expected).abs() < 0.0001, "{}", result ); +} diff --git a/src/main.rs b/src/main.rs index 15b0397..2d22452 100644 --- a/src/main.rs +++ b/src/main.rs @@ -237,7 +237,8 @@ fn main() { let tries_per_step : f64 = step_counter as f64 /accept_counter as f64; let acceptance_rate = 1.0/tries_per_step * 100.0; let avg_energy = energy_sum / step_counter as f64; - println_stderr!("Minim {:<10} Energy: {:<30.3} Accept.: {:<4.1}% dr: {:.3}", step+1, avg_energy, acceptance_rate, displacement); + let avg_virial = virial_sum / step_counter as f64; + println_stderr!("Minim {:<10} Energy: {:<30.3} Virial: {:<30.3} Accept.: {:<4.1}% dr: {:.3}", step+1, avg_energy, avg_virial, acceptance_rate, displacement); if SCALE { let scale_factor = (TRIES_INTENDED/tries_per_step * DISP_SCALE_FACTOR).abs(); diff --git a/src/surface_tension.rs b/src/surface_tension.rs index 8d80a37..ccb7890 100644 --- a/src/surface_tension.rs +++ b/src/surface_tension.rs @@ -8,19 +8,6 @@ const LJ_EPS : f64 = 1.0; const LJ_SIG : f64 = 1.0; -fn get_virial(distance: f64) -> f64 { - let r7 = (LJ_SIG/distance).powi(7); - let r13 = (LJ_SIG/distance).powi(13); - return 24.0 * LJ_EPS / LJ_SIG * ( r7-2.0*r13 ); -} - -#[test] -fn test_get_viral() { - let result = get_virial(2.5); - let expected = 0.038999477; - assert!( (result - expected) < 0.0001, "{}", result ); -} - fn get_distance_with_pbc(x1: f64, x2: f64, length: f64, half_length: f64) -> f64 { let mut d = (x1-x2).abs(); if d > half_length { d -= length } @@ -89,7 +76,7 @@ fn main() { 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); - let virial = get_virial(dist); + let virial = eval_virial(dist, LJ_EPS, LJ_SIG); trace_xy += (dx * dx + dy * dy) / dist * virial; trace_z += (dz * dz) / dist * virial; }