use core::panic;
use nalgebra::{Point3, Rotation3, Vector3};
use qc_file_parsers::{Xyz, XyzLine};
use rotate::rotator;
use std::env;
use std::error::Error;
use std::fs::File;
use std::io::{BufReader, BufWriter, Write};
fn print_usage() {
eprintln!("Usage:\nrotator [xyz-file] [index 1] [index 2] -p [[index 3]] -a [[y,z]] -u [[unit]]\n
molecular coordinates are rotated. \n
Rotation matrix is constructed from distance vector of atom [index 1] to atom [index 2] and target axis option -a [[y z]] (default x).
It is possible to define the units of the vector components in 'bohr' or 'ang' (option -u [[unit]] default 'ang' === angstroem)
Furthermore -p defines a plane using a third plane and rotates its normal vector to be || to -a.
");
}
fn get_molecule(parsed: Xyz) -> Vec<(String, nalgebra::OPoint<f32, nalgebra::Const<3>>)> {
let iter_lines = parsed.lines.iter();
let mut molecule: Vec<(_, Point3<f32>)> = Vec::new();
for l in iter_lines {
match l {
XyzLine::Numeric(n) => molecule.push((n.z_value.to_string(), n.xyz)),
XyzLine::Symbolic(s) => molecule.push((s.symbol.clone(), s.xyz)),
}
}
molecule
}
fn get_user_axis(a: String) -> char {
if a.eq("x") {
'i'
} else if a.eq("y") {
'j'
} else if a.eq("z") {
'k'
} else {
panic!("Unknown axis {}!", a)
}
}
fn get_user_unit(u: String) -> &'static str {
if u.eq("ang") {
"ang"
} else if u.eq("bohr") {
"bohr"
} else {
panic!("Unkown unit {}!", u)
}
}
fn main() -> Result<(), Box<dyn Error>> {
let mut uinp: env::Args = env::args();
if uinp.len() < 4 {
print_usage();
}
let _ = uinp.next();
let in_fname = uinp.next().unwrap();
let f = File::open(in_fname.clone())?;
let mut index1 = std::usize::MAX;
let mut index2 = std::usize::MAX;
match uinp.next().unwrap().parse::<usize>() {
Ok(i) => index1 = i - 1,
Err(i) => {
println!("Could not parse {} as unsigned integer.", i);
panic!("Cannot parse atomic index!");
}
}
match uinp.next().unwrap().parse::<usize>() {
Ok(i) => index2 = i - 1,
Err(i) => {
println!("Could not parse {} as unsigned integer.", i);
panic!("Cannot parse atomic index!");
}
}
let mut axis = 'i';
let mut unit = "ang";
let mut index3 = std::usize::MAX;
if let Some(s) = uinp.next() {
let options = uinp.next().unwrap();
if s.eq("-a") {
axis = get_user_axis(options);
} else if s.eq("-u") {
unit = get_user_unit(options);
} else if s.eq("-p") {
index3 = options.parse::<usize>().unwrap() - 1;
} else {
panic!("Unknown Option {}", s);
}
}
if let Some(s) = uinp.next() {
let options = uinp.next().unwrap();
if s.eq("-a") {
axis = get_user_axis(options);
} else if s.eq("-u") {
unit = get_user_unit(options);
} else if s.eq("-p") {
index3 = options.parse::<usize>().unwrap() - 1;
} else {
panic!("Unknown Option {}", s);
}
}
if let Some(s) = uinp.next() {
let options = uinp.next().unwrap();
if s.eq("-a") {
axis = get_user_axis(options);
} else if s.eq("-u") {
unit = get_user_unit(options);
} else if s.eq("-p") {
index3 = options.parse::<usize>().unwrap() - 1;
} else {
panic!("Unknown Option {}", s);
}
}
let mut reader = BufReader::new(f);
let parsed = Xyz::new(&mut reader, unit)?;
println!(
"Parsed xyz file {} with {} atoms",
in_fname, parsed.number_of_atoms
);
let mut molecule = get_molecule(parsed);
let point_a1 = molecule.get(index1).unwrap().1;
let point_a2 = molecule.get(index2).unwrap().1;
let distance_vector: Vector3<f32> = point_a1 - point_a2;
let mut rotate_to = distance_vector;
if index3.lt(&std::usize::MAX) {
let point_a3 = molecule.get(index3).unwrap().1;
let scnd_d_vector: Vector3<f32> = point_a1 - point_a3;
rotate_to = distance_vector.cross(&scnd_d_vector);
}
let target_vector = rotator::get_unit(axis);
let rot_mat = Rotation3::rotation_between(&rotate_to, &target_vector).unwrap();
println!(
"Rotation angle is: {} to align {:?} and {:?}",
rot_mat.angle() * 180.0 / std::f32::consts::PI,
distance_vector,
target_vector
);
let out_file = File::create("rotated.xyz")?;
let mut out_buffer = BufWriter::new(out_file);
out_buffer
.write_all(format!("{}\n\n", molecule.len()).as_bytes())
.unwrap();
for tuple in &mut molecule {
print!("{:?} -> ", tuple.1);
tuple.1 = rot_mat * tuple.1;
println!("{:?}", tuple.1);
out_buffer
.write_all(
format!(
"{} {:.5} {:.5} {:.5}\n",
tuple.0.to_uppercase(), tuple.1.x, tuple.1.y, tuple.1.z
)
.as_bytes(),
)
.unwrap();
}
out_buffer.flush().unwrap();
println!("Wrote rotated structure to rotated.xyz");
Ok(())
}
#[cfg(test)]
mod unit_tests {
use crate::{get_user_axis, get_user_unit};
#[test]
#[should_panic]
fn test_get_unknown_user_axis() {
get_user_axis(String::from("vorn"));
get_user_axis(String::from("f"));
get_user_axis(String::from("k"));
}
#[test]
fn test_get_user_axis() {
assert_eq!('i', get_user_axis("x".to_string()));
assert_eq!('j', get_user_axis("y".to_string()));
assert_eq!('k', get_user_axis("z".to_string()));
}
#[test]
#[should_panic]
fn test_get_unknown_user_unit() {
get_user_unit(String::from("frobel"));
get_user_unit(String::from("a"));
get_user_unit(String::from("an"));
get_user_unit(String::from("angh"));
get_user_unit(String::from("bor"));
}
#[test]
fn test_get_user_unit() {
assert_eq!("ang", get_user_unit(String::from("ang")));
assert_eq!("bohr", get_user_unit(String::from("bohr")));
}
}