LOOS  v2.3.2
hbonds.cpp
1 /*
2  hbonds
3 
4  Find putative hydrogen-bonds based on user-specified criteria (angle and distance)
5 */
6 
7 
8 
9 /*
10 
11  This file is part of LOOS.
12 
13  LOOS (Lightweight Object-Oriented Structure library)
14  Copyright (c) 2010, Tod D. Romo
15  Department of Biochemistry and Biophysics
16  School of Medicine & Dentistry, University of Rochester
17 
18  This package (LOOS) is free software: you can redistribute it and/or modify
19  it under the terms of the GNU General Public License as published by
20  the Free Software Foundation under version 3 of the License.
21 
22  This package is distributed in the hope that it will be useful,
23  but WITHOUT ANY WARRANTY; without even the implied warranty of
24  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  GNU General Public License for more details.
26 
27  You should have received a copy of the GNU General Public License
28  along with this program. If not, see <http://www.gnu.org/licenses/>.
29 */
30 
31 
32 #include <loos.hpp>
33 #include <boost/format.hpp>
34 #include <boost/program_options.hpp>
35 
36 #include "hcore.hpp"
37 
38 using namespace std;
39 using namespace loos;
40 using namespace HBonds;
41 namespace po = boost::program_options;
42 namespace opts = loos::OptionsFramework;
43 
44 
45 typedef vector<string> veString;
46 typedef vector<double> veDouble;
47 typedef vector<veDouble> veveDouble;
48 typedef vector<uint> veUint;
49 typedef vector<veUint> veveUint;
50 
51 
53 
54 
55 // --------------- GLOBALS
56 
57 double length_low, length_high;
58 double max_angle;
59 bool use_periodicity;
60 bool verbose;
61 bool use_stderr;
62 vector<string> acceptor_names;
63 vector<string> acceptor_selections;
64 string donor_selection;
65 string model_name;
66 vector<string> traj_names;
67 
68 uint skip = 0;
69 
70 
71 // ---------------
72 // @cond TOOLS_INTERNAL
73 
74 
75 string fullHelpMessage(void) {
76  string msg =
77  "\n"
78  "SYNOPSIS\n"
79  "\tHydrogen bond occupancy for a trajectory\n"
80  "\n"
81  "DESCRIPTION\n"
82  "\n"
83  "\tThis tool computes the occupancy of putative hydrogen bonds (defined by\n"
84  "a simple distance and angle criteria). The 'donor' selection must have one\n"
85  "hydrogen present and the 'acceptor' should have no hydrogens. Multiple acceptors\n"
86  "may be given on the command line. These are specified by using multiple sets of\n"
87  "options, i.e. -N name -S selection where name is the label for the acceptor and\n"
88  "selection is the corresponding LOOS selection string. There must be at least\n"
89  "one name/selection pair. The occupancy calculation can also be performed over\n"
90  "multiple trajectories by specifying more than one on the command line.\n"
91  "\n"
92  "EXAMPLES\n"
93  "\n"
94  "\thbonds -N 'Carbonyl' -S 'name == \"O1\" && resname == \"PALM\"' \\\n"
95  "\t 'resid == 4 && name == \"HE1\"' model.psf traj.dcd\n"
96  "This example uses the palmitoyl carbonyl oxygen as the acceptor and the HE1 hydrogen from\n"
97  "residue 4 as the donor.\n"
98  "\n"
99  "\thbonds -N 'Carbonyl' -S 'name == \"O1\" && resname == \"PALM\"' \\\n"
100  "\t -N 'Phosphate' -S 'name == \"OP1\" && resname == \"PEGL\"' \\\n"
101  "\t 'resid == 4 && name == \"HE1\"' model.psf traj.dcd\n"
102  "This example uses the palmitoyl carbonyl oxygen as above, but also looks for hydrogen\n"
103  "bonds with the OP1 phosphate oxygen in residue PEGL. The same donor as above is used.\n"
104  "\n"
105  "\thbonds --blow 2 --bhi 4 --angle 20 -N 'Carbonyl' \\\n"
106  "\t -S 'name == \"O1\" && resname == \"PALM\"' 'resid == 4 && name == \"HE1\"' \\\n"
107  "\t model.psf traj.dcd\n"
108  "This example is the same as the first, however the criteria for hydrogen bonds are now\n"
109  "that they cannot be shorter than 2 angstroms nor longer than 4 angstroms, and the angle\n"
110  "cannot be more than 20 degrees from linear.\n"
111  "\n"
112  "SEE ALSO\n"
113  "\thmatrix, hcorrelation\n";
114 
115  return(msg);
116 }
117 
118 
119 
120 class ToolOptions : public opts::OptionsPackage {
121 public:
122  void addGeneric(po::options_description& o) {
123  o.add_options()
124  ("skip,k", po::value<uint>(&skip)->default_value(0), "Number of frames to skip")
125  ("stderr", po::value<bool>(&use_stderr)->default_value(false), "Report stderr rather than stddev")
126  ("blow", po::value<double>(&length_low)->default_value(1.5), "Low cutoff for bond length")
127  ("bhi", po::value<double>(&length_high)->default_value(3.0), "High cutoff for bond length")
128  ("angle", po::value<double>(&max_angle)->default_value(30.0), "Max bond angle deviation from linear")
129  ("periodic", po::value<bool>(&use_periodicity)->default_value(false), "Use periodic boundary")
130  ("name,N", po::value< vector<string> >(&acceptor_names), "Name of an acceptor selection (required)")
131  ("acceptor,S", po::value< vector<string> >(&acceptor_selections), "Acceptor selection (required)");
132  }
133 
134  void addHidden(po::options_description& o) {
135  o.add_options()
136  ("donor", po::value<string>(&donor_selection), "donor selection")
137  ("model", po::value<string>(&model_name), "model")
138  ("trajs", po::value< vector<string> >(&traj_names), "Trajectories");
139  }
140 
141  void addPositional(po::positional_options_description& opts) {
142  opts.add("donor", 1);
143  opts.add("model", 1);
144  opts.add("trajs", -1);
145  }
146 
147  bool postConditions(po::variables_map& map) {
148  if (acceptor_selections.empty()) {
149  cerr << "Error- must provide at least one acceptor name and selection.\n";
150  return(false);
151  }
152  if (acceptor_selections.size() != acceptor_names.size()) {
153  cerr << "Error- must provide one name for each acceptor selection.\n";
154  return(false);
155  }
156 
157  return(true);
158  }
159 
160  string help() const {
161  return("donor model traj [traj ...]");
162  }
163 
164  string print() const {
165  ostringstream oss;
166  oss << boost::format("skip=%d,stderr=%d,blow=%f,bhi=%f,angle=%f,periodic=%d,names=\"%s\",acceptors=\"%s\",donor=\"%s\",model=\"%s\",trajs=\"%s\"")
167  % skip
168  % use_stderr
169  % length_low
170  % length_high
171  % max_angle
172  % use_periodicity
173  % vectorAsStringWithCommas(acceptor_names)
174  % vectorAsStringWithCommas(acceptor_selections)
175  % donor_selection
176  % model_name
177  % vectorAsStringWithCommas(traj_names);
178 
179  return(oss.str());
180  }
181 
182 };
183 
184 
185 
186 
187 // @endcond
188 
189 
190 veDouble rowAverage(const Matrix& M) {
191  uint m = M.rows();
192  uint n = M.cols();
193 
194  veDouble avg;
195 
196  for (uint j=0; j<m; ++j) {
197  double x = 0.0;
198  for (uint i=0; i<n; ++i)
199  x += M(j, i);
200  avg.push_back(x/n);
201  }
202 
203  return(avg);
204 }
205 
206 
207 veDouble rowStd(const Matrix& M, const veDouble& avg) {
208  uint m = M.rows();
209  uint n = M.cols();
210 
211 
212  if (n < 3)
213  return(veDouble(m, 0.0));
214 
215  veDouble stddev;
216 
217  for (uint j=0; j<m; ++j) {
218  double x = 0.0;
219  for (uint i=0; i<n; ++i) {
220  double d = M(j, i) - avg[j];
221  x += d*d;
222  }
223  stddev.push_back(sqrt(x/(n-1.0)));
224  }
225 
226  return(stddev);
227 }
228 
229 
230 
231 int main(int argc, char *argv[]) {
232  string hdr = invocationHeader(argc, argv);
233 
234  opts::BasicOptions* bopts = new opts::BasicOptions(fullHelpMessage());
235  ToolOptions* topts = new ToolOptions;
236 
237  opts::AggregateOptions options;
238  options.add(bopts).add(topts);
239  if (!options.parse(argc, argv))
240  exit(-1);
241 
242  SimpleAtom::innerRadius(length_low);
243  SimpleAtom::outerRadius(length_high);
244  SimpleAtom::maxDeviation(max_angle);
245 
246  cout << "# " << hdr << endl;
247  cout << "# " << vectorAsStringWithCommas(options.print()) << endl;
248 
249  AtomicGroup model = createSystem(model_name);
250 
251  SAGroup donors = SimpleAtom::processSelection(donor_selection, model, use_periodicity);
252 
253  vector<SAGroup> acceptors;
254  for (uint i=0; i<acceptor_selections.size(); ++i) {
255  SAGroup acceptor = SimpleAtom::processSelection(acceptor_selections[i], model, use_periodicity);
256  cout << boost::format("# Group %d size is %d\n") % i % acceptor.size();
257  acceptors.push_back(acceptor);
258  }
259 
260  acceptor_names.push_back("Unbound/Other");
261 
262  uint n = traj_names.size() * donors.size();
263  uint m = acceptor_selections.size();
264 
265 
266 
267  Matrix M(m+1, n);
268 
269  if (verbose)
270  cerr << "Processing- ";
271 
272  for (uint k = 0; k<traj_names.size(); ++k) {
273  if (verbose)
274  cerr << traj_names[k] << " ";
275 
276  pTraj traj = createTrajectory(traj_names[k], model);
277  if (skip >= traj->nframes()) {
278  cerr << boost::format("Error- trajectory '%s' only has %d frames in it, but we are skipping %d frames...\n")
279  % traj_names[k]
280  % traj->nframes()
281  % skip;
282  exit(-20);
283  }
284 
285  BondMatrix B(m, donors.size());
286 
287  for (uint t = skip; t<traj->nframes(); ++t) {
288  traj->readFrame(t);
289  traj->updateGroupCoords(model);
290 
291  for (uint i=0; i<donors.size(); ++i) {
292  for (uint j=0; j<acceptors.size(); ++j) {
293  AtomicGroup found = donors[i].findHydrogenBonds(acceptors[j], true);
294  if (! found.empty())
295  B(j, i) += 1;
296  }
297  }
298  }
299 
300  for (uint i=0; i<donors.size(); ++i) {
301  double sum = 0.0;
302  for (uint j=0; j<m; ++j) {
303  double fraction = static_cast<double>(B(j, i)) / traj->nframes();
304  sum += fraction;
305  M(j, k * donors.size() + i) = fraction;
306  }
307  M(m, k * donors.size() + i) = (sum > 1.0) ? 0.0 : (1.0 - sum);
308  }
309 
310  }
311 
312  if (verbose)
313  cerr << endl;
314 
315  veDouble averages = rowAverage(M);
316  veDouble standards = rowStd(M, averages);
317 
318  for (uint i=0; i<averages.size(); ++i)
319  cout << boost::format("%-3d %-20s %.4f %.4f\n")
320  % i
321  % acceptor_names[i]
322  % averages[i]
323  % (standards[i] / (use_stderr ? sqrt(donors.size() * traj_names.size()) : 1.0));
324 
325 }
Simple matrix template class using policy classes to determine behavior.
Definition: MatrixImpl.hpp:53
Namespace for encapsulating options processing.
STL namespace.
bool parse(int argc, char *argv[])
Parses a command line, returning true if parsing was ok.
Options common to all tools (including –fullhelp)
AtomicGroup createSystem(const std::string &filename)
Factory function for reading in structure files.
Definition: sfactories.cpp:115
std::string invocationHeader(int argc, char *argv[])
Create an invocation header.
Definition: utils.cpp:124
std::vector< std::string > print() const
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
Namespace for most things not already encapsulated within a class.
std::string vectorAsStringWithCommas(const std::vector< std::string > &v)
Specialization for strings that sanitizes the contained strings.
Definition: utils.cpp:415
Combines a set of OptionsPackages.
AggregateOptions & add(OptionsPackage *pack)
Add a pointer to an OptionsPackage that will be used for options.