LOOS  v2.3.2
hmatrix.cpp
1 /*
2  hmatrix.cpp
3 
4  Writes out a matrix representing hbonds over time
5 */
6 
7 
8 /*
9 
10  This file is part of LOOS.
11 
12  LOOS (Lightweight Object-Oriented Structure library)
13  Copyright (c) 2010, Tod D. Romo
14  Department of Biochemistry and Biophysics
15  School of Medicine & Dentistry, University of Rochester
16 
17  This package (LOOS) is free software: you can redistribute it and/or modify
18  it under the terms of the GNU General Public License as published by
19  the Free Software Foundation under version 3 of the License.
20 
21  This package is distributed in the hope that it will be useful,
22  but WITHOUT ANY WARRANTY; without even the implied warranty of
23  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24  GNU General Public License for more details.
25 
26  You should have received a copy of the GNU General Public License
27  along with this program. If not, see <http://www.gnu.org/licenses/>.
28 */
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 // --------------- GLOBALS
46 
47 double length_low, length_high;
48 double max_angle;
49 bool use_periodicity;
50 string donor_selection, acceptor_selection;
51 string model_name;
52 string traj_name;
53 
54 uint currentTimeStep = 0;
55 
56 
57 
58 // ---------------
59 // @cond TOOLS_INTERNAL
60 
61 
62 string fullHelpMessage(void) {
63  string msg =
64  "\n"
65  "SYNOPSIS\n"
66  "\tHydrogen bond state for a trajectory as a matrix\n"
67  "\n"
68  "DESCRIPTION\n"
69  "\n"
70  "\tThis tool creates a matrix representing the state of putative hydrogen bonds\n"
71  "(1 = present, 0 = absent). Each row of the matrix is one frame (time-point) of\n"
72  "the trajectory. Each column corresponds to a possible h-bond acceptor. Only\n"
73  "one donor may be specified. Note that the donor is specified by selecting the\n"
74  "donated hydrogen. Criteria for putative hydrogen-bonds are an inner and outer\n"
75  "distance cutoff and an angle deviation from linear (in degrees) cutoff.\n"
76  "\n"
77  "EXAMPLES\n"
78  "\n"
79  "\thmatrix model.psf sim.dcd 'segid == \"PE1\" && resid == 4 && name == \"HE1\"'\\\n"
80  "\t 'name == \"O1\" && resname == \"PALM\"'\n"
81  "This example looks for hbonds between the HE1 hydrogen of residue 4 in the PE1 segment and\n"
82  "any palmitoyl carbonyl oxygen, O1.\n"
83  "\n"
84  "\thmatrix --blow 2.0 --bhi 4.0 --angle 25.0 model.psf sim.dcd \\\n"
85  "\t 'segid == \"PE1\" && resid == 4 && name == \"HE1\"'\\\n"
86  "\t 'name == \"O1\" && resname == \"PALM\"'\n"
87  "This example is the same as the above one, but with the hydrogen bond criteria changed\n"
88  "to greater than or equal to 2.0 angstroms and less than or equal to 4.0 angstroms, with\n"
89  "an angle of less than or equal to 25.0 degrees.\n"
90  "\n"
91  "SEE ALSO\n"
92  "\thbonds, hcorrelation\n";
93 
94  return(msg);
95 }
96 
97 
98 
99 
100 class ToolOptions : public opts::OptionsPackage {
101 public:
102  void addGeneric(po::options_description& o) {
103  o.add_options()
104  ("blow", po::value<double>(&length_low)->default_value(1.5), "Low cutoff for bond length")
105  ("bhi", po::value<double>(&length_high)->default_value(3.0), "High cutoff for bond length")
106  ("angle", po::value<double>(&max_angle)->default_value(30.0), "Max bond angle deviation from linear")
107  ("periodic", po::value<bool>(&use_periodicity)->default_value(false), "Use periodic boundary");
108  }
109 
110  void addHidden(po::options_description& o) {
111  o.add_options()
112  ("donor", po::value<string>(&donor_selection), "donor selection")
113  ("acceptor", po::value<string>(&acceptor_selection), "acceptor selection");
114 
115  }
116 
117  void addPositional(po::positional_options_description& opts) {
118  opts.add("donor", 1);
119  opts.add("acceptor", 1);
120  }
121 
122 
123  string help() const {
124  return("donor-selection acceptor-selection");
125  }
126 
127  string print() const {
128  ostringstream oss;
129  oss << boost::format("blow=%f,bhi=%f,angle=%f,periodic=%d,acceptor=\"%s\",donor=\"%s\"")
130  % length_low
131  % length_high
132  % max_angle
133  % use_periodicity
134  % acceptor_selection
135  % donor_selection;
136 
137  return(oss.str());
138  }
139 
140 };
141 
142 
143 
144 
145 // @endcond
146 
147 
148 
149 
150 
151 int main(int argc, char *argv[]) {
152  string hdr = invocationHeader(argc, argv);
153 
154 
155  opts::BasicOptions* bopts = new opts::BasicOptions(fullHelpMessage());
157  ToolOptions* topts = new ToolOptions;
158 
159  opts::AggregateOptions options;
160  options.add(bopts).add(tropts).add(topts);
161  if (! options.parse(argc, argv))
162  exit(-1);
163 
164 
165  AtomicGroup model = tropts->model;
166  pTraj traj = tropts->trajectory;
167 
168  if (use_periodicity && !traj->hasPeriodicBox()) {
169  cerr << "Error- trajectory has no periodic box information\n";
170  exit(-1);
171  }
172 
173  SimpleAtom::innerRadius(length_low);
174  SimpleAtom::outerRadius(length_high);
175  SimpleAtom::maxDeviation(max_angle);
176 
177  SAGroup donors = SimpleAtom::processSelection(donor_selection, model, use_periodicity);
178  if (donors.size() != 1) {
179  cerr << "Error- only specify one donor atom (the attached hydrogen)\n";
180  exit(-1);
181  }
182 
183  SAGroup acceptors = SimpleAtom::processSelection(acceptor_selection, model, use_periodicity);
184  BondMatrix bonds = donors[0].findHydrogenBondsMatrix(acceptors, traj, model);
185  writeAsciiMatrix(cout, bonds, hdr);
186 }
187 
pTraj trajectory
The trajectory, primed by the –skip value (if specified)
Simple matrix template class using policy classes to determine behavior.
Definition: MatrixImpl.hpp:53
Namespace for encapsulating options processing.
STL namespace.
AtomicGroup model
Model that describes the trajectory.
bool parse(int argc, char *argv[])
Parses a command line, returning true if parsing was ok.
Options common to all tools (including –fullhelp)
std::string invocationHeader(int argc, char *argv[])
Create an invocation header.
Definition: utils.cpp:124
std::ostream & writeAsciiMatrix(std::ostream &os, const Math::Matrix< T, P, S > &M, const std::string &meta, const Math::Range &start, const Math::Range &end, const bool trans=false, F fmt=F())
Write a submatrix to a stream.
Basic trajectory with a –skip option.
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
Namespace for most things not already encapsulated within a class.
Combines a set of OptionsPackages.
AggregateOptions & add(OptionsPackage *pack)
Add a pointer to an OptionsPackage that will be used for options.