LOOS  v2.3.2
alignment.cpp
1 #include <loos_defs.hpp>
2 
3 #include <AtomicGroup.hpp>
4 #include <Trajectory.hpp>
5 
6 #include <ensembles.hpp>
7 #include <alignment.hpp>
8 
9 #include <cmath>
10 
11 
12 namespace loos {
13 
14 
15 
16  namespace alignment {
17 
18 
19  // Core aligmnent routine. Assumes input coord vectors are already centered.
20  // Returns the SVD results as a tuple
21  SVDTupleVec kabschCore(const vecDouble& u, const vecDouble& v) {
22  int n = u.size() / 3;
23 
24  // Compute correlation matrix...
25  vecDouble R(9);
26 
27 #if defined(__linux__) || defined(__CYGWIN__) || defined(__FreeBSD__)
28  char ta = 'N';
29  char tb = 'T';
30  f77int three = 3;
31  double one = 1.0;
32  double zero = 0.0;
33 
34  dgemm_(&ta, &tb, &three, &three, &n, &one, u.data(), &three, v.data(), &three, &zero, R.data(), &three);
35 
36 #else
37 
38  cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans, 3, 3, n, 1.0, u.data(), 3, v.data(), 3, 0.0, R.data(), 3);
39 
40 #endif
41 
42  // Now compute the SVD of R...
43  char joba='G';
44  char jobu = 'U', jobv = 'V';
45  int mv = 0;
46  f77int m = 3, lda = 3, ldv = 3, lwork=100, info;
47  double work[lwork];
48  f77int nn = 3;
49  vecDouble S(3);
50  vecDouble V(9);
51 
52  dgesvj_(&joba, &jobu, &jobv, &m, &nn, R.data(), &lda, S.data(), &mv, V.data(), &ldv, work, &lwork, &info);
53 
54  if (info > 0) {
55  char op = 'E';
56  double eps = dlamch_(&op);
57  std::cerr << boost::format("Warning- SVD in kabschCore() failed to converge with info=%d and TOL=%e\n") %
58  info % (eps * sqrt(3.0));
59  std::cerr << "\tThis may happen periodically, causing the output to be 'mostly' orthogonal.\n";
60  std::cerr << "\tThe residual is typically small and should not appreciably affect the resulting\n";
61  std::cerr << "\tsuperposition. If this warning appears frequently, then please notify the LOOS\n";
62  std::cerr << "\tdevelopers at loos.maintainer@gmail.com\n";
63  // throw(NumericalError("SVD in alignment::kabschCore returned an error", info));
64  } else if (info < 0)
65  throw(NumericalError("SVD in alignment::kabschCore returned an error", info));
66 
67 
68  double dR = R[0]*R[4]*R[8] + R[3]*R[7]*R[2] + R[6]*R[1]*R[5] -
69  R[0]*R[7]*R[5] - R[3]*R[1]*R[8] - R[6]*R[4]*R[2];
70 
71 
72  double dV = V[0]*V[4]*V[8] + V[3]*V[7]*V[2] + V[6]*V[1]*V[5] -
73  V[0]*V[7]*V[5] - V[3]*V[1]*V[8] - V[6]*V[4]*V[2];
74 
75 
76  if (dR * dV < 0.0) {
77  S[2] = -S[2];
78 
79  R[6] = -R[6];
80  R[7] = -R[7];
81  R[8] = -R[8];
82  }
83 
84  return SVDTupleVec(R, S, V);
85  }
86 
87 
88  GCoord centerAtOrigin(vecDouble& v) {
89  GCoord c;
90 
91  for (uint i=0; i<v.size(); i += 3) {
92  c.x() += v[i];
93  c.y() += v[i+1];
94  c.z() += v[i+2];
95  }
96 
97  for (uint i=0; i<3; ++i)
98  c[i] = 3*c[i]/v.size();
99 
100  for (uint i=0; i<v.size(); i += 3) {
101  v[i] -= c.x();
102  v[i+1] -= c.y();
103  v[i+2] -= c.z();
104  }
105 
106  return c;
107  }
108 
109 
110 
111  // Return the RMSD only for a kabsch alignment between U and V assuming
112  // both are centered
113  double centeredRMSD(const vecDouble& U, const vecDouble& V) {
114 
115  int n = U.size();
116 
117  double ssu[3] = {0.0, 0.0, 0.0};
118  double ssv[3] = {0.0, 0.0, 0.0};
119 
120  for (int j=0; j<n; j += 3) {
121  for (uint i=0; i<3; ++i) {
122  ssu[i] += U[j+i] * U[j+i];
123  ssv[i] += V[j+i] * V[j+i];
124  }
125  }
126 
127  n /= 3;
128 
129  double E0 = ssu[0] + ssu[1] + ssu[2] + ssv[0] + ssv[1] + ssv[2];
130 
131  SVDTupleVec svd = kabschCore(U, V);
132 
133  vecDouble S(boost::get<1>(svd));
134  double ss = S[0] + S[1] + S[2];
135  double rmsd = std::sqrt(std::abs(E0-2.0*ss)/n);
136 
137  return(rmsd);
138  }
139 
140 
141 
142  // Return the RMSD only for a kabsch alignment between U and V
143  // Both will be centered first.
144  double alignedRMSD(const vecDouble& U, const vecDouble& V) {
145 
146  int n = U.size();
147 
148  vecDouble cU(U);
149  vecDouble cV(V);
150 
151  centerAtOrigin(cU);
152  centerAtOrigin(cV);
153 
154  SVDTupleVec svd = kabschCore(cU, cV);
155 
156  double ssu[3] = {0.0, 0.0, 0.0};
157  double ssv[3] = {0.0, 0.0, 0.0};
158 
159  for (int j=0; j<n; j += 3) {
160  for (uint i=0; i<3; ++i) {
161  ssu[i] += cU[j+i] * cU[j+i];
162  ssv[i] += cV[j+i] * cV[j+i];
163  }
164  }
165 
166  n /= 3;
167 
168  double E0 = ssu[0] + ssu[1] + ssu[2] + ssv[0] + ssv[1] + ssv[2];
169 
170  vecDouble S(boost::get<1>(svd));
171  double ss = S[0] + S[1] + S[2];
172  double rmsd = std::sqrt(std::abs(E0-2.0*ss)/n);
173  return(rmsd);
174  }
175 
176 
177 
178 
179  // Kabsch alignment between U and V, assuming both are centered.
180  // Returns the tranformation matrix to align U onto V.
181  GMatrix kabschCentered(const vecDouble& U, const vecDouble& V) {
182  SVDTupleVec svd = kabschCore(U, V);
183 
184  vecDouble R(boost::get<0>(svd));
185  vecDouble VV(boost::get<2>(svd));
186 
187  double M[9];
188 
189 #if defined(__linux__) || defined(__CYGWIN__) || defined(__FreeBSD__)
190  char ta = 'N';
191  char tb = 'T';
192  f77int three = 3;
193  double one = 1.0;
194  double zero = 0.0;
195 
196  dgemm_(&ta, &tb, &three, &three, &three, &one, R.data(), &three, VV.data(), &three, &zero, M, &three);
197 
198 #else
199 
200  cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans, 3, 3, 3, 1.0, R.data(), 3, VV.data(), 3, 0.0, M, 3);
201 
202 #endif
203 
204  // Construct the new transformation matrix... (W = M')
205  GMatrix Z;
206  for (uint i=0; i<3; i++)
207  for (uint j=0; j<3; j++)
208  Z(i,j) = M[i*3+j];
209 
210  return Z;
211 
212  }
213 
214 
215  GMatrix kabsch(const vecDouble& U, const vecDouble& V) {
216 
217  vecDouble cU(U);
218  vecDouble cV(V);
219 
220  GCoord U_center = centerAtOrigin(cU);
221  GCoord V_center = centerAtOrigin(cV);
222  GMatrix M = kabschCentered(cU, cV);
223 
224  XForm W;
225  W.identity();
226  W.translate(V_center);
227  W.concat(M);
228  W.translate(-U_center);
229 
230  return W.current();
231  }
232 
233 
234  void applyTransform(const GMatrix& M, vecDouble& v) {
235  for (uint i=0; i<v.size(); i += 3) {
236  GCoord c(v[i],v[i+1],v[i+2]);
237 
238  c = M * c;
239  v[i] = c.x();
240  v[i+1] = c.y();
241  v[i+2] = c.z();
242  }
243  }
244 
245 
246  vecDouble averageCoords(const vecMatrix& ensemble) {
247  uint m = ensemble.size();
248  uint n = ensemble[0].size();
249 
250  vecDouble avg(n, 0.0);
251  for (uint j=0; j<m; ++j)
252  for (uint i=0; i<n; ++i)
253  avg[i] += ensemble[j][i];
254 
255  for (uint i=0; i<n; ++i)
256  avg[i] /= m;
257 
258  return avg;
259  }
260 
261 
262  double rmsd(const vecDouble& u, const vecDouble& v) {
263  double rms = 0.0;
264  for (uint i=0; i<u.size(); i += 3) {
265  double l = 0.0;
266  for (uint j=0; j<3; ++j) {
267  double d = u[i+j] - v[i+j];
268  l += d*d;
269  }
270  rms += l;
271  }
272  rms = std::sqrt(3.0 * rms / u.size());
273 
274  return rms;
275  }
276 
277 
278 
279  }
280 
281 
282  // The following are the routines most should use... They behave the same way as the old
283  // ones from ensembles.cpp
284 
285 
286  boost::tuple<std::vector<XForm>,greal,int> iterativeAlignment(alignment::vecMatrix& ensemble,
287  greal threshold, int maxiter) {
288  using namespace alignment;
289 
290  int n = ensemble.size();
291  std::vector<XForm> xforms(n);
292 
293  // Start by aligning against the first structure in the ensemble
294 
295  vecDouble target(ensemble[0]);
296  centerAtOrigin(target);
297 
298  double rms;
299  int iter = 0;
300 
301  do {
302  for (int i = 0; i<n; i++) {
303  XForm M(kabsch(ensemble[i], target));
304  applyTransform(M.current(), ensemble[i]);
305  xforms[i].premult(M.current());
306  }
307 
308  vecDouble avg = averageCoords(ensemble);
309  rms = rmsd(target, avg);
310  target = avg;
311  ++iter;
312  } while (rms > threshold && iter <= maxiter );
313 
314  boost::tuple<std::vector<XForm>, greal, int> res(xforms, rms, iter);
315  return(res);
316  }
317 
318 
319  boost::tuple<std::vector<XForm>,greal,int> iterativeAlignment(std::vector<AtomicGroup>& ensemble,
320  greal threshold, int maxiter) {
321  using namespace alignment;
322 
323  int n = ensemble.size();
324  std::vector<XForm> xforms(n);
325 
326  // Start by aligning against the first structure in the ensemble
327  vecDouble target(ensemble[0].coordsAsVector());
328  centerAtOrigin(target);
329 
330  double rms;
331  int iter = 0;
332  do {
333  for (int i = 0; i<n; i++) {
334  XForm M(kabsch(ensemble[i].coordsAsVector(), target));
335  ensemble[i].applyTransform(M);
336  xforms[i].premult(M.current());
337  }
338 
339  AtomicGroup avg_structure = averageStructure(ensemble);
340  vecDouble avg = avg_structure.coordsAsVector();
341  rms = rmsd(target, avg);
342  target = avg;
343  ++iter;
344  } while (rms > threshold && iter <= maxiter );
345 
346  boost::tuple<std::vector<XForm>, greal, int> res(xforms, rms, iter);
347  return(res);
348  }
349 
350 
351 
352 
353  boost::tuple<std::vector<XForm>, greal, int> iterativeAlignment(const AtomicGroup& g,
354  pTraj& traj,
355  const std::vector<uint>& frame_indices,
356  greal threshold, int maxiter) {
357 
358  using namespace alignment;
359 
360  // Must first prime the loop...
361  AtomicGroup frame = g.copy();
362  traj->readFrame(frame_indices[0]);
363  traj->updateGroupCoords(frame);
364 
365  int nf = frame_indices.size();
366 
367  int iter = 0;
368  greal rms;
369  std::vector<XForm> xforms(nf);
370  AtomicGroup avg = frame.copy();
371 
372  AtomicGroup target = frame.copy();
373  target.centerAtOrigin();
374 
375  do {
376  // Compute avg internally so we don't have to read traj twice...
377  for (uint j=0; j<avg.size(); ++j)
378  avg[j]->coords() = GCoord(0,0,0);
379 
380  for (int i=0; i<nf; ++i) {
381 
382  traj->readFrame(frame_indices[i]);
383  traj->updateGroupCoords(frame);
384 
385  GMatrix M = frame.alignOnto(target);
386  xforms[i].load(M);
387 
388  for (uint j=0; j<avg.size(); ++j)
389  avg[j]->coords() += frame[j]->coords();
390  }
391 
392  for (uint j=0; j<avg.size(); ++j)
393  avg[j]->coords() /= nf;
394 
395  rms = avg.rmsd(target);
396  target = avg.copy();
397  ++iter;
398  } while (rms > threshold && iter <= maxiter);
399 
400  boost::tuple<std::vector<XForm>, greal, int> res(xforms, rms, iter);
401  return(res);
402  }
403 
404 
405  boost::tuple<std::vector<XForm>, greal, int> iterativeAlignment(const AtomicGroup& g,
406  pTraj& traj,
407  greal threshold, int maxiter) {
408 
409  std::vector<uint> framelist(traj->nframes());
410  for (uint i=0; i<traj->nframes(); ++i)
411  framelist[i] = i;
412 
413  return(iterativeAlignment(g, traj, framelist, threshold, maxiter));
414 
415 
416  }
417 
418 }
419 
boost::tuple< RealMatrix, RealMatrix, RealMatrix > svd(std::vector< AtomicGroup > &ensemble, bool align)
Compute the SVD of an ensemble with optional alignment (note RSVs returned are transposed) ...
Definition: ensembles.cpp:239
Basic 3-D coordinates class.
Definition: Coord.hpp:36
def iterativeAlignment(ensemble, threshold=1e-8, maxiter=1000)
Iteratively align a loos.pyloos.Trajectory object (or a list of AtomicGroups)
Definition: alignment.py:12
greal rmsd(const AtomicGroup &)
Compute the RMSD between two groups.
Specialized 4x4 Matrix class for handling coordinate transforms.
Definition: Coord.hpp:37
GMatrix alignOnto(const AtomicGroup &)
Superimposes the current group onto the passed group.
Definition: AG_linalg.cpp:196
AtomicGroup copy(void) const
Creates a deep copy of this group.
Definition: AtomicGroup.cpp:56
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
Namespace for most things not already encapsulated within a class.
AtomicGroup averageStructure(const std::vector< AtomicGroup > &ensemble)
Compute the average structure of a set of AtomicGroup objects.
Definition: ensembles.cpp:36
GCoord centerAtOrigin(void)
Translates the group so that the centroid is at the origin.