LOOS  v2.3.2
hcore.cpp
1 /*
2  Core code for hbonds tools
3 */
4 
5 
6 /*
7 
8  This file is part of LOOS.
9 
10  LOOS (Lightweight Object-Oriented Structure library)
11  Copyright (c) 2010, Tod D. Romo
12  Department of Biochemistry and Biophysics
13  School of Medicine & Dentistry, University of Rochester
14 
15  This package (LOOS) is free software: you can redistribute it and/or modify
16  it under the terms of the GNU General Public License as published by
17  the Free Software Foundation under version 3 of the License.
18 
19  This package is distributed in the hope that it will be useful,
20  but WITHOUT ANY WARRANTY; without even the implied warranty of
21  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  GNU General Public License for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with this program. If not, see <http://www.gnu.org/licenses/>.
26 */
27 
28 
29 
30 
31 #include <boost/format.hpp>
32 
33 #include "hcore.hpp"
34 
35 using namespace loos;
36 using namespace HBonds;
37 
38 double SimpleAtom::inner = 0.0;
39 double SimpleAtom::outer = 3.5;
40 double SimpleAtom::deviation = 20.0;
41 
42 
43 // Reports distance^2 between hydrogen and heavy atom
44 // D-H ... X
45 // |-----|
46 double SimpleAtom::distance2(const SimpleAtom& s) const
47 {
48  double d;
49 
50  if (usePeriodicity)
51  d = atom->coords().distance2(s.atom->coords(), sbox.box());
52  else
53  d = atom->coords().distance2(s.atom->coords());
54 
55  return(d);
56 }
57 
58 
59 // Returns angle between atoms in degrees...
60 //
61 // D-H ... X
62 // \---/
63 double SimpleAtom::angle(const SimpleAtom& s) const {
64  loos::GCoord left, middle, right;
65 
66  if (isHydrogen) {
67  if (s.isHydrogen)
68  throw(std::runtime_error("Cannot take the angle between two hydrogens"));
69 
70  left = attached_to->coords();
71  middle = atom->coords();
72  right = s.atom->coords();
73 
74  } else {
75  if (!s.isHydrogen)
76  throw(std::runtime_error("Cannot take the angle between two non-hydrogens"));
77 
78  left = atom->coords();
79  middle = s.atom->coords();
80  right = s.attached_to->coords();
81  }
82 
83  if (usePeriodicity) {
84  loos::GCoord box = sbox.box();
85 
86  left.reimage(box);
87  middle.reimage(box);
88  right.reimage(box);
89  }
90 
91  return(loos::Math::angle(left, middle, right));
92 }
93 
94 
95 
96 
97 // Converts a selection into a vector of SimpleAtoms. All new
98 // SimpleAtom's get assigned the use_periodicity flag
99 
100 std::vector<SimpleAtom> SimpleAtom::processSelection(const std::string& selection, const loos::AtomicGroup& system, const bool use_periodicity) {
101  std::vector<SimpleAtom> processed;
102 
103  // We don't want to force model to be sorted (esp since it's const)
104  // So, create a lightweight copy but use the same pAtoms...
105  loos::AtomicGroup searchable(system);
106  searchable.sort();
107 
108  loos::AtomicGroup model = selectAtoms(system, selection);
109 
110  for (loos::AtomicGroup::const_iterator i= model.begin(); i != model.end(); ++i) {
111  SimpleAtom new_atom(*i, system.sharedPeriodicBox(), use_periodicity);
112  std::string n = (*i)->name();
113  if (n[0] == 'H') {
114  new_atom.isHydrogen = true;
115  std::vector<int> bond_list = (*i)->getBonds();
116  if (bond_list.size() == 0)
117  throw(ErrorWithAtom(*i, "Detected a hydrogen bond that has no connectivity"));
118  if (bond_list.size() != 1)
119  throw(ErrorWithAtom(*i, "Detected a hydrogen bond that has more than one atom bound"));
120 
121  loos::pAtom pa = searchable.findById(bond_list[0]);
122  if (pa == 0)
123  throw(ErrorWithAtom(*i, "Cannot find atom hydrogen is bound too"));
124 
125  new_atom.attached_to = pa;
126  }
127 
128  processed.push_back(new_atom);
129  }
130 
131  return(processed);
132 }
133 
134 
135 
136 
137 // Check for a hdyrogen bond between two SimpleAtom's
138 
139 bool SimpleAtom::hydrogenBond(const SAtom& o) const {
140  double dist = distance2(o);
141  double angl = angle(o);
142 
143  if (dist >= inner && dist <= outer && fmod(fabs(angl - 180.0), 360.0) <= deviation)
144  return(true);
145 
146  return(false);
147 }
148 
149 
150 
151 // Returns an AtomicGroup containing the atoms that are hydrogen
152 // bonded to self. If findFirstOnly is true, than the first hydrogen
153 // bond found causes the function to return... (i.e. it may be a
154 // small optimization in performance)
155 
156 loos::AtomicGroup SimpleAtom::findHydrogenBonds(const std::vector<SimpleAtom>& group, const bool findFirstOnly) {
157  loos::AtomicGroup results;
158 
159  for (SAGroup::const_iterator i = group.begin(); i != group.end(); ++i) {
160  bool b = hydrogenBond(*i);
161  if (b) {
162  results.append(i->atom);
163  if (findFirstOnly)
164  break;
165  }
166  }
167 
168  return(results);
169 }
170 
171 
172 
173 // Returns a vector of flags indicating which SimpleAtoms form a
174 // hydrogen bond to self.
175 
176 std::vector<uint> SimpleAtom::findHydrogenBondsVector(const std::vector<SimpleAtom>& group) {
177  std::vector<uint> results;
178 
179  for (SAGroup::const_iterator i = group.begin(); i != group.end(); ++i)
180  results.push_back(hydrogenBond(*i));
181 
182  return(results);
183 }
184 
185 
186 // Returns a matrix of flags indicating which SimpleAtoms form a
187 // hydrogen bond to self of a trajectory
188 
189 BondMatrix SimpleAtom::findHydrogenBondsMatrix(const std::vector<SimpleAtom>& group, loos::pTraj& traj, loos::AtomicGroup& model, const uint maxt) const {
190 
191  if (maxt > traj->nframes()) {
192  std::cerr << boost::format("Error- row clip (%d) exceeds trajectory size (%d)\n") % maxt % traj->nframes();
193  exit(-10);
194  }
195 
196  BondMatrix bonds(maxt, group.size());
197 
198  for (uint t = 0; t < maxt; ++t) {
199  traj->readFrame(t);
200  traj->updateGroupCoords(model);
201 
202  for (uint i = 0; i<group.size(); ++i)
203  bonds(t, i) = hydrogenBond(group[i]);
204 
205  }
206 
207  return(bonds);
208 }
209 
210 
211 
212 
213 bool SimpleAtom::divineHydrogen(const std::string& name) {
214  if (name[0] == 'H')
215  return(true);
216 
217  return(false);
218 }
Simple matrix template class using policy classes to determine behavior.
Definition: MatrixImpl.hpp:53
AtomicGroup & append(pAtom pa)
Append the atom onto the group.
loos::SharedPeriodicBox sharedPeriodicBox() const
Provide access to the underlying shared periodic box...
greal angle(const GCoord &, const GCoord &, const GCoord &, const GCoord *=NULL)
Compute the angle in degrees assuming the middle is the vertex.
Definition: Geometry.cpp:36
AtomicGroup selectAtoms(const AtomicGroup &source, const std::string selection)
Applies a string-based selection to an atomic group...
Definition: utils.cpp:195
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
Namespace for most things not already encapsulated within a class.
void reimage(const Coord< T > &box)
Handle coordinates with periodic boundary conditions.
Definition: Coord.hpp:404