LOOS  v2.3.2
AG_numerical.cpp
1 /*
2  AG_numerical.cpp
3 
4  Numerical methods for AtomicGroup class
5 */
6 
7 /*
8  This file is part of LOOS.
9 
10  LOOS (Lightweight Object-Oriented Structure library)
11  Copyright (c) 2008, Tod D. Romo, Alan Grossfield
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 #include <ios>
31 #include <sstream>
32 #include <iomanip>
33 
34 #include <assert.h>
35 #include <vector>
36 
37 #include <algorithm>
38 
39 #include <boost/random.hpp>
40 
41 #include <AtomicGroup.hpp>
42 #include <utils.hpp>
43 
44 
45 namespace loos {
46 
47  // Bounding box for all atoms in this group
48  // Returns a vector containing 2 GCoords, one containing
49  // (minx,miny,minz) and the other (maxx,maxy,maxz)
50  std::vector<GCoord> AtomicGroup::boundingBox(void) const {
51  greal min[3] = {0,0,0}, max[3] = {0,0,0};
52  const_iterator i;
53  int j;
54  std::vector<GCoord> res(2);
55  GCoord c;
56 
57  if (atoms.size() == 0) {
58  res[0] = c;
59  res[1] = c;
60  return(res);
61  }
62 
63  for (j=0; j<3; j++)
64  min[j] = max[j] = (atoms[0]->coords())[j];
65 
66  for (i=atoms.begin()+1; i != atoms.end(); i++)
67  for (j=0; j<3; j++) {
68  if (max[j] < ((*i)->coords())[j])
69  max[j] = ((*i)->coords())[j];
70  if (min[j] > ((*i)->coords())[j])
71  min[j] = ((*i)->coords())[j];
72  }
73 
74  c.set(min[0], min[1], min[2]);
75  res[0] = c;
76  c.set(max[0], max[1], max[2]);
77  res[1] = c;
78 
79  return(res);
80  }
81 
82  // Geometric center of the group
84  GCoord c(0,0,0);
85  const_iterator i;
86 
87  // Optimization for groups containing only one atom (such as a
88  // water heavy-atom)
89  if (atoms.size() == 1)
90  return(atoms[0]->coords());
91 
92  for (i = atoms.begin(); i != atoms.end(); i++)
93  c += (*i)->coords();
94 
95  c /= atoms.size();
96  return(c);
97  }
98 
99 
101  GCoord c(0,0,0);
102  const_iterator i;
103 
104  // Optimization for groups containing only one atom (such as a
105  // water heavy-atom)
106  if (atoms.size() == 1) {
107  return(atoms[0]->coords());
108  }
109 
110  for (i=atoms.begin(); i != atoms.end(); i++) {
111  c += (*i)->mass() * (*i)->coords();
112  }
113  c /= totalMass();
114  return(c);
115  }
116 
117  //* Note: this code implicitly assumes the group is neutral,
118  // since it sums based on atomic number
120  GCoord c(0,0,0);
121  const_iterator i;
122  int electrons = 0;
123 
124  for (i=atoms.begin(); i != atoms.end(); i++) {
125  int an = (*i)->atomic_number();
126  c += an * (*i)->coords();
127  electrons += an;
128  }
129  c /= electrons;
130  return(c);
131  }
132 
134  GCoord center = centroid();
135  GCoord moment(0,0,0);
136  const_iterator i;
137  for (i=atoms.begin(); i != atoms.end(); i++) {
138  moment += (*i)->charge() * ((*i)->coords() - center);
139  }
140  return(moment);
141  }
142 
143  greal AtomicGroup::totalCharge(void) const {
144  const_iterator i;
145  greal charge = 0.0;
146 
147  for (i = atoms.begin(); i != atoms.end(); i++)
148  charge += (*i)->charge();
149 
150  return(charge);
151  }
152 
153  greal AtomicGroup::totalMass(void) const {
154  const_iterator i;
155  greal mass = 0.0;
156 
157  for (i = atoms.begin(); i != atoms.end(); i++)
158  mass += (*i)->mass();
159 
160  return(mass);
161  }
162 
163  // Geometric max radius of the group (relative to the centroid)
164  greal AtomicGroup::radius(const bool use_atom_as_reference) const {
165  GCoord c;
166  if (use_atom_as_reference) {
167  c = atoms[0]->coords();
168  }
169  else {
170  c = centroid();
171  }
172  greal radius = 0.0;
173  const_iterator i;
174 
175  for (i=atoms.begin(); i != atoms.end(); i++) {
176  greal d = c.distance2((*i)->coords());
177  if (d > radius)
178  radius = d;
179  }
180 
181  radius = sqrt(radius);
182  return(radius);
183  }
184 
185 
186 
187  greal AtomicGroup::radiusOfGyration(void) const {
188  GCoord c = centerOfMass();
189  greal radius = 0;
190  const_iterator i;
191 
192  for (i = atoms.begin(); i != atoms.end(); i++)
193  radius += c.distance2((*i)->coords());
194 
195  radius = sqrt(radius / atoms.size());
196  return(radius);
197  }
198 
204  greal AtomicGroup::sphericalVariance(const pAtom target) const {
205  return sphericalVariance(target->coords());
206  }
207 
208 
209  greal AtomicGroup::sphericalVariance(const GCoord target) const {
210  GCoord var;
211  for (const_iterator i = atoms.begin(); i != atoms.end(); i++) {
212  GCoord vec = (*i)->coords() - target;
213  greal length = vec.length();
214  var += vec / length;
215  }
216 
217  greal variance = var.length() / atoms.size();
218  return variance;
219 
220  }
221 
222  greal AtomicGroup::rmsd(const AtomicGroup& v) {
223 
224  if (size() != v.size())
225  throw(LOOSError("Cannot compute RMSD between groups with different sizes"));
226 
227 
228  int n = size();
229  double d = 0.0;
230  for (int i = 0; i < n; i++) {
231  GCoord x = atoms[i]->coords();
232  GCoord y = v.atoms[i]->coords();
233  d += x.distance2(y);
234  }
235 
236  d = sqrt(d/n);
237 
238  return(d);
239  }
240 
241 
242 
243  std::vector<GCoord> AtomicGroup::getTransformedCoords(const XForm& M) const {
244  std::vector<GCoord> crds(atoms.size());
245  const_iterator i;
246  GMatrix W = M.current();
247  int j = 0;
248 
249  for (i = atoms.begin(); i != atoms.end(); i++) {
250  GCoord res = W * (*i)->coords();
251  crds[j++] = res;
252  }
253 
254  return(crds);
255  }
256 
257 
258  void AtomicGroup::translate(const GCoord & v) {
259  iterator i;
260  for (i = atoms.begin(); i != atoms.end(); i++)
261  (*i)->coords() += v;
262  }
263 
265  iterator i;
266  GMatrix W = M.current();
267 
268  for (i = atoms.begin(); i != atoms.end(); i++)
269  (*i)->coords() = W * (*i)->coords();
270 
271  }
272 
273 
274  void AtomicGroup::rotate(const GCoord& axis, const greal angle_in_degrees) {
275  XForm M;
276 
277  GCoord center = centroid();
278  M.translate(center);
279  M.rotate(axis, -angle_in_degrees);
280  M.translate(-center);
281  applyTransform(M);
282  }
283 
284 
285  std::vector<double> AtomicGroup::coordsAsVector() const {
286  std::vector<double> v(size() * 3);
287 
288  uint k = 0;
289  for (uint i=0; i<size(); ++i) {
290  v[k++] = atoms[i]->coords().x();
291  v[k++] = atoms[i]->coords().y();
292  v[k++] = atoms[i]->coords().z();
293  }
294 
295  return(v);
296  }
297 
298 
299  // Returns a newly allocated array of double coords in row-major
300  // order...
301  double* AtomicGroup::coordsAsArray(void) const {
302  double *A;
303  int n = size();
304 
305  A = new double[n*3];
306  int k = 0;
307  int i;
308  for (i=0; i<n; i++) {
309  A[k++] = atoms[i]->coords().x();
310  A[k++] = atoms[i]->coords().y();
311  A[k++] = atoms[i]->coords().z();
312  }
313 
314  return(A);
315  }
316 
317  // Returns a newly allocated array of double coords in row-major order
318  // transformed by the current transformation.
319  double* AtomicGroup::transformedCoordsAsArray(const XForm& M) const {
320  double *A;
321  GCoord x;
322  int n = size();
323  GMatrix W = M.current();
324 
325  A = new double[n*3];
326  int k = 0;
327  int i;
328  for (i=0; i<n; i++) {
329  x = W * atoms[i]->coords();
330  A[k++] = x.x();
331  A[k++] = x.y();
332  A[k++] = x.z();
333  }
334 
335  return(A);
336  }
337 
338 
339 
341  GCoord c = centroid();
342  iterator i;
343 
344  for (i = atoms.begin(); i != atoms.end(); i++)
345  (*i)->coords() -= c;
346 
347  return(c);
348  }
349 
350 
351 
352  void AtomicGroup::perturbCoords(const greal rms) {
353  int i, n = size();
354  GCoord r;
355 
356  for (i=0; i<n; i++) {
357  r.random();
358  r *= rms;
359  atoms[i]->coords() += r;
360  }
361  }
362 
363 
364 }
void rotate(const GCoord &, const greal)
Definition: XForm.cpp:64
void set(const T x, const T y, const T z)
Short-cut to set the cartesian coordinates...
Definition: Coord.hpp:144
void random(void)
Generate a random vector on a unit sphere.
Definition: Coord.hpp:461
std::vector< GCoord > boundingBox(void) const
Bounding box for the group...
double length(void) const
Length of the coordinate (as a vector)
Definition: Coord.hpp:423
Matrix class for handling coordinate transforms...
Definition: XForm.hpp:91
double distance2(const Coord< T > &o) const
Distance squared between two coordinates.
Definition: Coord.hpp:429
greal rmsd(const AtomicGroup &)
Compute the RMSD between two groups.
void rotate(const GCoord &axis, const greal angle_in_degrees)
Rotate group's coordinates (right-handed, about centroid)
GCoord centerOfMass(void) const
Center of mass of the group (in group coordinates)
Specialized 4x4 Matrix class for handling coordinate transforms.
Definition: Coord.hpp:37
void translate(const greal, const greal, const greal)
Translation matrix.
Definition: XForm.cpp:37
greal radius(const bool use_atom_as_reference=false) const
Maximum radius from centroid of all atoms (not gyration)
void translate(const GCoord &v)
Translate an atomic group by vector v.
std::vector< GCoord > getTransformedCoords(const XForm &) const
Returns a vector of coordinates transformed by the passed XForm.
GMatrix current(void) const
Get the current trasnformation.
Definition: XForm.cpp:110
GCoord centerOfElectrons(void) const
Analogous to center of mass.
Generic LOOS exception.
Definition: exceptions.hpp:40
GCoord dipoleMoment(void) const
Dipole moment, relative to group's centroid.
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
void applyTransform(const XForm &)
Apply the given transform to the group's coordinates...
Namespace for most things not already encapsulated within a class.
void perturbCoords(const greal)
Each atom is moved in a random direction by a vector of the passed size.
GCoord centerAtOrigin(void)
Translates the group so that the centroid is at the origin.
GCoord centroid(void) const
Centroid of atoms (ignores mass, operates in group coordinates)
greal sphericalVariance(const pAtom) const
Spherical variance of group with respect to target atom.