use std::io::BufReader;
use std::path::Path;
use std::{io, str::FromStr};
use glam::Vec3;
use crate::structure::{Atom, AtomName, AtomNum, BoxVecs, ResName, ResNum, Structure};
#[non_exhaustive]
#[derive(Debug)]
pub enum ParseGroError {
IOError(io::Error),
UnexpectedEOF(ExpectedItem),
UnexpectedEOL(ExpectedItem),
BadNAtoms(std::num::ParseIntError),
InsufficientAtoms { exp: usize, enc: usize },
BadResnum(std::num::ParseIntError),
BadAtomnum(std::num::ParseIntError),
BadPositionValue(std::num::ParseFloatError),
BadVelocityValue(std::num::ParseFloatError),
BadBoxVecsLength(usize),
BadBoxVecsValue(std::num::ParseFloatError),
}
#[non_exhaustive]
#[derive(Debug, Clone, Copy)]
pub enum ExpectedItem {
Title, NAtoms, BoxVecs,
Resnum, Resname, Atomname, Atomnum, Position, Velocity, }
impl std::fmt::Display for ParseGroError {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
match self {
Self::IOError(err) => write!(f, "I/O error: {err:?}"),
Self::UnexpectedEOF(expected) => {
write!(f, "unexpected end of file: expected {expected}")
}
Self::UnexpectedEOL(expected) => {
write!(f, "unexpected end of line: while reading {expected}")
}
Self::BadNAtoms(err) => write!(f, "could not read the number of atoms: {err}"),
Self::InsufficientAtoms { exp, enc } => write!(
f,
"could not read the specified number of atoms: \
expected {exp} but could only read {enc}"
),
Self::BadResnum(err) => write!(f, "could not read the residue number: {err}"),
Self::BadAtomnum(err) => write!(f, "could not read the atom number: {err}"),
Self::BadPositionValue(err) => write!(f, "could not read position value: {err}"),
Self::BadVelocityValue(err) => write!(f, "could not read velocity value: {err}"),
Self::BadBoxVecsLength(n) => {
write!(
f,
"could not read box vectors: expected 3 or 9 values, found {n}"
)
}
Self::BadBoxVecsValue(err) => write!(f, "could not read box vectors value: {err}"),
}
}
}
impl std::error::Error for ParseGroError {}
impl std::fmt::Display for ExpectedItem {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
f.write_str(match self {
Self::Title => "title",
Self::NAtoms => "number of atoms",
Self::BoxVecs => "box vectors",
Self::Resnum => "resnum",
Self::Resname => "resname",
Self::Atomname => "atomname",
Self::Atomnum => "atomnum",
Self::Position => "position",
Self::Velocity => "velocity",
})
}
}
impl From<ParseGroError> for io::Error {
fn from(err: ParseGroError) -> Self {
match err {
ParseGroError::IOError(err) => err,
other => io::Error::other(other),
}
}
}
impl From<io::Error> for ParseGroError {
fn from(err: io::Error) -> Self {
ParseGroError::IOError(err)
}
}
pub type Result<T> = std::result::Result<T, ParseGroError>;
pub mod ranges {
use std::ops::Range;
pub const RESNUM: Range<usize> = 0..5;
pub const RESNAME: Range<usize> = 5..10;
pub const ATOMNAME: Range<usize> = 10..15;
pub const ATOMNUM: Range<usize> = 15..20;
pub const POSITION: Range<usize> = 20..44;
pub const VELOCITY: Range<usize> = 44..68;
pub const MIN_LINE_LEN: usize = POSITION.end;
}
pub struct ParseList {
pub resnum: bool,
pub resname: bool,
pub atomname: bool,
pub atomnum: bool,
pub position: bool,
pub velocity: bool,
}
impl ParseList {
pub const ALL: Self = Self {
resnum: true,
resname: true,
atomname: true,
atomnum: true,
position: true,
velocity: true,
};
}
pub trait ReadGro<A>: Sized {
const PARSE_LIST: ParseList = ParseList::ALL;
fn build_atom(
resnum: Option<ResNum>,
resname: Option<ResName>,
atomname: Option<AtomName>,
atomnum: Option<AtomNum>,
position: Option<[f32; 3]>,
velocity: Option<[f32; 3]>,
) -> A;
fn build_structure(title: String, atoms: Vec<A>, boxvecs: BoxVecs) -> Self;
#[inline]
fn parse_resnum(line: &str) -> Result<ResNum> {
line.get(ranges::RESNUM)
.ok_or(ParseGroError::UnexpectedEOL(ExpectedItem::Resnum))?
.trim()
.parse()
.map_err(ParseGroError::BadResnum)
}
#[inline]
fn parse_resname(line: &str) -> Result<ResName> {
line.get(ranges::RESNAME)
.ok_or(ParseGroError::UnexpectedEOL(ExpectedItem::Resname))
.map(|s| s.trim().into())
}
#[inline]
fn parse_atomname(line: &str) -> Result<AtomName> {
line.get(ranges::ATOMNAME)
.ok_or(ParseGroError::UnexpectedEOL(ExpectedItem::Atomname))
.map(|s| s.trim().into())
}
#[inline]
fn parse_atomnum(line: &str) -> Result<AtomNum> {
line.get(ranges::ATOMNUM)
.ok_or(ParseGroError::UnexpectedEOL(ExpectedItem::Atomnum))?
.trim()
.parse()
.map_err(ParseGroError::BadAtomnum)
}
#[inline]
fn parse_position(line: &str) -> Result<[f32; 3]> {
parse_floats(
transpose_option_array([20..28, 28..36, 36..44].map(|r| line.get(r)))
.ok_or(ParseGroError::UnexpectedEOL(ExpectedItem::Position))?
.map(str::trim),
)
.map_err(ParseGroError::BadPositionValue)
}
#[inline]
fn parse_velocity(line: &str) -> Result<[f32; 3]> {
parse_floats(
transpose_option_array([44..52, 52..60, 60..68].map(|r| line.get(r)))
.ok_or(ParseGroError::UnexpectedEOL(ExpectedItem::Velocity))?
.map(str::trim),
)
.map_err(ParseGroError::BadVelocityValue)
}
#[rustfmt::skip]
#[inline]
fn parse_atom_line(line: &str) -> Result<A> {
assert!(line.len() >= ranges::MIN_LINE_LEN);
let resnum = if Self::PARSE_LIST.resnum { Some(Self::parse_resnum(line)?) } else { None };
let resname = if Self::PARSE_LIST.resname { Some(Self::parse_resname(line)?) } else { None };
let atomname = if Self::PARSE_LIST.atomname { Some(Self::parse_atomname(line)?) } else { None };
let atomnum = if Self::PARSE_LIST.atomnum { Some(Self::parse_atomnum(line)?) } else { None };
let position = if Self::PARSE_LIST.position { Some(Self::parse_position(line)?) } else { None };
let velocity = if Self::PARSE_LIST.velocity {
if line.len() > 44 {
Some(Self::parse_velocity(line)?)
} else {
None
}
} else {
None
};
Ok(Self::build_atom(
resnum, resname, atomname, atomnum, position, velocity,
))
}
fn parse_boxvecs(line: &str) -> Result<BoxVecs> {
let vs: Vec<&str> = line.split_ascii_whitespace().collect();
match vs.len() {
3 => Ok(BoxVecs::Short(
parse_floats(vs.try_into().unwrap()).map_err(ParseGroError::BadBoxVecsValue)?,
)),
9 => Ok(BoxVecs::Full(
parse_floats(vs.try_into().unwrap()).map_err(ParseGroError::BadBoxVecsValue)?,
)),
n => Err(ParseGroError::BadBoxVecsLength(n)),
}
}
fn read(reader: impl io::BufRead) -> Result<Self> {
let mut lines = reader.lines();
let title = String::from(
lines
.next()
.ok_or(ParseGroError::UnexpectedEOF(ExpectedItem::Title))??
.trim(),
);
let natoms = lines
.next()
.ok_or(ParseGroError::UnexpectedEOF(ExpectedItem::NAtoms))??
.trim()
.parse()
.map_err(ParseGroError::BadNAtoms)?;
let mut atoms = Vec::with_capacity(natoms);
for idx in 0..natoms {
let atom = Self::parse_atom_line(&lines.next().ok_or(
ParseGroError::InsufficientAtoms {
exp: natoms,
enc: idx,
},
)??)?;
atoms.push(atom);
}
let boxvecs = Self::parse_boxvecs(
&lines
.next()
.ok_or(ParseGroError::UnexpectedEOF(ExpectedItem::BoxVecs))??,
)?;
Ok(Self::build_structure(title, atoms, boxvecs))
}
fn read_from_file(file: std::fs::File) -> Result<Self> {
let reader = BufReader::new(file);
Self::read(reader)
}
fn open_gro<P: AsRef<Path>>(path: P) -> Result<Self> {
let file = std::fs::File::open(path)?;
Self::read_from_file(file)
}
}
impl ReadGro<Atom> for Structure {
#[inline]
fn build_atom(
resnum: Option<ResNum>,
resname: Option<ResName>,
atomname: Option<AtomName>,
atomnum: Option<AtomNum>,
position: Option<[f32; 3]>,
velocity: Option<[f32; 3]>,
) -> Atom {
Atom {
resnum: resnum.unwrap(),
resname: resname.unwrap(),
atomname: atomname.unwrap(),
atomnum: atomnum.unwrap(),
position: Vec3::from_array(position.unwrap()),
velocity: Vec3::from_array(velocity.unwrap_or_default()),
}
}
fn build_structure(title: String, atoms: Vec<Atom>, boxvecs: BoxVecs) -> Self {
Self {
title,
atoms,
boxvecs,
}
}
}
impl FromStr for Structure {
type Err = ParseGroError;
fn from_str(gro: &str) -> Result<Self> {
Self::read(gro.as_bytes())
}
}
fn transpose_option_array<const N: usize, T>(vs: [Option<T>; N]) -> Option<[T; N]> {
if vs.iter().any(Option::is_none) {
None
} else {
Some(vs.map(Option::unwrap))
}
}
#[inline]
fn parse_floats<const N: usize>(
vs: [&str; N],
) -> std::result::Result<[f32; N], std::num::ParseFloatError> {
vs.iter()
.map(|v| v.parse())
.collect::<std::result::Result<Vec<_>, _>>()
.map(|vs| vs.try_into().unwrap()) }
#[cfg(test)]
mod tests {
use std::io;
use super::*;
const EPS: f32 = 0.0001;
#[test]
fn open_gro() -> io::Result<()> {
let structure = Structure::open_gro(crate::tests::PATH)?;
assert_eq!(structure.title, "cg protein in water");
assert_eq!(structure.natoms(), 2869);
let center = Vec3::new(3.9875, 3.9760, 2.7035);
assert!(structure.center().abs_diff_eq(center, EPS));
assert_eq!(
structure.atoms[123].position,
Vec3::new(5.589, 5.256, 3.264)
);
assert_eq!(
structure.atoms[456].velocity,
Vec3::new(-0.144, -0.2159, -0.1276)
);
assert_eq!(
structure.boxvecs,
BoxVecs::Full([7.59661, 7.59661, 5.37162, 0.0, 0.0, 0.0, 0.0, 3.79831, 3.79831])
);
Ok(())
}
}