heme/
io.rs

1// File input and output
2// WARNING: contains hand-rolled PDB parser
3
4use conformation::{XYZ, Atom, Pose};
5use sampling::get_protocol;
6
7use std::fs::File;
8use std::error::Error;
9use std::io::prelude::*;
10
11pub struct Config {
12    // a Config objects holds the  config options
13    pub protocol: String,
14    pub filename: String,
15}
16
17impl Config {
18    pub fn new(args: &[String]) -> Result<Config, &'static str> {
19        if args.len() < 3 {
20            return Err("not enough arguments");
21        }
22
23        let protocol = args[1].clone();
24        let filename = args[2].clone();
25
26        Ok(Config { protocol, filename })
27    }
28}
29
30// PDB parsing
31
32// This time, and this one time *only*
33// we are going to implement a PDB parser.
34// Literally as soon as we have Internet, we
35// will look for an existing PDB parser.
36// I promise. -- Alex (Jul 30, 2017 at 7:30 PM)
37
38// 0         1         2         3         4         5         6
39// 01234567890123456789012345678901234567890123456789012345678901234567890123456789
40// ATOM      1  N   ASP A   1      27.405  -7.086  18.389  1.00  0.00           N
41// ATOM      2  CA  ASP A   1      26.221  -7.728  18.989  1.00  0.00           C
42//                               ||||||||        ||||||||
43//                                       ||||||||
44//                               30..38  38..46  46..44
45
46pub struct Record {
47    // PDB atom record class
48    // and this is the only struct
49    // we're gonna make
50
51    pub xyz: XYZ,
52    pub charge: f64,
53    pub element: String,
54    pub residue_index: i32,
55    pub residue_name: String,
56}
57
58impl Record {
59    pub fn new(args: &str) -> Record {
60
61        // coords
62        let x: f64 = args[30..38].trim().parse().expect("X not a number");
63        let y: f64 = args[38..46].trim().parse().expect("Y not a number");
64        let z: f64 = args[46..54].trim().parse().expect("Z not a number");
65        let xyz = XYZ::new(x, y, z);
66
67        // residue-level information
68        let residue_index: i32 = args[23..30].trim().parse().expect("Residue index not a number");
69        let residue_name = args[17..20].to_string(); //String::from(args[17..20]);
70
71        // charge and radius
72        let charge: f64 = -0.69; // actually look up
73        let element = args[11..17].trim().to_string();
74        // let charge = look_up_charge(atom_name); or something
75        //println!("Creating atom with {}, {}, {}, charge {}", x, y, z, charge);
76
77        Record { xyz, charge, element, residue_index, residue_name }
78    }
79}
80
81pub fn parse_pdb(contents: &str) -> Vec<Atom> {
82    // Parse a PDB file
83
84    let mut atoms = Vec::new();
85    // let mut residues = Vec::new();
86    for line in contents.lines() {
87        if line.starts_with("ATOM") {
88            let record = Record::new(&line);
89            let atom = Atom::new(&record);
90            atoms.push(atom);
91        }
92    }
93
94    atoms
95}
96
97pub fn run(config: Config) -> Result<(), Box<Error>> {
98
99    // read input files from Config object
100    let mut f = File::open(config.filename)?;
101    let mut contents = String::new();
102    f.read_to_string(&mut contents)?;
103
104    // create pose object by parsing the PDB
105    let atoms = parse_pdb(&contents);
106    let mut pose = Pose::from_atoms(atoms);
107
108    // Apply a protocol to the Pose
109    let protocol = get_protocol(&config.protocol);
110    protocol.run(&mut pose);
111
112    // println!("Total score: {}", total_score);
113
114    Ok(())
115}
116
117#[cfg(test)]
118mod tests {
119    use super::*;
120
121    #[test]
122    fn test_new_atoms_testy() {
123        let record = Record::new("ATOM      8  OD2 ASP A   1      28.042 -10.262  20.858  1.00  0.00           O ");
124        let atom = Atom::new(&record);
125        assert_eq!(atom.xyz.x, 28.042);
126    }
127
128    #[test]
129    fn test_pdb_parsing() {
130        let contents = "ATOM      1  N   ASP A   1      27.405  -7.086  18.389  1.00  0.00           N
131ATOM      2  CA  ASP A   1      26.221  -7.728  18.989  1.00  0.00           C
132ATOM      3  C   ASP A   1      25.150  -6.679  19.328  1.00  0.00           C
133ATOM      4  O   ASP A   1      24.099  -6.710  18.722  1.00  0.00           O
134ATOM      5  CB  ASP A   1      26.607  -8.504  20.250  1.00  0.00           C
135ATOM      6  CG  ASP A   1      27.438  -9.745  19.951  1.00  0.00           C
136ATOM      7  OD1 ASP A   1      27.457 -10.165  18.817  1.00  0.00           O
137ATOM      8  OD2 ASP A   1      28.042 -10.262  20.858  1.00  0.00           O
138ATOM      9 1H   ASP A   1      28.091  -7.781  18.175  1.00  0.00           H
139ATOM     10 2H   ASP A   1      27.138  -6.611  17.550  1.00  0.00           H
140";
141        let atoms = parse_pdb(&contents);
142        let pose = Pose::from_atoms(atoms);
143
144        assert_eq!(pose.atoms.len(), 10);
145
146    }
147
148    #[test]
149    fn test_entropy_function() {
150
151        // entropy function
152        fn H(l: &f64) -> f64 {
153
154            // bounds of N choose lN
155            // can be found by 2^(n*H(l))
156            // where H is this function
157
158            // 0 < l < 0
159
160            let value = 1.0 - l;
161            l.ln() - value * value.ln()
162        }
163
164        let test_float: f64 = 18.11;
165        let answer: f64 = H(&test_float);
166        println!(">>> {}", answer);
167
168        // assert_eq!(4, 4);
169    }
170
171    #[test]
172    fn test_lennard_jones_atom_loading() {
173
174
175    }
176
177}
heme/io.rs

heme/
io.rs
