LOOS  v2.3.2
water-sides.cpp
1 /*
2  water-sides.cpp
3 
4  Desc:
5  Classifies a water as being on one side of the membrane or the
6  other or inside the membrane (1 = upper, 0 = inside, -1 = lower).
7 */
8 
9 
10 /*
11  This file is part of LOOS.
12 
13  LOOS (Lightweight Object-Oriented Structure library)
14  Copyright (c) 2008, Tod D. Romo, Alan Grossfield
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 <limits>
34 #include <boost/format.hpp>
35 #include <DensityTools.hpp>
36 
37 
38 using namespace std;
39 using namespace loos;
40 
41 
42 
43 // @cond TOOLS_INTERNAL
44 
45 typedef std::pair<double,double> Range;
47 
48 Range membrane(0.0, 0.0);
49 string model_name, traj_name, selection_string;
50 
51 enum Location { UPPER = 1, MEMBRANE = 0, LOWER = -1 };
52 
53 
54 string fullHelpMessage(void) {
55  string msg =
56  "\n"
57  "SYNOPSIS\n"
58  "\n"
59  "\tClassify waters as above, below, or inside a membrane (based on Z-coordinate)\n"
60  "\n"
61  "DESCRIPTION\n"
62  "\n"
63  "\twater-sides constructs a TxN matrix where T is the size of the trajectory (# of frames)\n"
64  "and N is the number of water molecules. Each element has a value of 1 (above membrane),\n"
65  "0 (inside membrane), or -1 (below membrane). The classification of the water is based\n"
66  "solely on it's z-coordinate and the range specified on the command line.\n"
67  "\n"
68  "\nEXAMPLES\n"
69  "\twater-sides foo.pdb foo.dcd -15 15\n"
70  "This example uses the default water selection (\"name == 'OH2'\") and places the\n"
71  "membrane at -15 <= Z <= 15\n"
72  "\n"
73  "\twater-sides --selection 'name == \"HOH\"' foo.pdb foo.dcd -25 20\n"
74  "This example picks all atoms called \"HOH\" as waters and places the membrane\n"
75  "at -25 <= Z <= 20\n"
76  ;
77 
78  return(msg);
79 }
80 
81 
82 class WaterSidesOptions : public opts::OptionsPackage {
83 public:
84 
85  WaterSidesOptions() : lower_bounds(0.0), upper_bounds(0.0) { }
86 
87  void addHidden(po::options_description& options) {
88  options.add_options()
89  ("lower", po::value<double>(&lower_bounds), "Lower leaflet bounds")
90  ("upper", po::value<double>(&upper_bounds), "Upper leaflet bounds");
91  }
92 
93  void addPositional(po::positional_options_description& options) {
94  options.add("lower", 1);
95  options.add("upper", 1);
96  }
97 
98  bool check(po::variables_map& map) {
99  return(! (map.count("lower") && map.count("upper")) );
100  }
101 
102  string help() const { return("membrane-lower-bounds membrane-upper-bounds"); }
103  string print() const {
104  ostringstream oss;
105  oss << boost::format("lower=%f, upper=%f") % lower_bounds % upper_bounds;
106  return(oss.str());
107  }
108 
109  double lower_bounds, upper_bounds;
110 };
111 
112 // @endcond
113 
114 
115 Range parseRange(const string& s) {
116  double a, b;
117  int i = sscanf(s.c_str(), "%lf:%lf", &a, &b);
118  if (i != 2) {
119  cerr << "Parse error with " << s << endl;
120  exit(-1);
121  }
122  return(Range(a,b));
123 }
124 
125 
126 
127 
128 int main(int argc, char *argv[]) {
129  string hdr = invocationHeader(argc, argv);
130  opts::BasicOptions *basic_opts = new opts::BasicOptions(fullHelpMessage());
131  opts::BasicSelection *basic_selection = new opts::BasicSelection("name == 'OH2'");
133  WaterSidesOptions *my_opts = new WaterSidesOptions;
134 
135  opts::AggregateOptions options;
136  options.add(basic_opts).add(basic_selection).add(basic_traj).add(my_opts);
137  if (!options.parse(argc, argv)) {
138  exit(0);
139  }
140 
141  AtomicGroup model = basic_traj->model;
142  pTraj traj = basic_traj->trajectory;
143  AtomicGroup subset = selectAtoms(model, basic_selection->selection);
144  vector<uint> frames = basic_traj->frameList();
145 
146  uint n = subset.size();
147  uint m = frames.size();
148 
149  Matrix M(m,n+1);
150  for (uint j=0; j<m; ++j) {
151  traj->readFrame(frames[j]);
152  traj->updateGroupCoords(subset);
153  M(j, 0) = frames[j];
154  for (uint i=0; i<n; ++i) {
155  GCoord c = subset[i]->coords();
156  Location l;
157  if (c[2] > membrane.second)
158  l = UPPER;
159  else if (c[2] >= membrane.first)
160  l = MEMBRANE;
161  else
162  l = LOWER;
163  M(j, i+1) = l;
164  }
165  }
166 
167  writeAsciiMatrix(cout, M, hdr);
168 }
Trajectory with either a –range or –skip.
Simple matrix template class using policy classes to determine behavior.
Definition: MatrixImpl.hpp:53
Provides a single LOOS selection (–selection)
STL namespace.
bool parse(int argc, char *argv[])
Parses a command line, returning true if parsing was ok.
AtomicGroup model
Model that describes the trajectory.
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.
AtomicGroup selectAtoms(const AtomicGroup &source, const std::string selection)
Applies a string-based selection to an atomic group...
Definition: utils.cpp:195
std::pair< uint, uint > Range
Specify a range for columns/rows [first,second)
Definition: MatrixImpl.hpp:49
std::vector< uint > frameList() const
Returns the list of frames the user requested.
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.