LOOS  v2.3.2
amber.cpp
1 /*
2  This file is part of LOOS.
3 
4  LOOS (Lightweight Object-Oriented Structure library)
5  Copyright (c) 2008, Tod D. Romo, Alan Grossfield
6  Department of Biochemistry and Biophysics
7  School of Medicine & Dentistry, University of Rochester
8 
9  This package (LOOS) is free software: you can redistribute it and/or modify
10  it under the terms of the GNU General Public License as published by
11  the Free Software Foundation under version 3 of the License.
12 
13  This package is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with this program. If not, see <http://www.gnu.org/licenses/>.
20 */
21 
22 
23 #include <amber.hpp>
24 #include <exceptions.hpp>
25 #include <stdlib.h>
26 #include <ctype.h>
27 
28 #include <iomanip>
29 
30 namespace loos {
31 
32 
33  // Parse simple Fortran format specifications, extracted from a %FORMAT tag...
34  // Takes a string of expected format types (characters) that the extracted format
35  // is compared against. For example, to parse floats, expected could be any of "FEG".
36  // If the format type does not match (for example, say "%FORMAT (20I8)" when expected
37  // is "FEG"), then an error will be thrown.
38  //
39  // Format specs are converted to upper-case for validation
40 
41  Amber::FormatSpec Amber::parseFormat(const std::string& expected_types, const std::string& where) {
42 
43  reader.getNext();
44  char period;
45  FormatSpec fmt;
46 
47  // Verify line has a %FORMAT tag
48  if (reader.line().compare(0, 7, "%FORMAT") != 0)
49  throw(FileReadErrorWithLine(reader.name(), "Expected format for " + where, reader.lineNumber()));
50 
51  // Extract format spec between parens...
52  boost::char_separator<char> sep("()");
53  std::string input_line(reader.line()); // Ubuntu 16.04 apparently requires that the string in tokenizer
54  tokenizer tokens(input_line, sep); // be non-const
55  tokenizer::iterator toks = tokens.begin();
56 
57  ++toks;
58  std::istringstream iss(*toks);
59 
60  // try nXw.d first
61  if (! (iss >> fmt.repeat >> fmt.type >> fmt.width >> period >> fmt.precision) ) {
62  iss.clear();
63  iss.seekg(0);
64  // Now try nXw
65  if (! (iss >> fmt.repeat >> fmt.type >> fmt.width) ) {
66  iss.clear();
67  iss.seekg(0);
68  // Now try Xw
69  if (! (iss >> fmt.type >> fmt.width) ) {
70  iss.clear();
71  iss.seekg(0);
72  // And finally just try X
73  if (! (iss >> fmt.type) )
74  throw(FileReadErrorWithLine(reader.name(), "Cannot parse format for " + where, reader.lineNumber()));
75  }
76  }
77  }
78 
79  // Compare against expected types
80  std::string expected_types_UC = boost::to_upper_copy(expected_types);
81 
82  if (expected_types_UC.find_first_of(toupper(fmt.type)) == std::string::npos)
83  throw(FileReadErrorWithLine(reader.name(), "Invalid format type for " + where, reader.lineNumber()));
84 
85  return(fmt);
86 
87  }
88 
89 
90  void Amber::parseCharges() {
91  FormatSpec fmt = parseFormat("EFG", "charges");
92 
93  std::vector<double> charges = readBlock<double>(fmt.width);
94  if (charges.size() != atoms.size())
95  throw(FileReadErrorWithLine(reader.name(), "Error parsing charges from amber file", reader.lineNumber()));
96 
97  for (uint i=0; i<charges.size(); ++i)
98  atoms[i]->charge(charges[i]);
99  }
100 
101 
102  void Amber::parseMasses() {
103  FormatSpec fmt = parseFormat("EFG", "masses");
104 
105  std::vector<double> masses = readBlock<double>(fmt.width);
106  if (masses.size() != atoms.size())
107  throw(FileReadErrorWithLine(reader.name(), "Error parsing masses from amber file", reader.lineNumber()));
108 
109  for (uint i=0; i<masses.size(); ++i)
110  atoms[i]->mass(masses[i]);
111 
112  }
113 
114 
115 
116  void Amber::parseResidueLabels() {
117  FormatSpec fmt = parseFormat("a", "residue labels");
118 
119  std::vector<std::string> labels = readBlock<std::string>(fmt.width);
120  if (labels.size() != nres)
121  throw(FileReadErrorWithLine(reader.name(), "Error parsing residue labels from amber file", reader.lineNumber()));
122 
123  residue_labels = labels;
124  }
125 
126 
127  void Amber::parseResiduePointers() {
128  FormatSpec fmt = parseFormat("I", "residue pointers");
129 
130  std::vector<uint> pointers = readBlock<uint>(fmt.width);
131  if (pointers.size() != nres)
132  throw(FileReadErrorWithLine(reader.name(), "Error parsing residue pointers from amber file", reader.lineNumber()));
133  residue_pointers = pointers;
134  }
135 
136 
137  void Amber::assignResidues(void) {
138  if (!(residue_pointers.size() == nres && residue_labels.size() == nres))
139  throw(std::runtime_error("Unable to assign residues."));
140 
141  int curresid = 0;
142  uint i = 0;
143  std::string curresname;
144 
145  for (i=0; i<nres-1; i++) {
146  ++curresid;
147  curresname = residue_labels[i];
148  for (uint j = residue_pointers[i]; j<residue_pointers[i+1]; j++) {
149  atoms[j-1]->resid(curresid);
150  atoms[j-1]->resname(curresname);
151  }
152  }
153 
154  // Fix the end...
155  ++curresid;
156  curresname = residue_labels[i];
157  for (uint j = residue_pointers[i]-1; j<natoms; j++) {
158  atoms[j]->resid(curresid);
159  atoms[j]->resname(curresname);
160  }
161  }
162 
163 
164 
165 
166  void Amber::parseBonds(const uint n) {
167  FormatSpec fmt = parseFormat("I", "bonds");
168 
169 
170  std::vector<int> bond_list = readBlock<int>(fmt.width);
171 
172  if (bond_list.size() != 3*n)
173  throw(FileReadErrorWithLine(reader.name(), "Error parsing bonds in amber file", reader.lineNumber()));
174 
175  for (uint i=0; i<bond_list.size(); i += 3) {
176  if (bond_list[i] == bond_list[i+1])
177  continue;
178 
179  pAtom aatom = atoms[bond_list[i]/3];
180  pAtom batom = atoms[bond_list[i+1]/3];
181 
182  // Amber bond lists are not symmetric, so make sure we add both pairs...
183  if (!(aatom->isBoundTo(batom)))
184  aatom->addBond(batom);
185  if (!(batom->isBoundTo(aatom)))
186  batom->addBond(aatom);
187 
188  }
189 
190  }
191 
192 
193  void Amber::parsePointers() {
194  FormatSpec fmt = parseFormat("I", "pointers");
195 
196  std::vector<uint> pointers = readBlock<uint>(fmt.width);
197 
198  // Now build up the atomic-group...
199  if (atoms.size() != 0)
200  throw(std::logic_error("Internal error: trying to read in an amber parmtop into a non-empty group!"));
201 
202 
203  natoms = pointers[0];
204  nbonh = pointers[2];
205  mbona = pointers[3];
206  nres = pointers[11];
207 
208  for (uint i=0; i<natoms; i++) {
209  pAtom pa(new Atom);
210  pa->id(i+1);
211  pa->index(i);
212  atoms.push_back(pa);
213  }
214 
215  }
216 
217 
218  // Simply slurp up the title (for now)
219  void Amber::parseTitle() {
220 
221  FormatSpec fmt = parseFormat("a", "title");
222  std::vector<std::string> titles = readBlock<std::string>(fmt.width);
223  for (std::vector<std::string>::const_iterator i = titles.begin(); i != titles.end(); ++i)
224  _title += *i;
225  }
226 
227 
228  void Amber::parseAtomNames() {
229  FormatSpec fmt = parseFormat("a", "atom names");
230 
231  std::vector<std::string> names = readBlock<std::string>(fmt.width);
232  if (names.size() != natoms)
233  throw(FileReadErrorWithLine(reader.name(), "Error parsing atom names", reader.lineNumber()));
234  for (uint i=0; i<names.size(); ++i)
235  atoms[i]->name(names[i]);
236  }
237 
238  void Amber::parseAmoebaRegularBondNumList() {
239  FormatSpec fmt = parseFormat("I", "amoeba_regular_num_bond_list");
240  reader.getNext();
241  std::istringstream iss(reader.line());
242 
243  if (! (iss >> std::setw(fmt.width) >> _amoeba_regular_bond_num_list))
244  throw(FileReadErrorWithLine(reader.name(), "Error parsing amoeba_regular_bond_num_list", reader.lineNumber()));
245  }
246 
247 
248 
249  void Amber::parseAmoebaRegularBondList(const uint n) {
250  FormatSpec fmt = parseFormat("I", "amoeba_regular_bond_list");
251 
252 
253  std::vector<int> bond_list = readBlock<int>(fmt.width);
254 
255  if (bond_list.size() != 3*n)
256  throw(FileReadErrorWithLine(reader.name(), "Error parsing amoeba bonds in amber file", reader.lineNumber()));
257 
258  for (uint i=0; i<bond_list.size(); i += 3) {
259  if (bond_list[i] == bond_list[i+1])
260  continue;
261 
262  // Amoeba bond indices appear not to be /3 as regular amber bonds are...
263  // Are we sure???
264  pAtom aatom = atoms[bond_list[i]-1];
265  pAtom batom = atoms[bond_list[i+1]-1];
266 
267  // Amber bond lists are not symmetric, so make sure we add both pairs...
268  if (!(aatom->isBoundTo(batom)))
269  aatom->addBond(batom);
270  if (!(batom->isBoundTo(aatom)))
271  batom->addBond(aatom);
272 
273  }
274 
275  }
276 
277 
278 
279  void Amber::read(std::istream& ifs) {
280  reader.stream(ifs);
281 
282  while (reader.getNext()) {
283 
284  boost::char_separator<char> sep(" \t");
285  std::string input_line(reader.line()); // See parseFormat() for explanation...required by Ubuntu 16
286  tokenizer tokens(input_line, sep);
287  tokenizer::iterator toks = tokens.begin();
288 
289  if (*toks != "%FLAG")
290  continue;
291 
292  ++toks;
293  if (*toks == "TITLE")
294  parseTitle();
295  else if (*toks == "POINTERS")
296  parsePointers();
297  else if (*toks == "ATOM_NAME")
298  parseAtomNames();
299  else if (*toks == "CHARGE")
300  parseCharges();
301  else if (*toks == "MASS")
302  parseMasses();
303  else if (*toks == "RESIDUE_LABEL")
304  parseResidueLabels();
305  else if (*toks == "RESIDUE_POINTER")
306  parseResiduePointers();
307  else if (*toks == "BONDS_INC_HYDROGEN")
308  parseBonds(nbonh);
309  else if (*toks == "BONDS_WITHOUT_HYDROGEN")
310  parseBonds( mbona);
311  else if (*toks == "AMOEBA_REGULAR_BOND_NUM_LIST")
312  parseAmoebaRegularBondNumList();
313  else if (*toks == "AMOEBA_REGULAR_BOND_LIST")
314  parseAmoebaRegularBondList(_amoeba_regular_bond_num_list);
315 
316  }
317 
318  assignResidues();
320  setGroupConnectivity();
321  }
322 
323 
324 }
325 
void read(std::istream &ifs)
Parse the parmtop file.
Definition: amber.cpp:279
uint deduceAtomicNumberFromMass(const double tol=0.1)
Deduce atomic number from mass (if present), returning number of atoms assigned.
Namespace for most things not already encapsulated within a class.