LOOS  v2.3.2
psf.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 #include <psf.hpp>
23 #include <exceptions.hpp>
24 
25 
26 namespace loos {
27 
28 
29  PSF* PSF::clone(void) const {
30  return(new PSF(*this));
31  }
32 
33  PSF PSF::copy(void) const {
34  AtomicGroup grp = this->AtomicGroup::copy();
35  PSF p(grp);
36 
37  // Add PSF specific member data copies here...
38  return(p);
39  }
40 
41 
42 
43  void PSF::read(std::istream& is) {
44  std::string input;
45 
46  // first line is the PSF header
47  if (!getline(is, input))
48  throw(FileReadError(_filename, "Failed reading first line of psf"));
49  if (input.substr(0,3) != "PSF")
50  throw(FileReadError(_filename, "PSF detected a non-PSF file"));
51 
52  // second line is blank
53  if (!getline(is, input))
54  throw(FileReadError(_filename, "PSF failed reading first header blank"));
55 
56  // third line is title header
57  getline(is, input);
58  int num_title_lines;
59  if (!(std::stringstream(input) >> num_title_lines))
60  throw(FileReadError(_filename, "PSF has malformed title header"));
61 
62  // skip the rest of the title
63  for (int i=0; i<num_title_lines; i++)
64  getline(is, input);
65 
66  // verify nothing went wrong
67  if (!(is.good()))
68  // Yes, I know, I should figure out what went wrong instead
69  // of running home crying. Sorry, Tod...
70  throw(FileReadError(_filename, "PSF choked reading the header"));
71 
72  // next line is blank
73  if (!getline(is, input))
74  throw(FileReadError(_filename, "PSF failed reading second header blank"));
75 
76  // next line is the number of atoms
77 
78  if (!getline(is, input))
79  throw(FileReadError(_filename, "PSF failed reading natom line"));
80  int num_atoms;
81  if (!(std::stringstream(input) >> num_atoms))
82  throw(FileReadError(_filename, "PSF has malformed natom line"));
83 
84  for (int i=0; i<num_atoms; i++) {
85  if (!getline(is, input)) {
86  std::ostringstream oss;
87  oss << "Failed reading PSF atom line for atom #" << (i+1);
88  throw(FileReadError(_filename, oss.str()));
89  }
90  parseAtomRecord(input);
91  }
92 
93  // next line is blank
94  if (!getline(is, input))
95  throw(FileReadError(_filename, "PSF failed reading blank after atom lines"));
96 
97  // next block of lines is the list of bonds
98  // Bond title line
99  if (!getline(is, input))
100  throw(FileReadError(_filename, "PSF failed reading nbond line"));
101  int num_bonds;
102  if (!(std::stringstream(input) >> num_bonds))
103  throw(FileReadError(_filename, "PSF has malformed nbond line"));
104 
105  int bonds_found = 0;
106  getline(is, input);
107  while (input.size() > 1) { // end of the block is marked by a blank line
108  // Note: >1 to handle \r in files that came from windows...
109  int ind1, ind2;
110  std::stringstream s(input);
111 
112  while (s.good()) {
113  if (!(s >> ind1 >> ind2))
114  throw(FileReadError(_filename, "PSF error parsing bonds.\n> " + input));
115 
116  if (ind1 > num_atoms || ind2 > num_atoms)
117  throw(FileReadError(_filename, "PSF bond error: bound atomid exceeds number of atoms.\n> " + input));
118 
119  ind1--; // our indices are 1 off from the numbering in the pdb/psf file
120  ind2--;
121  pAtom pa1 = getAtom(ind1);
122  pAtom pa2 = getAtom(ind2);
123  pa1->addBond(pa2);
124  pa2->addBond(pa1);
125  bonds_found++;
126 
127  // Catch returns in files that came from windows...
128  if (s.peek() == '\r')
129  break;
130  }
131  getline(is, input);
132  }
133  // sanity check
134  if (bonds_found != num_bonds)
135  throw(FileReadError(_filename, "PSF number of bonds disagrees with number found"));
136 
138  setGroupConnectivity();
139  }
140 
141 
142 
143  void PSF::parseAtomRecord(const std::string s) {
144  gint index;
145  std::string segname;
146  gint resid;
147  std::string resname;
148  std::string atomname;
149  std::string atomtype;
150  greal charge;
151  greal mass;
152  gint fixed;
153 
154  pAtom pa(new Atom);
155  pa->index(_max_index++);
156 
157  std::stringstream ss(s);
158 
159  if (!(ss >> index))
160  throw(FileReadError(_filename, "PSF parse error.\n> " + s));
161  pa->id(index);
162 
163  if (!(ss >> segname))
164  throw(FileReadError(_filename, "PSF parse error.\n> " + s));
165  pa->segid(segname);
166 
167 
168  if (!(ss >> resid))
169  throw(FileReadError(_filename, "PSF parse error.\n> " + s));
170  pa->resid(resid);
171 
172  if (!(ss >> resname))
173  throw(FileReadError(_filename, "PSF parse error.\n> " + s));
174  pa->resname(resname);
175 
176  if (!(ss >> atomname))
177  throw(FileReadError(_filename, "PSF parse error.\n> " + s));
178  pa->name(atomname);
179 
180  // If this is a charmm psf, the atomtype will be an integer.
181  // NAMD/XPLOR psfs use the symbolic atomtype, which must start with a letter
182  // At the moment, the Atom class doesn't care about this value (it's mostly
183  // used in charmm and namd as a means to look up parameters), so we're going to
184  // discard it. However, if we ever decide we're going to use this, we'll need
185  // to keep track of the distinction between charmm and namd usage.
186  if (!(ss >> atomtype))
187  throw(FileReadError(_filename, "PSF parse error.\n> " + s));
188 
189  if (!(ss >> charge))
190  throw(FileReadError(_filename, "PSF parse error.\n> " + s));
191  pa->charge(charge);
192 
193  if (!(ss >> mass))
194  throw(FileReadError(_filename, "PSF parse error.\n> " + s));
195  pa->mass(mass);
196 
197  // Is the atom fixed or mobile?
198  // for now, we're going to silently drop this
199  ss >> fixed;
200 
201  append(pa);
202  }
203 
204 }
pAtom getAtom(const int i) const
Get the ith atom from this group.
Definition: AtomicGroup.cpp:83
Errors that occur while reading a file.
Definition: exceptions.hpp:180
AtomicGroup & append(pAtom pa)
Append the atom onto the group.
virtual PSF * clone(void) const
Clones an object for polymorphism (see AtomicGroup::clone() for more info)
Definition: psf.cpp:29
PSF copy(void) const
Creates a deep copy (see AtomicGroup::copy() for more info)
Definition: psf.cpp:33
AtomicGroup copy(void) const
Creates a deep copy of this group.
Definition: AtomicGroup.cpp:56
uint deduceAtomicNumberFromMass(const double tol=0.1)
Deduce atomic number from mass (if present), returning number of atoms assigned.
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
Namespace for most things not already encapsulated within a class.
Class for reading a subset of the PSF format.
Definition: psf.hpp:60