LOOS  v2.3.2
water-hist-lib.cpp
1 // -------------------------------------------------
2 // Water Histogram Library
3 // -------------------------------------------------
4 
5 /*
6  This file is part of LOOS.
7 
8  LOOS (Lightweight Object-Oriented Structure library)
9  Copyright (c) 2009 Tod D. Romo, Alan Grossfield
10  Department of Biochemistry and Biophysics
11  School of Medicine & Dentistry, University of Rochester
12 
13  This package (LOOS) is free software: you can redistribute it and/or modify
14  it under the terms of the GNU General Public License as published by
15  the Free Software Foundation under version 3 of the License.
16 
17  This package is distributed in the hope that it will be useful,
18  but WITHOUT ANY WARRANTY; without even the implied warranty of
19  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  GNU General Public License for more details.
21 
22  You should have received a copy of the GNU General Public License
23  along with this program. If not, see <http://www.gnu.org/licenses/>.
24 */
25 
26 
27 
28 #include <water-hist-lib.hpp>
29 
30 
31 namespace loos {
32  namespace DensityTools {
33 
34 
35  void ZClipEstimator::reinitialize(pTraj& traj, const std::vector<uint>& frames) {
36  std::vector<GCoord> bdd = getBounds(traj, water_, frames);
37 
38  bdd[0] -= 1;
39  if (bdd[0][2] > zclip_)
40  throw(std::runtime_error("Bulk zclip is too small"));
41  bdd[0][2] = zclip_;
42  bdd[1] += 1;
43 
44  GCoord gridsize = bdd[1] - bdd[0] + 1;
45  gridsize /= gridres_;
46  DensityGridpoint dims;
47  for (int i=0; i<3; ++i)
48  dims[i] = static_cast<int>(floor(gridsize[i] + 0.5));
49 
50  thegrid.resize(bdd[0], bdd[1], dims);
51  }
52 
53  void ZClipEstimator::operator()(const double density) {
54  for (AtomicGroup::const_iterator i = water_.begin(); i != water_.end(); ++i) {
55  GCoord c = (*i)->coords();
56  if (c.z() >= zclip_)
57  thegrid((*i)->coords()) += density;
58  }
59  }
60 
61  double ZClipEstimator::bulkDensity(void) const {
62  double mean = 0.0;
63  long n = 0;
64  for (long i = 0; i < thegrid.maxGridIndex(); ++i) {
65  if (!count_zero && thegrid(i) == 0.0)
66  continue;
67  mean += thegrid(i);
68  ++n;
69  }
70 
71  return(mean/n);
72  }
73 
74  double ZClipEstimator::stdDev(const double mean) const {
75  double std = 0.0;
76  long n = 0;
77  for (long i = 0; i < thegrid.maxGridIndex(); ++i)
78  if (thegrid(i) != 0.0) {
79  double d = thegrid(i) - mean;
80  std += d*d;
81  ++n;
82  }
83 
84  std = sqrt(std/(n-1));
85  return(std);
86  }
87 
88 
89  // ------------------------------------------------------------------------
90 
91  // Note: no checks on whether the z-slice is sensible...
92 
93  void ZSliceEstimator::reinitialize(pTraj& traj, const std::vector<uint>& frames) {
94  std::vector<GCoord> bdd = getBounds(traj, water_, frames);
95 
96  bdd[0] -= 1;
97  bdd[0][2] = zmin_;
98  bdd[1] += 1;
99  bdd[1][2] = zmax_;
100 
101  GCoord gridsize = bdd[1] - bdd[0] + 1;
102  gridsize /= gridres_;
103  DensityGridpoint dims;
104  for (int i=0; i<3; ++i)
105  dims[i] = static_cast<int>(floor(gridsize[i] + 0.5));
106 
107  thegrid.resize(bdd[0], bdd[1], dims);
108  }
109 
110  void ZSliceEstimator::operator()(const double density) {
111  for (AtomicGroup::const_iterator i = water_.begin(); i != water_.end(); ++i) {
112  GCoord c = (*i)->coords();
113  if (c.z() >= zmin_ && c.z() < zmax_)
114  thegrid((*i)->coords()) += density;
115  }
116  }
117 
118  double ZSliceEstimator::bulkDensity(void) const {
119  double mean = 0.0;
120  long n = 0;
121  for (long i = 0; i < thegrid.maxGridIndex(); ++i) {
122  if (!count_zero && thegrid(i) == 0.0)
123  continue;
124  mean += thegrid(i);
125  ++n;
126  }
127 
128  return(mean/n);
129  }
130 
131  double ZSliceEstimator::stdDev(const double mean) const {
132  double std = 0.0;
133  long n = 0;
134  for (long i = 0; i < thegrid.maxGridIndex(); ++i)
135  if (thegrid(i) != 0.0) {
136  double d = thegrid(i) - mean;
137  std += d*d;
138  ++n;
139  }
140 
141  std = sqrt(std/(n-1));
142  return(std);
143  }
144 
145 
146  // ------------------------------------------------------------------------
147 
148 
149  void WaterHistogrammer::setGrid(const GCoord& min, const GCoord& max, const double resolution) {
150  GCoord gridsize = max - min + 1;
151  gridsize /= resolution;
152  DensityGridpoint dims;
153  for (int i=0; i<3; ++i)
154  dims[i] = static_cast<int>(floor(gridsize[i] + 0.5));
155 
156  grid_.resize(min, max, dims);
157  }
158 
159 
160  void WaterHistogrammer::setGrid(pTraj& traj, const std::vector<uint>& frames, const double resolution, const double pad) {
161 
162  std::vector<GCoord> bdd(2);
163  for (uint i=0; i<frames.size(); ++i) {
164  traj->readFrame(i);
165  traj->updateGroupCoords(protein_);
166  std::vector<GCoord> fbdd = the_filter->boundingBox(protein_);
167  for (uint j=0; j<3; ++j) {
168  if (fbdd[0][j] < bdd[0][j])
169  bdd[0][j] = fbdd[0][j];
170  if (fbdd[1][j] > bdd[1][j])
171  bdd[1][j] = fbdd[1][j];
172  }
173  }
174 
175 
176  setGrid(bdd[0] - pad, bdd[1] + pad, resolution);
177  }
178 
179 
180  void WaterHistogrammer::accumulate(const double density) {
181  std::vector<int> picks = the_filter->filter(water_, protein_);
182  for (uint i = 0; i<picks.size(); ++i)
183  if (picks[i]) {
184  GCoord c = water_[i]->coords();
185 
186  if (!grid_.inRange(grid_.gridpoint(c)))
187  ++out_of_bounds;
188  else
189  grid_(c) += density;
190  }
191  (*estimator_)(density);
192  }
193 
194  void WaterHistogrammer::accumulate(pTraj& traj, const std::vector<uint>& frames) {
195  estimator_->reinitialize(traj, frames);
196  double density = 1.0 / frames.size();
197  for (std::vector<uint>::const_iterator i = frames.begin(); i != frames.end(); ++i) {
198  traj->readFrame(*i);
199  traj->updateGroupCoords(protein_);
200  traj->updateGroupCoords(water_);
201 
202  accumulate(density);
203  }
204  }
205 
206  };
207 };
DensityGridpoint gridpoint(const loos::GCoord &x) const
Converts a real-space coordinate into grid coords.
bool inRange(const DensityGridpoint &g) const
Checks to make sure the gridpoint lies within the grid boundaries.
std::vector< GCoord > getBounds(pTraj &traj, AtomicGroup &g, const std::vector< uint > &indices)
Get the max bounding box for a group over the trajectory.
Definition: water-lib.cpp:35
STL namespace.
virtual std::vector< int > filter(const loos::AtomicGroup &, const loos::AtomicGroup &)=0
Given a molecule and a set of waters, pick which waters are inside.
virtual std::vector< loos::GCoord > boundingBox(const loos::AtomicGroup &)=0
Calculate the appropriate bounding box (given the molecule)
Namespace for most things not already encapsulated within a class.