LOOS  v2.3.2
AG_linalg.cpp
1 /*
2  AG_linalg.cpp
3 
4  Linear-algebra routines (requiring ATLAS or vecLib) for AtomicGroup
5  class
6 */
7 
8 /*
9  This file is part of LOOS.
10 
11  LOOS (Lightweight Object-Oriented Structure library)
12  Copyright (c) 2008, Tod D. Romo, Alan Grossfield
13  Department of Biochemistry and Biophysics
14  School of Medicine & Dentistry, University of Rochester
15 
16  This package (LOOS) is free software: you can redistribute it and/or modify
17  it under the terms of the GNU General Public License as published by
18  the Free Software Foundation under version 3 of the License.
19 
20  This package is distributed in the hope that it will be useful,
21  but WITHOUT ANY WARRANTY; without even the implied warranty of
22  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23  GNU General Public License for more details.
24 
25  You should have received a copy of the GNU General Public License
26  along with this program. If not, see <http://www.gnu.org/licenses/>.
27 */
28 
29 
30 
31 #include <ios>
32 #include <sstream>
33 #include <iomanip>
34 
35 #include <assert.h>
36 #include <algorithm>
37 
38 #include <boost/random.hpp>
39 
40 #include <AtomicGroup.hpp>
41 #include <alignment.hpp>
42 
43 
44 
45 namespace loos {
46 
47 
48  std::vector<GCoord> AtomicGroup::momentsOfInertia(void) const {
49  Math::Matrix<double, Math::ColMajor> I(3, 3); // This gets initialized to zero...
50  GCoord c = centerOfMass();
51 
52  for (uint i = 0; i < atoms.size(); ++i) {
53 
54  GCoord u = atoms[i]->coords() - c;
55  double m = atoms[i]->mass();
56  I(0,0) += m * (u.y() * u.y() + u.z() * u.z());
57  I(1,0) += m * u.x() * u.y();
58  I(2,0) += m * u.x() * u.z();
59  I(1,1) += m * (u.x() * u.x() + u.z() * u.z());
60  I(2,1) += m * u.y() * u.z();
61  I(2,2) += m * (u.x() * u.x() + u.y() * u.y());
62  }
63 
64  I(1,0) = I(0,1) = -I(1,0);
65  I(2,0) = I(0,2) = -I(2,0);
66  I(2,1) = I(1,2) = -I(2,1);
67 
68 
69  // Now compute the eigen-decomp...
70  char jobz = 'V', uplo = 'U';
71  f77int nn;
72  f77int lda = 3;
73  double W[3], work[128];
74  f77int lwork = 128; // ??? Just a guess for sufficient storage
75  f77int info;
76  nn = 3;
77 
78  dsyev_(&jobz, &uplo, &nn, I.get(), &lda, W, work, &lwork, &info);
79  if (info < 0)
80  throw(NumericalError("dsyev_ reported an argument error.", info));
81 
82  if (info > 0)
83  throw(NumericalError("dsyev_ failed to converge.", info));
84 
85  std::vector<GCoord> results(4);
86 
87  for (int i=0; i<3; i++)
88  results[2-i] = GCoord( I(0,i), I(1,i), I(2,i) );
89 
90 
91  // Now push the eigenvalues on as a GCoord...
92  c[0] = W[2];
93  c[1] = W[1];
94  c[2] = W[0];
95 
96  c /= size();
97 
98  results[3] = c;
99 
100  return(results);
101  }
102 
103 
104  std::vector<GCoord> AtomicGroup::principalAxes(void) const {
105  // Extract out the group's coordinates...
106  int i;
107  int n = size();
108  double M[3] = {0.0, 0.0, 0.0};
109  int k = 0;
110 
111  double *A = coordsAsArray();
112  for (i=0; i<n; i++) {
113  M[0] += A[k++];
114  M[1] += A[k++];
115  M[2] += A[k++];
116  }
117 
118  M[0] /= n;
119  M[1] /= n;
120  M[2] /= n;
121 
122  // Subtract off the mean...
123  for (i=k=0; i<n; i++) {
124  A[k++] -= M[0];
125  A[k++] -= M[1];
126  A[k++] -= M[2];
127  }
128 
129  // Multiply A*A'...
130  double C[9];
131 #if defined(__linux__) || defined(__CYGWIN__) || defined(__FreeBSD__)
132  char ta = 'N';
133  char tb = 'T';
134  f77int three = 3;
135  double zero = 0.0;
136  double one = 1.0;
137 
138  dgemm_(&ta, &tb, &three, &three, &n, &one, A, &three, A, &three, &zero, C, &three);
139 
140 #else
141  cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans,
142  3, 3, n, 1.0, A, 3, A, 3, 0.0, C, 3);
143 #endif
144 
145  delete[] A;
146 
147  // Now compute the eigen-decomp...
148  char jobz = 'V', uplo = 'U';
149  f77int nn;
150  f77int lda = 3;
151  double W[3], work[128];
152  f77int lwork = 128; // ??? Just a guess for sufficient storage to be
153  // efficient...
154  f77int info;
155  nn = 3;
156 
157  dsyev_(&jobz, &uplo, &nn, C, &lda, W, work, &lwork, &info);
158  if (info < 0)
159  throw(NumericalError("dsyev_ reported an argument error.", info));
160 
161  if (info > 0)
162  throw(NumericalError("dsyev_ failed to converge.", info));
163 
164  std::vector<GCoord> results(4);
165  GCoord c;
166 
167  k = 0;
168  for (i=0; i<3; i++) {
169  c[0] = C[k++];
170  c[1] = C[k++];
171  c[2] = C[k++];
172  results[2-i] = c;
173  }
174 
175  // Now push the eigenvalues on as a GCoord...
176  c[0] = W[2];
177  c[1] = W[1];
178  c[2] = W[0];
179  c /= size(); // Scale by # of atoms...
180 
181  results[3] = c;
182 
183  return(results);
184  }
185 
186 
187 
189  alignment::vecDouble u = coordsAsVector();
190  alignment::vecDouble v = grp.coordsAsVector();
191 
192  return(alignment::kabsch(u, v));
193  }
194 
195 
197  XForm W;
198  GMatrix M = superposition(grp);
199 
200  W.load(M);
201  applyTransform(W);
202 
203  return(M);
204  }
205 
206 }
Simple matrix template class using policy classes to determine behavior.
Definition: MatrixImpl.hpp:53
GMatrix superposition(const AtomicGroup &)
Calculates the transformation matrix for superposition of groups.
Definition: AG_linalg.cpp:188
Matrix class for handling coordinate transforms...
Definition: XForm.hpp:91
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
std::vector< GCoord > momentsOfInertia(void) const
Computes the moments of inertia for a group.
Definition: AG_linalg.cpp:48
GMatrix alignOnto(const AtomicGroup &)
Superimposes the current group onto the passed group.
Definition: AG_linalg.cpp:196
Exception caused by a blas/atlas error.
Definition: exceptions.hpp:91
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 load(const GMatrix &)
Load a matrix onto the current transform.
Definition: XForm.cpp:28
std::vector< GCoord > principalAxes(void) const
Compute the principal axes of a group.
Definition: AG_linalg.cpp:104