Rosetta
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
bridgeObjects.cc
Go to the documentation of this file.
1 // -*- mode:c++;tab-width:2;indent-tabs-mode:t;show-trailing-whitespace:t;rm-trailing-spaces:t -*-
2 // vi: set ts=2 noet:
3 //
4 // (c) Copyright Rosetta Commons Member Institutions.
5 // (c) This file is part of the Rosetta software suite and is made available under license.
6 // (c) The Rosetta software is developed by the contributing members of the Rosetta Commons.
7 // (c) For more information, see http://www.rosettacommons.org. Questions about this can be
8 // (c) addressed to University of Washington UW TechTransfer, email: license@u.washington.edu.
9 
10 /// @file bridgeObjects.cc
11 /// @brief Solves the triaxial loop closure problem for a system of atoms
12 /// @author Evangelos A. Coutsias
13 /// @author Daniel J. Mandell
14 
15 // C++ headers
16 #include <cmath>
17 // Utility headers
18 #include <utility/vector1.hh>
19 // Rosetta Headers
20 #include <numeric/types.hh>
21 //#include <basic/Tracer.hh> // tracer output
22 #include <numeric/kinematic_closure/dixon.hh>
23 #include <numeric/kinematic_closure/kinematic_closure_helpers.hh>
24 
25 #include <iostream>
26 
27 // Constants **DJM: replace with standard from mini or rosetta++**
28 #define rad2deg 57.2957795130823208767981548141051703
29 #define deg2rad 0.0174532925199432957692369076848861271
30 #define pi 3.14159265358979323846264338327950288
31 #define pidegs 180.0
32 
33 using numeric::Real;
34 //static basic::Tracer std::cout( "numeric.kinematic_closure.bridgeObjects" );
35 
36 namespace numeric {
37 namespace kinematic_closure {
38 
39 void triangle (const utility::vector1<Real>& vbond, utility::vector1<Real>& calpha, utility::vector1<Real>& salpha) {
40  Real d11, d12, d13, d22, d23, d33;
41  d11 = vbond[1]*vbond[1];
42  d22 = vbond[2]*vbond[2];
43  d33 = vbond[3]*vbond[3];
44  d12 = 2*vbond[1]*vbond[2];
45  d13 = 2*vbond[1]*vbond[3];
46  d23 = 2*vbond[2]*vbond[3];
47  calpha[1]=(d22-d33-d11)/d13; // exterior angle
48  salpha[1]=std::sqrt(1-calpha[1]*calpha[1]);
49  calpha[2]=(d33-d11-d22)/d12; // exterior angle
50  salpha[2]=std::sqrt(1-calpha[2]*calpha[2]);
51  calpha[3]=(d11-d22-d33)/d23; // exterior angle
52  salpha[3]=std::sqrt(1-calpha[3]*calpha[3]);
53 }
54 
55 void sincos (const utility::vector1<Real>& theta,
56  const int flag,
57  utility::vector1<Real>& cosine,
58  utility::vector1<Real>& sine)
59 {
60  utility::vector1<Real> a, aa, a1;
61  if ( flag == 1 ) {
62  for ( int i=1; i<=3; i++ ) {
63  cosine[i]=std::cos(theta[i]);
64  sine[i]=std::sin(theta[i]);
65  }
66  } else {
67  a.resize(3);
68  aa.resize(3);
69  a1.resize(3);
70  for ( int i=1; i<=3; i++ ) {
71  a[i]=std::tan(theta[i]/2);
72  aa[i]=a[i]*a[i];
73  a1[i]=1+aa[i];
74  sine[i]=2*a[i]/a1[i];
75  cosine[i]=(1-aa[i])/a1[i];
76  }
77  }
78 }
79 
80 /*
81 * triaxialCoefficients
82 * --------------------
83 * Sets up characteristic polynomial matrices A,B,C,D from tripeptide parameters
84 */
85 
86 void triaxialCoefficients(
87  const utility::vector1<Real>& vb,
88  const utility::vector1<Real>& xi,
89  const utility::vector1<Real>& eta,
90  const utility::vector1<Real>& delta,
91  const utility::vector1<Real>& theta,
92  const utility::vector1<int>& order,
93  utility::vector1<utility::vector1<Real> >& A,
94  utility::vector1<utility::vector1<Real> >& B,
95  utility::vector1<utility::vector1<Real> >& C,
96  utility::vector1<utility::vector1<Real> >& D,
97  utility::vector1<Real>& cal,
98  utility::vector1<Real>& sal,
99  int& f) {
100 
101  utility::vector1<Real> ctheta (4), calpha (4), salpha (4), cxi (3), sxi (3), ceta (4), seta (4), cdelta (3), sdelta (3), caceta (4), caseta (4), saceta (4), saseta (4), capeta (4), sapeta (4), cameta (4), sameta (4);
102  utility::vector1<utility::vector1<Real> > L, M, N;
103  utility::vector1<utility::vector1<utility::vector1<Real> > > p (4);
104 
105  // allocate p
106  for ( int i=1; i<=4; i++ ) {
107  p[i].resize(3);
108  for ( int j=1; j<=3; j++ ) {
109  p[i][j].resize(3);
110  }
111  }
112 
113  // allocate polynomial matrices
114  A.resize(3);
115  B.resize(3);
116  C.resize(3);
117  D.resize(3);
118  for ( int i=1; i<=3; i++ ) {
119  for ( int j=1; j<=3; j++ ) {
120  A[i].resize(3);
121  B[i].resize(3);
122  C[i].resize(3);
123  D[i].resize(3);
124  }
125  }
126 
127  // if triangle is feasible, compute the polynomial matrices
128  if ( (std::abs(vb[1] - vb[2]) <= vb[3]) && (vb[1] + vb[2] >= vb[3]) ) {
129  f=1; // we have a feasible triangle
130 
131  for ( int i=1; i<=3; i++ ) {
132  ctheta[i]=std::cos(theta[i]);
133  }
134  triangle(vb, calpha, salpha);
135  cal.resize(3);
136  sal.resize(3);
137  for ( int i=1; i<=3; i++ ) {
138  cal[i]=calpha[i];
139  sal[i]=salpha[i];
140  }
141  sincos(xi, 2, cxi, sxi);
142  sincos(eta, 2, ceta, seta);
143  sincos(delta, 2, cdelta, sdelta);
144  calpha[4] = calpha[1];
145  salpha[4] = salpha[1];
146  ceta[4] = ceta[1];
147  seta[4] = seta[1];
148  ctheta[4] = ctheta[1];
149  for ( int i=1; i<=4; i++ ) {
150  caceta[i] = calpha[i] * ceta[i];
151  caseta[i] = calpha[i] * seta[i];
152  saceta[i] = salpha[i] * ceta[i];
153  saseta[i] = salpha[i] * seta[i];
154  capeta[i] = caceta[i] - saseta[i];
155  sapeta[i] = saceta[i] + caseta[i];
156  cameta[i] = caceta[i] + saseta[i];
157  sameta[i] = saceta[i] - caseta[i];
158  }
159  for ( int i=2; i<=4; i++ ) {
160  Real am, ap, cdsx, sdsx, b, bm, bp, d, dp, dm;
161  am = cxi[i-1]*cameta[i] + ctheta[i];
162  ap = cxi[i-1]*capeta[i] + ctheta[i];
163  cdsx = cdelta[i-1]*sxi[i-1];
164  sdsx = 2*sdelta[i-1]*sxi[i-1];
165  b = 4*cdsx*seta[i];
166  bm = cdsx*sameta[i];
167  bp = cdsx*sapeta[i];
168  d = sdsx*seta[i];
169  dp = sdsx*sapeta[i];
170  dm = sdsx*sameta[i];
171  p[i][3][3]= -am -bm;
172  p[i][3][2]= -d;
173  p[i][3][1]= -ap -bp;
174  p[i][2][3]= -dm;
175  p[i][2][2]= b;
176  p[i][2][1]= -dp;
177  p[i][1][3]= -am +bm;
178  p[i][1][2]= d;
179  p[i][1][1]= -ap +bp;
180  }
181  for ( int i=1; i<=3; i++ ) {
182  for ( int j=1; j<=3; j++ ) {
183  p[1][i][j]=p[4][i][j];
184  }
185  }
186 
187  L.resize(3);
188  M.resize(3);
189  N.resize(3);
190  for ( int i=1; i<=3; i++ ) {
191  L[i].resize(3);
192  M[i].resize(3);
193  N[i].resize(3);
194  }
195  for ( int k=1; k<=3; k++ ) {
196  for ( int i=1; i<=3; i++ ) {
197  L[k][i] = p[order[1]][i][k];
198  M[k][i] = p[order[2]][k][i];
199  N[k][i] = p[order[3]][k][i];
200  }
201  }
202 
203  for ( int i=1; i<=3; i++ ) {
204  for ( int j=1; j<=3; j++ ) {
205  // not using transpose indexing so we don't take transpose in dixon
206  A[i][j] = M[i][2]*N[1][j] - M[i][1]*N[2][j];
207  B[i][j] = M[i][3]*N[1][j] - M[i][1]*N[3][j];
208  C[i][j] = M[i][3]*N[2][j] - M[i][2]*N[3][j];
209  D[i][j] = L[i][j];
210  }
211  }
212  } else {
213  f=0; // we have an unfeasible triangle
214  for ( int i=1; i<=3; i++ ) {
215  for ( int j=1; j<=3; j++ ) {
216  A[i][j]=0;
217  B[i][j]=0;
218  C[i][j]=0;
219  D[i][j]=0;
220  }
221  }
222  }
223 }
224 
225 /*
226 * cross
227 * -----
228 * Output r is cross product of vectors L and r0
229 */
230 void cross(const utility::vector1<Real>& L, const utility::vector1<Real>& r0, utility::vector1<Real>& r) {
231  r.resize(3);
232  r[1]=L[2]*r0[3]-L[3]*r0[2];
233  r[2]=L[3]*r0[1]-L[1]*r0[3];
234  r[3]=L[1]*r0[2]-L[2]*r0[1];
235 }
236 
237 /*
238 * frame
239 * -----
240 * construct body frame for r1 origin
241 * r2 r1-->r2: positive x axis
242 * r3 on pos. y side of xy plane
243 */
244 void frame(const utility::vector1<utility::vector1<Real> >& R, utility::vector1<utility::vector1<Real> >& U) {
245  Real norm1, norm3;
246 
247  // utility::vector1<Real> cross1 (3), cross3 (3), cross11 (3), cross12 (3), cross31 (3), cross32 (3);
248  utility::vector1<Real> dR3R2 (3);
249  U.resize(3);
250  for ( int i=1; i<=3; i++ ) {
251  U[i].resize(3);
252  }
253  for ( int i=1; i<=3; i++ ) {
254  U[1][i]=R[2][i]-R[1][i];
255  }
256  norm1=std::sqrt(U[1][1]*U[1][1] + U[1][2]*U[1][2] + U[1][3]*U[1][3]);
257  for ( int i=1; i<=3; i++ ) {
258  U[1][i] /= norm1;
259  dR3R2[i]=R[3][i]-R[2][i];
260  }
261  cross(U[1], dR3R2, U[3]);
262  norm3=std::sqrt(U[3][1]*U[3][1] + U[3][2]*U[3][2] + U[3][3]*U[3][3]);
263  for ( int i=1; i<=3; i++ ) {
264  U[3][i] /= norm3;
265  }
266  cross(U[3], U[1], U[2]);
267 }
268 
269 Real eucDistance(const utility::vector1<Real>& a, const utility::vector1<Real>& b) {
270  return std::sqrt(std::pow(a[1]-b[1],2) + std::pow(a[2]-b[2],2) + std::pow(a[3]-b[3],2));
271 }
272 
273 Real scpn(const utility::vector1<Real>& a, const utility::vector1<Real>& b, const utility::vector1<Real>& c) {
274  Real d=0;
275  for ( int i=1; i<=3; i++ ) {
276  d += (a[i]-b[i]) * (c[i]-b[i]);
277  }
278  return d;
279 }
280 
281 Real bondangle(const utility::vector1<Real>& a, const utility::vector1<Real>& b, const utility::vector1<Real>& c) {
282  Real r = scpn(a,b,c) / (eucDistance(a,b) * eucDistance(b,c));
283  Real ang = std::atan2(sqrt(1-r*r), (Real) r) * rad2deg;
284  return ang;
285 }
286 
287 Real torsion(const utility::vector1<Real>& a, const utility::vector1<Real>& b, const utility::vector1<Real>& c, const utility::vector1<Real>& d) {
288  Real chi;
289  utility::vector1<Real> r (3), sc1 (3), sc2 (3), sc3 (3), cs12, cs31;
290  utility::vector1<utility::vector1<Real> > s (3);
291  for ( int i=1; i<=3; i++ ) { // init s
292  s[i].resize(3);
293  }
294  if ( b==c ) {
295  chi=pidegs;
296  } else {
297  for ( int i=1; i<=3; i++ ) {
298  r[i] = d[i] - b[i];
299  s[i][1] = c[i] - b[i];
300  s[i][2] = a[i] - b[i];
301  sc1[i] = s[i][1];
302  sc2[i] = s[i][2];
303  }
304  cross(sc1, sc2, cs12);
305  for ( int i=1; i<=3; i++ ) {
306  s[i][3] = cs12[i];
307  sc3[i] = s[i][3];
308  }
309  cross(sc3, sc1, cs31);
310  for ( int i=1; i<=3; i++ ) {
311  s[i][2] = cs31[i];
312  }
313  for ( int i=2; i<=3; i++ ) {
314  Real f = std::sqrt( pow(s[1][i],2) + pow(s[2][i],2) + pow(s[3][i],2));
315  s[1][i] /= f;
316  s[2][i] /= f;
317  s[3][i] /= f;
318  }
319  Real y = r[1]*s[1][2] + r[2]*s[2][2] + r[3]*s[3][2];
320  Real z = r[1]*s[1][3] + r[2]*s[2][3] + r[3]*s[3][3];
321  chi = atan2(z,y) * rad2deg;
322  if ( chi < 0 ) {
323  chi += 360.0;
324  }
325  }
326  return chi;
327 }
328 
329 void chainParams(const int& n, const utility::vector1<utility::vector1<Real> >& atoms, Real& vbond, Real& xi, Real& eta, Real& delta, utility::vector1<Real>& R0, utility::vector1<utility::vector1<Real> >& Q) {
330  //utility::vector1<Real> ac1 (3), ac2 (3), acn (3), acneg1 (3);
331  utility::vector1<utility::vector1<Real> > a1n2 (3);
332  // R0.resize(3);
333  // for (int i=1; i<=3; i++) {
334  // ac1[i] = atoms[i][1];
335  // ac2[i] = atoms[i][2];
336  // acn[i] = atoms[i][n];
337  // acneg1[i]= atoms[i][n-1];
338  // R0[i] = ac1[i];
339  // }
340  a1n2[1]=atoms[1];
341  a1n2[2]=atoms[n];
342  a1n2[3]=atoms[2];
343 
344  // printVector(acn);
345  // printVector(ac1);
346 
347  R0=atoms[1];
348  frame(a1n2, Q);
349  vbond=eucDistance(atoms[n], atoms[1]);
350  xi=bondangle(atoms[n-1], atoms[n], atoms[1]);
351  eta=bondangle(atoms[n], atoms[1], atoms[2]);
352  delta=torsion(atoms[n-1], atoms[n], atoms[1], atoms[2]);
353 
354  // frame(a1n2, Q);
355  // vbond=eucDistance(acn, ac1);
356  // xi=bondangle(acneg1, acn, ac1);
357  // eta=bondangle(acn, ac1, ac2);
358  // delta=torsion(acneg1, acn, ac1, ac2);
359 }
360 
361 
362 /* Generates a set of cartesian coordinates (atoms) given the supplied bond angles, bond lengths, and torsions
363 * of a peptide chain
364 */
365 void chainXYZ(const int& n,
366  const utility::vector1<Real>& b_len1,
367  const utility::vector1<Real>& b_ang1,
368  const utility::vector1<Real>& t_ang1,
369  const bool space,
370  const utility::vector1<Real>& R0,
371  const utility::vector1<utility::vector1<Real> >& Q,
372  utility::vector1<utility::vector1<Real> >& atoms
373 )
374 {
375  Real ca, sa, ct, st; // cos(b_ang[2]), sin(b_ang[2]), cos(t_ang[2]), sin(t_ang[2]),
376  Real ba;
377  int n1, n2, j1, j2;
378  utility::vector1<Real> b_len, b_ang (n), t_ang (n);//, cos_a (n), cos_t (n), sin_a (n), sin_t (n);
379  utility::vector1<Real> v(3), s(3);
380  utility::vector1<utility::vector1<Real> > U (3); // rotation matrix
381 
382  //std::cout << "about to allocate atoms" << std::endl;
383  atoms.resize(n);
384  b_len = b_len1; // only defined for consistent naming
385  // convert angles to radians and allocate atom vectors
386  for ( int i=1; i<=n; i++ ) {
387  t_ang[i] = deg2rad*t_ang1[i];
388  b_ang[i] = deg2rad*b_ang1[i];
389  atoms[i].resize(3);
390  }
391  //std::cout << "atoms allocated." << std::endl;
392  // place the first atom at the origin
393  for ( int i=1; i<=n; i++ ) {
394  for ( int j=1; j<=3; j++ ) {
395  atoms[i][j]=0.0;
396  }
397  }
398  //std::cout << "placed first atom" << std::endl;
399  //printMatrix(atoms);
400  // place the second atom at bond length along x-axis from first atom
401  atoms[2][1] = b_len[1];
402  // third atom shows that b_ang is interior
403  ca = std::cos(b_ang[2]);
404  sa = std::sin(b_ang[2]);
405  atoms[3][1] = b_len[1] - b_len[2] * ca;
406  atoms[3][2] = b_len[2] * sa;
407  //std::cout << "placed second and third atoms" << std::endl;
408  //printMatrix(atoms);
409 
410  // initialize rotation matrix
411  for ( int i=1; i<=3; i++ ) {
412  U[i].resize(3);
413  }
414  U[1][1] = -ca;
415  U[2][1] = -sa;
416  U[1][2] = sa;
417  U[2][2] = -ca;
418  for ( int i=1; i<=2; i++ ) {
419  U[3][i] = 0;
420  U[i][3] = 0;
421  }
422  U[3][3] = 1;
423 
424  // all other atoms
425  n1 = n-1;
426  n2 = n-2;
427  j1 = 3;
428  j2 = 2;
429  //std::cout << "preparing to place all other atoms" << std::endl;
430  //std::cout << "U before for loop: " << std::endl;
431  //printMatrix(U);
432  //exit(0);
433  for ( int j=4; j<=n1; j++ ) {
434  ca = std::cos(b_ang[j1]);
435  sa = std::sin(b_ang[j1]);
436  ct = std::cos(t_ang[j2]);
437  st = std::sin(t_ang[j2]);
438  for ( int i=1; i<=3; i++ ) {
439  s[i] = U[2][i];
440  }
441  for ( int i=1; i<=3; i++ ) {
442  U[2][i] = s[i]*ct + U[3][i]*st;
443  }
444  for ( int i=1; i<=3; i++ ) {
445  U[3][i] = -s[i]*st + U[3][i]*ct;
446  }
447  for ( int i=1; i<=3; i++ ) {
448  s[i] = U[2][i];
449  }
450  for ( int i=1; i<=3; i++ ) {
451  U[2][i] = -s[i]*ca - U[1][i]*sa;
452  }
453  for ( int i=1; i<=3; i++ ) {
454  U[1][i] = s[i]*sa - U[1][i]*ca;
455  }
456  for ( int i=1; i<=3; i++ ) {
457  atoms[j][i] = atoms[j1][i] + b_len[j1] * U[1][i];
458  }
459  j2 = j1;
460  j1 = j;
461  //std::cout << "at end of j=" << j << std::endl;
462  //printMatrix(atoms);
463  }
464  //std::cout << "ended atom placement for loop. printing U" << std::endl;
465  //printMatrix(U);
466  //exit(0);
467  ca = std::cos(b_ang[n1]);
468  sa = std::sin(b_ang[n1]);
469  ct = std::cos(t_ang[n2]);
470  st = std::sin(t_ang[n2]);
471  v[1] = -b_len[n1] * ca;
472  ba = b_len[n1] * sa;
473  v[2] = ba * ct;
474  v[3] = ba * st;
475  //std::cout << "about to calculate s = U*v'. first printing v" << std::endl;
476  //printVector(v);
477  // implementing s = U*v'
478  for ( int i=1; i<=3; i++ ) {
479  s[i] = (U[1][i] * v[1]) + (U[2][i] * v[2]) + (U[3][i] * v[3]);
480  //std::cout << "check 1" << std::endl;
481  }
482  for ( int i=1; i<=3; i++ ) {
483  atoms[n][i] = atoms[n1][i] + s[i];
484  //std::cout << "check 2" << std::endl;
485  }
486  //printMatrix(atoms);
487  // check atoms here but still need to implement space and frame
488  if ( space == true ) {
489  // implement:
490  // for i = 1:n
491  // atoms(:,i) = R0 + Q*atoms(:,i);
492  // end
493  utility::vector1<Real> spaceatoms (3); // holder for inner-loop coordinates
494  for ( int ind=1; ind<=n; ind++ ) {
495  for ( int i=1; i<=3; i++ ) {
496  Real inner_sum=0.0;
497  for ( int j=1; j<=3; j++ ) {
498  inner_sum += (Q[i][j] * atoms[ind][j]);
499  }
500  spaceatoms[i]=R0[i]+inner_sum;
501  }
502  for ( int k=1; k<=3; k++ ) {
503  atoms[ind][k]=spaceatoms[k];
504  }
505  }
506  }
507  /*
508  else if (index[1] != 1 || index[2] != 2 || index[3] != 3) {
509  std::cout << "setting to body frame based on index" << std::endl;
510  utility::vector1<Real> Ra (3);
511  utility::vector1<utility::vector1<Real> > RR (3), Qa (3), body_atoms(n);
512  //std::cout << "atoms is " << std::endl;
513  //printMatrix(atoms);
514  for (int i=1; i<=3; i++) {
515  RR[i].resize(3);
516  Ra[i] = -atoms[index[1]][i];
517  }
518  //std::cout << "Ra is " << std::endl;
519  //printVector(Ra);
520  for (int i=1; i<=3; i++) {
521  for (int j=1; j<=3; j++) {
522  RR[j][i]=atoms[index[j]][i];
523  }
524  }
525  //std::cout << "RR is " << std::endl;
526  //printMatrix(RR);
527  frame(RR, Qa);
528  //std::cout << "Qa is " << std::endl;
529  //printMatrix(Qa);
530  for (int ind=1; ind<=n; ind++) {
531  // implement atoms(:,i) = Qa'*(atoms(:,i)+Ra);
532  body_atoms[ind].resize(3);
533  for (int i=1; i<=3; i++) {
534  Real inner_sum=0.0;
535  for (int j=1; j<=3; j++) {
536  inner_sum += Qa[i][j] * (atoms[ind][j] + Ra[j]);
537  //std::cout << "Qa[" << i << "][" << j << "]: " << Qa[i][j] << std::endl;
538  //std::cout << "atoms[" << ind << "][" << j << "]: " << atoms[ind][j] << std::endl;
539  }
540  body_atoms[ind][i] = inner_sum;
541  }
542  }
543  atoms = body_atoms;
544  //exit(0);
545  }
546  */
547 }
548 
549 void chainXYZ(
550  const int& n,
551  const utility::vector1<Real>& b_len,
552  const utility::vector1<Real>& b_ang,
553  const utility::vector1<Real>& t_ang,
554  utility::vector1<utility::vector1<Real> >& atoms) {
555 
556  utility::vector1<Real> dummy_R0;
557  utility::vector1<utility::vector1<Real> > dummy_Q;
558 
559  chainXYZ(n, b_len, b_ang, t_ang, false, dummy_R0, dummy_Q, atoms);
560 }
561 
562 /* chainTORS
563 // ---------
564 // calculates torsions, bond angles, bond lengths from atoms
565 // for a linear chain
566 // calculates angles in degrees
567 //
568 // (j-1) -------- (j)
569 // b_ang(j) \
570 // t_ang(j) \b_len(j)
571 // \
572 // (j+1) ------ (j+2)
573 //
574 */
575 
576 void chainTORS (const int& n, const utility::vector1<utility::vector1<Real> >& atoms, utility::vector1<Real>& t_ang, utility::vector1<Real>& b_ang, utility::vector1<Real>& b_len, utility::vector1<Real>& R0, utility::vector1<utility::vector1<Real> >& Q) {
577  t_ang.resize(n);
578  b_ang.resize(n);
579  b_len.resize(n);
580  utility::vector1<utility::vector1<Real> > framein (3);
581  for ( int i=1; i<=3; i++ ) {
582  framein[i]=atoms[i];
583  }
584  R0=atoms[1];
585  frame(framein,Q);
586  for ( int j=1; j <= n-1; j++ ) {
587  b_len[j] = eucDistance( atoms[j], atoms[j+1] );
588  }
589  b_len[n] = eucDistance( atoms[n], atoms[1] );
590  for ( int j=2; j <= n-1; j++ ) {
591  b_ang[j] = bondangle( atoms[j-1], atoms[j], atoms[j+1] );
592  }
593  b_ang[n] = bondangle( atoms[n-1], atoms[n], atoms[1] );
594  b_ang[1] = bondangle( atoms[n], atoms[1], atoms[2] );
595 
596  if ( n >= 4 ) {
597  for ( int j=2; j <= n-2; j++ ) {
598  t_ang[j] = torsion( atoms[j-1], atoms[j], atoms[j+1], atoms[j+2] );
599  }
600  t_ang[n-1] = torsion( atoms[n-2], atoms[n-1], atoms[n], atoms[1] );
601  t_ang[n] = torsion( atoms[n-1], atoms[n], atoms[1], atoms[2] );
602  t_ang[1] = torsion( atoms[n] , atoms[1], atoms[2], atoms[3] );
603  } else {
604  for ( int j=1; j<=n; j++ ) {
605  t_ang[j]=0.0;
606  }
607  }
608 }
609 
610 void rotateX(const utility::vector1<utility::vector1<Real> >& R, const Real& c, const Real& s, utility::vector1<utility::vector1<Real> >& Rx) {
611  int dim1=R.size();
612  Rx.resize(dim1);
613  for ( int i=1; i<=dim1; i++ ) {
614  Rx[i].resize(R[i].size());
615  Rx[i][1] = R[i][1];
616  Rx[i][2] = c*R[i][2] - s*R[i][3];
617  Rx[i][3] = s*R[i][2] + c*R[i][3];
618  }
619 }
620 
621 void rotateY(const utility::vector1<utility::vector1<Real> >& R, const Real& c, const Real& s, utility::vector1<utility::vector1<Real> >& Ry) {
622  int dim1=R.size();
623  Ry.resize(dim1);
624  for ( int i=1; i<=dim1; i++ ) {
625  Ry[i].resize(R[i].size());
626  Ry[i][1] = c*R[i][1] + s*R[i][3];
627  Ry[i][2] = R[i][2];
628  Ry[i][3] = -s*R[i][1] + c*R[i][3];
629  }
630 }
631 
632 void rotateZ(const utility::vector1<utility::vector1<Real> >& R, const Real& c, const Real& s, utility::vector1<utility::vector1<Real> >& Rz) {
633  int dim1=R.size();
634  Rz.resize(dim1);
635  for ( int i=1; i<=dim1; i++ ) {
636  Rz[i].resize(R[i].size());
637  Rz[i][1] = c*R[i][1] - s*R[i][2];
638  Rz[i][2] = s*R[i][1] + c*R[i][2];
639  Rz[i][3] = R[i][3];
640  }
641 }
642 
643 
644 ////////////////////////////////////////////////////////////////////////////////
645 ///
646 /// @brief
647 /// Solve the triaxial loop closure problem for a system of atoms
648 ///
649 /// @details
650 ///
651 /// @param[in] atoms - matrix of cartesian coordiantes of N-CA-C atoms indexed as atoms[atom][dimension]
652 /// @param[in] dt - desired torsions for each atom
653 /// @param[in] da - desired bond angle for each atom
654 /// @param[in] db - desired bond length for each atom
655 /// @param[in] pivots - 3 indices (base 1) of atoms to be used as loop closure pivots
656 /// @param[in] order - length 3 vector giving order to solve for the tau parameters (use 1,2,3 if unsure)
657 /// @param[out] t_ang - matrix giving torsion angles for each atom for each solution, indexed as t_ang[solution][atom]
658 /// @param[out] b_ang - matrix giving bond angles for each atom for each solution, indexed as b_ang[solution][atom]
659 /// @param[out] b_len - matrix giving bond lengths for each atom for each solution, indexed as b_len[solution][atom]
660 /// @param[out] nsol - number of solutions found
661 ///
662 /// @global_read
663 ///
664 /// @global_write
665 ///
666 /// @return
667 ///
668 /// @remarks
669 /// dt, da, and db are cast to Real precision when placed in t_ang1,2, b_ang1,2 and b_len1,2.
670 /// Solution is carried out in Real precision.
671 /// Solutions are cast back to Reals when placed into t_ang.
672 ///
673 /// @references
674 ///
675 /// @author Evangelos A. Coutsias
676 /// @author Daniel J. Mandell
677 ///
678 ////////////////////////////////////////////////////////////////////////////////
679 
680 void bridgeObjects (const utility::vector1<utility::vector1<Real> >& atoms,
681  const utility::vector1<Real> & dt,
682  const utility::vector1<Real> & da,
683  const utility::vector1<Real> & db,
684  const utility::vector1<int>& pivots,
685  const utility::vector1<int>& order,
686  utility::vector1<utility::vector1<Real> >& t_ang,
687  utility::vector1<utility::vector1<Real> >& b_ang,
688  utility::vector1<utility::vector1<Real> >& b_len,
689  int& nsol)
690 {
691  utility::vector1<Real> R0 (3) /* for chainXYZ interface */, R1 (3), R2 (3), R3 (3); // origins of triangle pieces
692  utility::vector1<utility::vector1<Real> > Q0; // for chainXYZ interface
693  utility::vector1<utility::vector1<Real> > chain1, chain2, chain3, fchain1, fchain2, fchain3; // 3 triangle pieces
694  utility::vector1<utility::vector1<Real> > chain1a, chain1b, chain2a, chain2b, chain12, chain12rX; // reconstructions of the molecular chain
695  utility::vector1<Real> t_ang1, b_ang1, b_len1, t_ang2, b_ang2, b_len2; // torsions, angles, lengths of chains 1,2
696  utility::vector1<Real> vbond (3), xi (3), eta (3), delta (3), theta (3); // triaxial parameters for each pivot
697  utility::vector1<Real> cal, sal;
698  utility::vector1<utility::vector1<Real> > frame1, frame2, frame3;
699  utility::vector1<utility::vector1<Real> > A, B, C, D; // dixon matrices
700  utility::vector1<utility::vector1<Real> > cosines, sines, taus;
701  utility::vector1<utility::vector1<utility::vector1<Real> > > loop; // the solutions
702  int k1, k2, k3, l1, l2, l3, l3a, l3b, f;
703  int ind;
704  int N=atoms.size(); // number of atoms in chain
705 
706  k1 = pivots[1];
707  k2 = pivots[2];
708  k3 = pivots[3];
709  l1 = k2-k1+1;
710  l2 = k3-k2+1;
711  l3a = N-k3+1;
712  l3b = k1;
713  l3 = l3a + l3b;
714  chain1.resize(l1);
715  chain2.resize(l2);
716  chain3.resize(l3);
717  t_ang1.resize(l1);
718  t_ang2.resize(l2);
719  b_ang1.resize(l1);
720  b_ang2.resize(l2);
721  b_len1.resize(l1);
722  b_len2.resize(l2);
723 
724  // break the designed torsions, angles, and lengths into 3 chains (3rd chain's coords set directly)
725  ind=1;
726  for ( int i=pivots[1]; i<=pivots[2]; i++ ) {
727  t_ang1[ind++]=dt[i];
728  }
729  ind=1;
730  for ( int i=pivots[2]; i<=pivots[3]; i++ ) {
731  t_ang2[ind++]=dt[i];
732  }
733  ind=1;
734  for ( int i=pivots[1]; i<=pivots[2]; i++ ) {
735  b_ang1[ind++]=da[i];
736  }
737  ind=1;
738  for ( int i=pivots[2]; i<=pivots[3]; i++ ) {
739  b_ang2[ind++]=da[i];
740  }
741  ind=1;
742  for ( int i=pivots[1]; i<=pivots[2]; i++ ) {
743  b_len1[ind++] = db[i];
744  }
745  ind=1;
746  for ( int i=pivots[2]; i<=pivots[3]; i++ ) {
747  b_len2[ind++] = db[i];
748  }
749 
750  chainXYZ( l1, b_len1, b_ang1, t_ang1, false, R0, Q0, chain1 );
751  chainXYZ( l2, b_len2, b_ang2, t_ang2, false, R0, Q0, chain2 );
752 
753  // set chain 3 directly from original atoms
754  ind=1;
755  for ( int i=k3; i<=N; i++ ) {
756  chain3[ind++] = atoms[i];
757  }
758  for ( int i=1; i<=k1; i++ ) {
759  chain3[ind++] = atoms[i];
760  }
761 
762  /*
763  std::cout << "chain 1: " << std::endl;
764  printMatrix(chain1);
765  std::cout << "chain 2: " << std::endl;
766  printMatrix(chain2);
767  std::cout << "chain 3: " << std::endl;
768  printMatrix(chain3);
769  */
770 
771  // compute pivot traingle parameters
772  chainParams( l1, chain1, vbond[1], xi[1], eta[1], delta[1], R1, frame1 );
773  for ( int i=1; i<=l1; i++ ) {
774  for ( int j=1; j<=3; j++ ) {
775  chain1[i][j] -= R1[j];
776  }
777  }
778  multTransMatrix(frame1, chain1, fchain1);
779  chainParams( l2, chain2, vbond[2], xi[2], eta[2], delta[2], R2, frame2 );
780  for ( int i=1; i<=l2; i++ ) {
781  for ( int j=1; j<=3; j++ ) {
782  chain2[i][j] -= R2[j];
783  }
784  }
785  multTransMatrix(frame2, chain2, fchain2);
786  chainParams( l3, chain3, vbond[3], xi[3], eta[3], delta[3], R3, frame3 );
787  for ( int i=1; i<=3; i++ ) {
788  theta[i] = deg2rad * da[pivots[i]]; // use the designed bond angles
789  eta[i] *= deg2rad;
790  xi[i] *= deg2rad;
791  delta[i] *= deg2rad;
792  }
793  fchain3=chain3; // for naming consistency
794 
795  // check feasibility of arrangement and define polynomial coefficients
796  triaxialCoefficients( vbond, xi, eta, delta, theta, order, A, B, C, D, cal, sal, f );
797 
798  if ( f == 0 ) {
799  nsol = 0;
800  return;
801  }
802 
803  // find the solutions
804  dixon_eig( A, B, C, D, order, cosines, sines, taus, nsol );
805 
806  // reconstruct the molecular chain with new torsions
807  loop.resize(nsol);
808  t_ang.resize(nsol);
809  b_ang.resize(nsol);
810  b_len.resize(nsol);
811  for ( int j=1; j<=nsol; j++ ) {
812  rotateX(fchain1, cosines[j][1], sines[j][1], chain1a);
813  rotateZ(chain1a, cal[1], sal[1], chain1b);
814  for ( unsigned k=1; k<=chain1b.size(); k++ ) {
815  chain1b[k][1] += vbond[3];
816  }
817  rotateX(fchain2, cosines[j][2], sines[j][2], chain2a);
818  rotateZ(chain2a, cal[3], -sal[3], chain2b);
819  for ( unsigned k=1; k<=chain2b.size(); k++ ) {
820  chain2b[k][1] -= chain2b[l2][1];
821  chain2b[k][2] -= chain2b[l2][2];
822  }
823  int chain12len = (l1-1) + (l2-2); // implementing chain12 = [chain1b(:,2:l1), chain2b(:,2:l2-1)];
824  chain12.resize(chain12len);
825  int ind = 1;
826  for ( int k=2; k<=l1; k++ ) {
827  chain12[ind].resize(3);
828  for ( int n=1; n<=3; n++ ) {
829  chain12[ind][n] = chain1b[k][n];
830  }
831  ind++;
832  }
833  for ( int k=2; k<=l2-1; k++ ) {
834  chain12[ind].resize(3);
835  for ( int n=1; n<=3; n++ ) {
836  chain12[ind][n] = chain2b[k][n];
837  }
838  ind++;
839  }
840  rotateX(chain12, cosines[j][3], -sines[j][3], chain12rX);
841  multMatrix(frame3, chain12rX, chain12);
842  for ( int k=1; k<=chain12len; k++ ) {
843  for ( int n=1; n<=3; n++ ) {
844  chain12[k][n] += R3[n];
845  }
846  }
847 
848  // populate the loop table with the cartesian solutions
849  loop[j].resize(N);
850  for ( int k=1; k<=l3b; k++ ) {
851  loop[j][k].resize(3);
852  for ( int kk=1; kk<=3; kk++ ) {
853  loop[j][k][kk] = chain3[l3a+k][kk];
854  }
855  }
856  for ( int k=1; k<= l1 + l2 - 3; k++ ) {
857  loop[j][l3b+k].resize(3);
858  for ( int kk=1; kk<=3; kk++ ) {
859  loop[j][l3b+k][kk] = chain12[k][kk];
860  }
861  }
862  for ( int k=1; k<=l3a; k++ ) {
863  loop[j][l3b+l1+l2-3+k].resize(3);
864  for ( int kk=1; kk<=3; kk++ ) {
865  loop[j][l3b+l1+l2-3+k][kk] = chain3[k][kk];
866  }
867  }
868 
869  // calculate the 6 closure torsions and place them (with Real precision) into the output vector
870  // together with the designed torsions, angles, and lengths (which are copied directly from the input)
871 
872  int pivots2m1 = pivots[2]-1, pivots2p1 = pivots[2]+1,
873  pivots3m1 = pivots[3]-1, pivots3p1 = pivots[3]+1;
874 
875  t_ang[j] = dt;
876  b_ang[j] = da;
877  b_len[j] = db;
878  t_ang[j][4] = static_cast<Real> (torsion(loop[j][3],loop[j][4],loop[j][5],loop[j][6])); // pivot 1 always atom 5
879  t_ang[j][5] = static_cast<Real> (torsion(loop[j][4],loop[j][5],loop[j][6],loop[j][7]));
880  t_ang[j][pivots2m1] = static_cast<Real> (torsion(loop[j][pivots[2]-2],loop[j][pivots2m1],
881  loop[j][pivots[2]],loop[j][pivots2p1]));
882  t_ang[j][pivots[2]] = static_cast<Real> (torsion(loop[j][pivots2m1],loop[j][pivots[2]],
883  loop[j][pivots2p1],loop[j][pivots[2]+2]));
884  t_ang[j][pivots3m1] = static_cast<Real> (torsion(loop[j][pivots[3]-2],loop[j][pivots3m1],
885  loop[j][pivots[3]],loop[j][pivots3p1]));
886  t_ang[j][pivots[3]] = static_cast<Real> (torsion(loop[j][pivots3m1],loop[j][pivots[3]],
887  loop[j][pivots3p1],loop[j][pivots[3]+2]));
888  }
889 }
890 
891 
892 void test_bridgeObjects() {
893  // input arguments
894  int N=8;
895  utility::vector1<Real> t_ang0 (N), b_ang0 (N), b_len0 (N), dt (N), db (N), da (N);
896  utility::vector1<int> pivots (3), order (3);
897  // input argument values
898  // DJM [ if inputting torsions, uncomment next 4 lines ]
899  Real t_ang0vals[] = {36.0213, 0.0, 253.8843, 52.8441, 36.0366, 359.8985, 253.9602, 52.8272};
900  Real b_ang0vals[] = {114.9908, 114.9805, 115.0112, 114.9923, 114.9582, 115.0341, 114.9920, 115.0594};
901  Real b_len0vals[] = {1.5200, 1.5202, 1.5197, 1.5205, 1.5198, 1.5201, 1.5197, 1.5196};
902  Real dtvals[] = {180.0000, 180.0000, -5.3651, 180.0000, 180.0000, 180.0000, 180.0000, -5.0393};
903 
904  utility::vector1<utility::vector1<Real> > atoms (N);
905  // output arguments
906  int nsol;
907  utility::vector1<utility::vector1<Real> > t_ang, b_ang, b_len;
908  /*
909  * DJM uncomment if inputting torsions
910  */
911  int ind=0;
912  for ( int i=1; i<=8; i++ ) {
913  t_ang0[i]=t_ang0vals[ind++];
914  }
915 
916  ind=0;
917  for ( int i=1; i<=8; i++ ) {
918  b_ang0[i]=b_ang0vals[ind++];
919  }
920 
921  ind=0;
922  for ( int i=1; i<=8; i++ ) {
923  b_len0[i]=b_len0vals[ind++];
924  }
925 
926  ind=0;
927  for ( int i=1; i<=8; i++ ) {
928  dt[i]=dtvals[ind++];
929  }
930 
931  for ( int i=1; i<=8; i++ ) {
932  db[i]=1.52;
933  }
934 
935  for ( int i=1; i<=8; i++ ) {
936  da[i]=115;
937  }
938  //*/
939  pivots[1] = 2;
940  pivots[2] = 5;
941  pivots[3] = 7;
942  order[1] = 1;
943  order[2] = 2;
944  order[3] = 3;
945 
946  //bridgeObjects(N, t_ang0, b_ang0, b_len0, dt, db, da, pivots, order, t_ang, b_ang, b_len, nsol);
947  bridgeObjects(atoms, dt, db, da, pivots, order, t_ang, b_ang, b_len, nsol);
948  //bridgeObjects(atoms, pivots, order, t_ang, b_ang, b_len, nsol);
949  //printMatrix(t_ang);
950  //printMatrix(b_ang);
951  //printMatrix(b_len);
952  ////std::cout << nsol << std::endl;
953 }
954 
955 /* // DJM: this was for benchmarking current version of bridgeObjects against earlier (slower) version.
956 // v2 is now current version. 6-8-08
957 void test_bridgeObjects_v2() {
958 // input arguments
959 int N=27;
960 utility::vector1<Real> dt (N), db (N), da (N);
961 utility::vector1<int> pivots (3), order (3);
962 // input argument values
963 // DJM [ if inputting torsions, uncomment next 4 lines ]
964 Real dt_ang_vals[] = {348.798, 337.654, 181.609, -170, 140, 179.743, -140, 80, 181.769, -130, 10, 180.419, -90, -30, 177.171, -90, -20, 184.154, -90, -40, 182.472, -80, 170, 184.776, 310.274, 150.675, 290.245};
965 Real db_ang_vals[] = {88.2267, 112.669, 115.978, 122.109, 110.442, 117.607, 121.606, 110.194, 115.415, 120.286, 110.888, 117.314, 125.907, 110.287, 117.773, 121.692, 109.691, 115.437, 122.169, 109.827, 119.212, 120.735, 111.479, 115.289, 121.921, 114.467, 85.1688};
966 Real db_len_vals[] = {1.41632, 1.48472, 1.32645, 1.46204, 1.54295, 1.31506, 1.45302, 1.51076, 1.28081, 1.46153, 1.51705, 1.33131, 1.46801, 1.5172, 1.34797, 1.45556, 1.49902, 1.31172, 1.41681, 1.49235, 1.36072, 1.4561, 1.52821, 1.31126, 1.44372, 1.48442, 9.90623};
967 // initialize the atoms vector
968 utility::vector1<utility::vector1<Real> > atoms (N);
969 Real atoms_vals[] = {13.92, 13.881, 12.607, 11.925, 10.671, 9.481, 9.657, 8.597, 8.771, 8.01, 8.116, 6.871, 5.803, 4.504, 3.703, 2.616, 1.752, 0.851, 0.206, -0.747, -0.251, 0.965, 1.458, 2.652, 2.912, 4.066, 4.232, 43.318, 44.382, 44.429, 43.293, 43.179, 43.781, 44.046, 44.619, 46.119, 46.785, 48.241, 48.954, 48.205, 48.696, 47.567, 47.894, 46.878, 46.229, 47.052, 46.598, 45.73, 45.989, 45.197, 44.279, 44.041, 43.277, 41.984, 25.637, 24.703, 23.942, 23.88, 23.137, 23.913, 25.189, 26.001, 26.047, 25.261, 25.191, 25.684, 25.95, 26.426, 27.047, 27.774, 28.357, 27.35, 26.558, 25.613, 24.505, 23.952, 22.834, 23.093, 24.356, 24.767, 24.057};
970 for (int i=1; i<=N; i++) {
971 atoms[i].resize(3);
972 }
973 int ind=0;
974 for (int j=1; j<=3; j++) {
975 for (int i=1; i<=N; i++) {
976 atoms[i][j] = atoms_vals[ind++];
977 }
978 }
979 ind=0;
980 // initialize the designed torsions, angles, and length
981 for (int i=1; i<=N; i++) {
982 dt[i]=dt_ang_vals[ind++];
983 }
984 ind=0;
985 for (int i=1; i<=N; i++) {
986 da[i]=db_ang_vals[ind++];
987 }
988 ind=0;
989 for (int i=1; i<=N; i++) {
990 db[i]=db_len_vals[ind++];
991 }
992 // these pivots specifically correspond to the designed inputs to produce the expected outputs
993 pivots[1] = 5;
994 pivots[2] = 14;
995 pivots[3] = 23;
996 order[1] = 1;
997 order[2] = 2;
998 order[3] = 3;
999 
1000 // DJM: debug
1001 //std::cout << "dt: " << std::endl;
1002 //printVector(dt);
1003 //std::cout << "db: " << std::endl;
1004 //printVector(db);
1005 //std::cout << "da: " << std::endl;
1006 //printVector(da);
1007 // output arguments
1008 int nits=100000; // number times to call bridgeObjects for performance benchmarking
1009 int nsol;
1010 utility::vector1<utility::vector1<Real> > t_ang, b_ang, b_len;
1011 
1012 std::cout << "calling bridge2" << std::endl;
1013 std::time_t const start_time_v2( std::time( NULL ) );
1014 for (int d=1; d<=nits; d++) {
1015 bridgeObjects_v2(atoms, dt, da, db, pivots, order, t_ang, b_ang, b_len, nsol);
1016 }
1017 std::time_t const end_time_v2( std::time( NULL ) );
1018 std::cout << "start time: " << start_time_v2 << std::endl;
1019 std::cout << "end time: " << end_time_v2 << std::endl;
1020 std::cout << "start time minus end time v2: " << end_time_v2 - start_time_v2 << std::endl;
1021 std::cout << "from bridgeObjects v2: " << std::endl;
1022 printMatrix(t_ang);
1023 
1024 nsol=0;
1025 
1026 std::cout << "calling bridge1" << std::endl;
1027 std::time_t const start_time_v1( std::time( NULL ) );
1028 for (int d=1; d<=nits; d++) {
1029 bridgeObjects(atoms, dt, da, db, pivots, order, t_ang, b_ang, b_len, nsol);
1030 }
1031 std::time_t const end_time_v1( std::time( NULL ) );
1032 std::cout << "start time: " << start_time_v1 << std::endl;
1033 std::cout << "end time: " << end_time_v1 << std::endl;
1034 std::cout << "start time minus end time v1: " << end_time_v1 - start_time_v1 << std::endl;
1035 std::cout << "from bridgeObjects v1: " << std::endl;
1036 printMatrix(t_ang);
1037 
1038 //std::cout << "output t_ang: " << std::endl;
1039 //printMatrix(t_ang);
1040 //std::cout << "output b_ang: " << std::endl;
1041 //printMatrix(b_ang);
1042 //std::cout << "output b_len: " << std::endl;
1043 //printMatrix(b_len);
1044 ////std::cout << nsol << std::endl;
1045 }
1046 */
1047 
1048 void test_rotateX() {
1049  Real cos1 = -0.852235910816473;
1050  Real sin1 = -0.523157674429820;
1051  Real chain1vals[] = { 0.0, 0.0, 0.0, 0.644662961550365, 1.376520855637543, 0.0, 2.163059129530747, 1.374750785448438, -0.069784983442129, 2.807722091081112, 0.0, 0.0 };
1052  utility::vector1<utility::vector1<Real> > chain1 (4), chain1a;
1053  // populate chain1
1054  int ind=0;
1055  for ( int i=1; i<=4; i++ ) {
1056  chain1[i].resize(3);
1057  for ( int j=1; j<=3; j++ ) {
1058  chain1[i][j]=chain1vals[ind++];
1059  }
1060  }
1061  rotateX(chain1, cos1, sin1, chain1a);
1062  printMatrix(chain1a);
1063 }
1064 
1065 void test_chainTORS () {
1066  // inputs
1067  int n=8;
1068  utility::vector1<utility::vector1<Real> > atoms (n);
1069  Real atomVals[] = {2.1636, 1.3766, 0, 1.6816, 2.3109, 1.0978, 1.1202, 3.6372, 0.6119, -0.3753, 3.6409, 0.3402, -1.0556, 2.2881, 0.4729, -0.6424, 1.2567, -0.5643, 0.0, 0.0, 0.0, 1.5200, 0.0, 0.0}; // for populating the atoms input
1070 
1071  //outputs
1072  utility::vector1<Real> t_ang, b_ang, b_len, R0;
1073  utility::vector1<utility::vector1<Real> > Q;
1074 
1075  // allocate and populate the atoms
1076  int ind=0;
1077  for ( int i=1; i<=n; i++ ) {
1078  atoms[i].resize(3);
1079  for ( int j=1; j<=3; j++ ) {
1080  atoms[i][j]=atomVals[ind++];
1081  }
1082  }
1083 
1084  chainTORS(n, atoms, t_ang, b_ang, b_len, R0, Q);
1085  printVector(t_ang);
1086  printVector(b_ang);
1087  printVector(b_len);
1088  printVector(R0);
1089  printMatrix(Q);
1090 
1091  // output should be:
1092  //
1093  // [t_ang, b_ang, b_len]
1094  // 236.3251 114.9908 1.5200
1095  // 90.4073 115.0000 1.5200
1096  // 354.6349 115.0000 1.5200
1097  // 293.3675 115.0000 1.5200
1098  // 119.1794 115.0000 1.5200
1099  // 263.4791 115.0000 1.5200
1100  // 24.1804 115.0000 1.5200
1101  // 52.8272 115.0594 1.5196
1102  //
1103  // R0: [2.1636, 1.3766, 0]'
1104  // Q:
1105  // -0.3172 -0.2596 -0.9121
1106  // 0.6147 0.6762 -0.4062
1107  // 0.7222 -0.6895 -0.0549
1108 }
1109 
1110 void test_chainXYZ() {
1111  utility::vector1<utility::vector1<Real> > atoms, Q(3);
1112  int n = 4;
1113  utility::vector1<Real> b_len (n), b_ang(n), t_ang(n), R0 (3);
1114 
1115  // setup bond lengths, bond angles, and torsions
1116  for ( int i=1; i<=n; i++ ) {
1117  b_len[i]=1.52;
1118  b_ang[i]=115.0;
1119  }
1120  t_ang[1] = 334.5902397604552;
1121  t_ang[2] = 52.8271503523038;
1122  t_ang[3] = 334.6192090171498;
1123  t_ang[4] = 180.0;
1124 
1125  // Q is the identity matrix and R0 is the origin
1126  for ( int i=1; i<=3; i++ ) {
1127  R0[i]=0.0;
1128  Q[i].resize(3);
1129  for ( int j=1; j<=3; j++ ) {
1130  (i==j) ? Q[i][j]=1.0 : Q[i][j]=0.0;
1131  }
1132  }
1133 
1134  chainXYZ(n, b_len, b_ang, t_ang, false, R0, Q, atoms);
1135  printMatrix(atoms);
1136 
1137  // output should be:
1138  // 0 1.520000000000000 2.163648967954409 1.681566958628807
1139  // 0 0 1.376581274771391 2.310875247341517
1140  // 0 0 0 1.097766691562338
1141 }
1142 
1143 /* // DJM: this was for benchmarking current version of chainXYZ against earlier (slower) version.
1144 // v2 is now current version. 6-8-08
1145 void test_chainXYZv2() {
1146 utility::vector1<utility::vector1<Real> > atoms_v1, atoms_v2, Q(3);
1147 int n = 12;
1148 utility::vector1<Real> b_len (n), b_ang(n), t_ang(n), R0 (3);
1149 utility::vector1<int> index(3);
1150 int nits = 1000000; // number of test iterations for time measurement
1151 
1152 // setup bond lengths, bond angles, and torsions
1153 for (int i=1; i<=n; i++) {
1154 b_len[i]=1.52;
1155 b_ang[i]=115.0;
1156 }
1157 t_ang[1] = 334.5902397604552;
1158 t_ang[2] = 52.8271503523038;
1159 t_ang[3] = 334.6192090171498;
1160 t_ang[4] = 180.0;
1161 t_ang[5] = 334.5902397604552;
1162 t_ang[6] = 52.8271503523038;
1163 t_ang[7] = 334.6192090171498;
1164 t_ang[8] = 180.0;
1165 t_ang[9] = 334.5902397604552;
1166 t_ang[10] = 52.8271503523038;
1167 t_ang[11] = 334.6192090171498;
1168 t_ang[12] = 180.0;
1169 
1170 // Q is the identity matrix and R0 is the origin
1171 for (int i=1; i<=3; i++) {
1172 R0[i]=0.0;
1173 Q[i].resize(3);
1174 for (int j=1; j<=3; j++) {
1175 (i==j) ? Q[i][j]=1.0 : Q[i][j]=0.0;
1176 }
1177 }
1178 
1179 // set up the index vector to select which atoms will establish the frame
1180 index[1] = 1;
1181 //index[2] = n;
1182 //index[3] = 2;
1183 index[2] = 2;
1184 index[3] = 3;
1185 
1186 std::time_t const start_time_v1( std::time( NULL ) );
1187 for (int d=1; d<=nits; d++) {
1188 chainXYZ(n, b_len, b_ang, t_ang, false, R0, Q, atoms_v1);
1189 }
1190 std::time_t const end_time_v1( std::time( NULL ) );
1191 std::cout << "start time: " << start_time_v1 << std::endl;
1192 std::cout << "end time: " << end_time_v1 << std::endl;
1193 std::cout << "start time minus end time v1: " << end_time_v1 - start_time_v1 << std::endl;
1194 
1195 //std::cout << "chainXYZ v1 says " << std::endl;
1196 //printMatrix(atoms_v1);
1197 std::time_t const start_time_v2( std::time( NULL ) );
1198 for (int d=1; d<=nits; d++) {
1199 chainXYZv2(n, b_len, b_ang, t_ang, false, R0, Q, atoms_v2);
1200 }
1201 std::time_t const end_time_v2( std::time( NULL ) );
1202 std::cout << "start time: " << start_time_v2 << std::endl;
1203 std::cout << "end time: " << end_time_v2 << std::endl;
1204 std::cout << "start time minus end time v2: " << end_time_v2 - start_time_v2 << std::endl;
1205 //std::cout << "chainXYZ v2 says " << std::endl;
1206 //printMatrix(atoms_v2);
1207 
1208 // output should be:
1209 // 0 1.520000000000000 2.163648967954409 1.681566958628807
1210 // 0 0 1.376581274771391 2.310875247341517
1211 // 0 0 0 1.097766691562338
1212 }
1213 */
1214 
1215 void test_chainParams() {
1216  int n=4;
1217  Real vbond, xi, eta, delta;
1218  Real atomVals[] = { -0.160000000000000, 1.225000000000001, -2.247000000000002, \
1219  -0.642082009325601, 1.113981114034975, -0.809755503317661, \
1220  0.000939232598462, -0.000706370692062, -0.000808351874906, \
1221  -0.481142776727140, -0.111725256657088, 1.436436144807435 };
1222  utility::vector1<Real> R0;
1223  utility::vector1<utility::vector1<Real> > atoms (4), Q;
1224  int ind=0;
1225  for ( int i=1; i<=4; i++ ) {
1226  atoms[i].resize(3);
1227  for ( int j=1; j<=3; j++ ) {
1228  atoms[i][j]=atomVals[ind++];
1229  }
1230  }
1231  //printMatrix(atoms);
1232  chainParams(n, atoms, vbond, xi, eta, delta, R0, Q);
1233  //std::cout << "vbond: " << vbond << std::endl;
1234  //std::cout << "xi: " << xi << std::endl;
1235  //std::cout << "eta: " << eta << std::endl;
1236  //std::cout << "delta: " << delta << std::endl;
1237  //std::cout << "R0: " << std::endl;
1238  printVector(R0);
1239  //std::cout << "Q: " << std::endl;
1240  printMatrix(Q);
1241  // output should be:
1242  // vbond: 3.93162
1243  // xi: 20.511
1244  // eta: 20.511
1245  // delta: 180
1246 
1247 }
1248 
1249 void test_torsion() {
1250  utility::vector1<Real> a (3), b(3), c(3), d(3);
1251  a[1]= 0.000939232598461515;
1252  a[2]= -0.000706370692062031;
1253  a[3]= -0.000808351874906243;
1254  b[1]= -0.481142776727140;
1255  b[2]= -0.111725256657088;
1256  b[3]= 1.436436144807435;
1257  c[1]= -0.160000000000000;
1258  c[2]= 1.225000000000001;
1259  c[3]= -2.247000000000002;
1260  d[1]= -0.642082009325601;
1261  d[2]= 1.113981114034975;
1262  d[3]= -0.809755503317661;
1263  Real tors=torsion(a,b,c,d);
1264  std::cout << tors << std::endl;
1265  // output should be 180.0
1266 }
1267 
1268 void test_bondangle() {
1269  utility::vector1<Real> a (3), b (3), c (3);
1270  a[1]= 0.000939232598461515;
1271  a[2]= -0.000706370692062031;
1272  a[3]= -0.000808351874906243;
1273  b[1]= -0.481142776727140;
1274  b[2]= -0.111725256657088;
1275  b[3]= 1.436436144807435;
1276  c[1]= -0.160000000000000;
1277  c[2]= 1.225000000000001;
1278  c[3]= -2.247000000000002;
1279  Real angle=bondangle(a,b,c);
1280  std::cout << angle << std::endl;
1281  // output should be 20.510953731276413
1282 }
1283 
1284 void test_scpn() {
1285  utility::vector1<Real> a (3), b (3), c (3);
1286  a[1]= 0.000939232598461515;
1287  a[2]= -0.000706370692062031;
1288  a[3]= -0.000808351874906243;
1289  b[1]= -0.481142776727140;
1290  b[2]= -0.111725256657088;
1291  b[3]= 1.436436144807435;
1292  c[1]= -0.160000000000000;
1293  c[2]= 1.225000000000001;
1294  c[3]= -2.247000000000002;
1295  Real scpn_val=scpn(a,b,c);
1296  std::cout << scpn_val << std::endl;
1297  // output should be 5.597217231925712
1298 }
1299 
1300 void test_eucDistance() {
1301  utility::vector1<Real> a (3), b (3);
1302  a[1]= -0.481142776727140;
1303  a[2]= -0.111725256657088;
1304  a[3]= 1.436436144807435;
1305  b[1]= -0.160000000000000;
1306  b[2]= 1.225000000000001;
1307  b[3]= -2.247000000000002;
1308  Real dist=eucDistance(a,b);
1309  std::cout << dist << std::endl;
1310  // output should be 3.931624209878514
1311 }
1312 
1313 void test_frame() {
1314  utility::vector1<utility::vector1<Real> > R (3), U;
1315  Real Rvals[] = { -0.160000000000000, 1.225000000000001, -2.247000000000002, \
1316  -0.481142776727140, -0.111725256657088, 1.436436144807435, \
1317  -0.642082009325601, 1.113981114034975, -0.809755503317661 };
1318 
1319  // intialize R
1320  int ind=0;
1321  for ( int i=1; i<=3; i++ ) {
1322  R[i].resize(3);
1323  for ( int j=1; j<=3; j++ ) {
1324  R[i][j]=Rvals[ind++];
1325  }
1326  }
1327  frame(R, U);
1328  printMatrix(U);
1329 }
1330 
1331 void test_cross() {
1332  utility::vector1<Real> L (3), r0 (3), r;
1333  L[1]= -0.081681961343163;
1334  L[2]= -0.339993139043773;
1335  L[3]= 0.936873909656095;
1336  r0[1]= -0.160939232598461;
1337  r0[2]= 1.225706370692063;
1338  r0[3]= -2.246191648125096;
1339  cross(L, r0, r);
1340  printVector(r);
1341 }
1342 
1343 void test_triaxialCoefficients() {
1344  // inputs
1345  utility::vector1<Real> vb (3), xi (3), eta (3), delta (3), theta (3);
1346  utility::vector1<int> order (3);
1347 
1348  // outputs
1349  utility::vector1<Real> cal, sal;
1350  utility::vector1<utility::vector1<Real> > A, B, C, D;
1351  int f;
1352 
1353  // dixon outputs (if testing dixon as well)
1354  utility::vector1<utility::vector1<Real> > sines, cosines, tau;
1355  int nsol;
1356 
1357  // initialize inputs
1358  vb[1] = 2.807722146231753;
1359  vb[2] = 2.563909995271172;
1360  vb[3] = 2.807373715195975;
1361  xi[1] = 1.132805948744302;
1362  xi[2] = 0.567232006898157;
1363  xi[3] = 1.133000954486070;
1364  eta[1] = 1.132805948744302;
1365  eta[2] = 0.567232006898157;
1366  eta[3] = 1.133000954486070;
1367  delta[1]= 6.232466452684100;
1368  delta[2]= 0.0;
1369  delta[3]= 6.235542719576548;
1370  theta[1]= 2.007128639793479;
1371  theta[2]= 2.007128639793479;
1372  theta[3]= 2.007128639793479;
1373  order[1]= 1;
1374  order[2]= 2;
1375  order[3]= 3;
1376 
1377  triaxialCoefficients(vb, xi, eta, delta, theta, order, A, B, C, D, cal, sal, f);
1378  printMatrix(B);
1379  dixon_sturm(A,B,C,D,order,sines,cosines,tau,nsol);
1380  //std::cout << nsol << std::endl;
1381  //printMatrix(sines);
1382 
1383  // cal and sal output should be
1384  //cal[1]= -0.583014272556966;
1385  //cal[2]= -0.456717749075032;
1386  //cal[3]= -0.456502620110421;
1387  //sal[1]= 0.812461911719480;
1388  //sal[2]= 0.889611655544056;
1389  //sal[3]= 0.889722067744934;
1390 }
1391 
1392 void test_sincos() {
1393  utility::vector1<Real> theta (3), cosine, sine;
1394  int flag=2;
1395  theta[1]=1.132805948744302;
1396  theta[2]=0.567232006898157;
1397  theta[3]=1.133000954486070;
1398 
1399  sincos(theta, flag, cosine, sine);
1400 }
1401 
1402 void test_triangle() {
1403  utility::vector1<Real> vbond (3), calpha, salpha;
1404  // initialize vbond
1405  vbond[1]=2.807722146231753;
1406  vbond[2]=2.563909995271172;
1407  vbond[3]=2.807373715195975;
1408 
1409  triangle(vbond, calpha, salpha);
1410  printVector(calpha);
1411  //std::cout << std::endl;
1412  printVector(salpha);
1413 }
1414 
1415 } // end namespace kinematic_closure
1416 } // end namespace numeric
1417 
1418 /*
1419 int main(int argc, char** argv)
1420 {
1421 //test_triangle();
1422 //test_sincos();
1423 //test_triaxialCoefficients();
1424 //test_cross();
1425 //test_frame();
1426 //test_eucDistance();
1427 //test_scpn();
1428 //test_bondangle();
1429 //test_torsion();
1430 //test_chainParams();
1431 //test_chainXYZ();
1432 //numeric::kinematic_closure::test_chainXYZv2();
1433 //test_chainTORS();
1434 //test_rotateX();
1435 //test_bridgeObjects();
1436 numeric::kinematic_closure::test_bridgeObjects_v2();
1437 }
1438 */
numeric::kinematic_closure::multTransMatrix
void multTransMatrix(const utility::vector1< utility::vector1< Real > > &A, const utility::vector1< utility::vector1< Real > > &B, utility::vector1< utility::vector1< Real > > &C)
Definition: kinematic_closure_helpers.cc:83
numeric::kinematic_closure::triaxialCoefficients
void triaxialCoefficients(const utility::vector1< Real > &vb, const utility::vector1< Real > &xi, const utility::vector1< Real > &eta, const utility::vector1< Real > &delta, const utility::vector1< Real > &theta, const utility::vector1< int > &order, utility::vector1< utility::vector1< Real > > &A, utility::vector1< utility::vector1< Real > > &B, utility::vector1< utility::vector1< Real > > &C, utility::vector1< utility::vector1< Real > > &D, utility::vector1< Real > &cal, utility::vector1< Real > &sal, int &f)
Definition: bridgeObjects.cc:86
numeric::kinematic_closure::test_frame
void test_frame()
Definition: bridgeObjects.cc:1313
numeric::kinematic_closure::test_eucDistance
void test_eucDistance()
Definition: bridgeObjects.cc:1300
numeric::kinematic_closure::eucDistance
Real eucDistance(const utility::vector1< Real > &a, const utility::vector1< Real > &b)
Definition: bridgeObjects.cc:269
numeric::kinematic_closure::test_torsion
void test_torsion()
Definition: bridgeObjects.cc:1249
amino_acids.size
dictionary size
Definition: amino_acids.py:44
numeric::kinematic_closure::chainParams
void chainParams(const int &n, const utility::vector1< utility::vector1< Real > > &atoms, Real &vbond, Real &xi, Real &eta, Real &delta, utility::vector1< Real > &R0, utility::vector1< utility::vector1< Real > > &Q)
Definition: bridgeObjects.cc:329
numeric::kinematic_closure::rotateZ
void rotateZ(const utility::vector1< utility::vector1< Real > > &R, const Real &c, const Real &s, utility::vector1< utility::vector1< Real > > &Rz)
Definition: bridgeObjects.cc:632
numeric::kinematic_closure::test_rotateX
void test_rotateX()
Definition: bridgeObjects.cc:1048
numeric::kinematic_closure::test_triangle
void test_triangle()
Definition: bridgeObjects.cc:1402
test.T007_TracerIO.M
tuple M
Definition: T007_TracerIO.py:38
test.T850_SubClassing.v
tuple v
Definition: T850_SubClassing.py:44
deg2rad
#define deg2rad
Definition: bridgeObjects.cc:29
numeric::kinematic_closure::bondangle
Real bondangle(const utility::vector1< Real > &a, const utility::vector1< Real > &b, const utility::vector1< Real > &c)
Definition: bridgeObjects.cc:281
numeric::kinematic_closure::triangle
void triangle(const utility::vector1< Real > &vbond, utility::vector1< Real > &calpha, utility::vector1< Real > &salpha)
Definition: bridgeObjects.cc:39
Equations.angle
def angle
Definition: Equations.py:36
ObjexxFCL::space
char space()
Space Character.
Definition: char.functions.hh:141
test.G200_Module_Region.loop
tuple loop
Definition: G200_Module_Region.py:48
numeric::kinematic_closure::rotateY
void rotateY(const utility::vector1< utility::vector1< Real > > &R, const Real &c, const Real &s, utility::vector1< utility::vector1< Real > > &Ry)
Definition: bridgeObjects.cc:621
numeric::kinematic_closure::test_chainTORS
void test_chainTORS()
Definition: bridgeObjects.cc:1065
numeric::kinematic_closure::test_sincos
void test_sincos()
Definition: bridgeObjects.cc:1392
rad2deg
#define rad2deg
Definition: bridgeObjects.cc:28
numeric::kinematic_closure::scpn
Real scpn(const utility::vector1< Real > &a, const utility::vector1< Real > &b, const utility::vector1< Real > &c)
Definition: bridgeObjects.cc:273
numeric::kinematic_closure::printMatrix
void printMatrix(const utility::vector1< utility::vector1< Real > > &M)
prints the matrix
Definition: kinematic_closure_helpers.cc:42
C
#define C(a, b)
Definition: functions.cc:27
numeric::kinematic_closure::test_bridgeObjects
void test_bridgeObjects()
Definition: bridgeObjects.cc:892
docking.p
tuple p
Definition: docking.py:9
numeric::kinematic_closure::chainXYZ
void chainXYZ(const int &n, const utility::vector1< Real > &b_len1, const utility::vector1< Real > &b_ang1, const utility::vector1< Real > &t_ang1, const bool space, const utility::vector1< Real > &R0, const utility::vector1< utility::vector1< Real > > &Q, utility::vector1< utility::vector1< Real > > &atoms)
Definition: bridgeObjects.cc:365
numeric::kinematic_closure::sincos
void sincos(const utility::vector1< Real > &theta, const int flag, utility::vector1< Real > &cosine, utility::vector1< Real > &sine)
Definition: bridgeObjects.cc:55
pidegs
#define pidegs
Definition: bridgeObjects.cc:31
z
def z
Definition: .pymolSurfaceMode.py:215
numeric::kinematic_closure::test_chainXYZ
void test_chainXYZ()
Definition: bridgeObjects.cc:1110
ObjexxFCL::abs
T abs(T const &x)
std::abs( x ) == | x |
Definition: Fmath.hh:295
dixon.hh
Header file for dixon code for ALC.
test.G201_Module_Structure.l1
tuple l1
Definition: G201_Module_Structure.py:27
numeric::kinematic_closure::test_bondangle
void test_bondangle()
Definition: bridgeObjects.cc:1268
ObjexxFCL::format::L
std::string L(int const w, bool const &t)
Logical Format.
Definition: format.cc:289
types.hh
rosetta project type declarations. Should be kept updated with core/types.hh. This exists because num...
numeric::kinematic_closure::dixon_eig
void dixon_eig(PseudoMatrix const &A, PseudoMatrix const &B, PseudoMatrix const &C, PseudoMatrix const &D, vector1< int > const &, PseudoMatrix &cos, PseudoMatrix &sin, PseudoMatrix &u, int &num_solutions)
Definition: dixon.cc:431
numeric::kinematic_closure::test_triaxialCoefficients
void test_triaxialCoefficients()
Definition: bridgeObjects.cc:1343
numeric::Real
double Real
Definition: types.hh:39
numeric::kinematic_closure::multMatrix
void multMatrix(const utility::vector1< utility::vector1< Real > > &A, const utility::vector1< utility::vector1< Real > > &B, utility::vector1< utility::vector1< Real > > &C)
Definition: kinematic_closure_helpers.cc:62
app.pyrosetta_toolkit.modules.SQLPDB.s
tuple s
Definition: SQLPDB.py:309
numeric::kinematic_closure::test_cross
void test_cross()
Definition: bridgeObjects.cc:1331
ObjexxFCL::pow
DimensionExpressionPow pow(Dimension const &dim1, Dimension const &dim2)
pow( Dimension, Dimension )
Definition: DimensionExpressions.hh:738
test.T850_SubClassing.a
tuple a
Definition: T850_SubClassing.py:43
numeric::kinematic_closure::test_scpn
void test_scpn()
Definition: bridgeObjects.cc:1284
numeric::kinematic_closure::test_chainParams
void test_chainParams()
Definition: bridgeObjects.cc:1215
utility::io::oc::cout
ocstream cout(std::cout)
Wrapper around std::cout.
Definition: ocstream.hh:287
vector1.hh
vector1: std::vector with 1-based indexing
numeric::kinematic_closure::printVector
void printVector(const utility::vector1< Real > &V)
helper functions ///
Definition: kinematic_closure_helpers.cc:30
coordlib.dist
def dist
Definition: coordlib.py:16
numeric::kinematic_closure::cross
void cross(const utility::vector1< Real > &L, const utility::vector1< Real > &r0, utility::vector1< Real > &r)
Definition: bridgeObjects.cc:230
numeric::kinematic_closure::torsion
Real torsion(const utility::vector1< Real > &a, const utility::vector1< Real > &b, const utility::vector1< Real > &c, const utility::vector1< Real > &d)
Definition: bridgeObjects.cc:287
numeric::kinematic_closure::chainTORS
void chainTORS(const int &n, const utility::vector1< utility::vector1< Real > > &atoms, utility::vector1< Real > &t_ang, utility::vector1< Real > &b_ang, utility::vector1< Real > &b_len, utility::vector1< Real > &R0, utility::vector1< utility::vector1< Real > > &Q)
Definition: bridgeObjects.cc:576
numeric::kinematic_closure::frame
void frame(const utility::vector1< utility::vector1< Real > > &R, utility::vector1< utility::vector1< Real > > &U)
Definition: bridgeObjects.cc:244
test.T850_SubClassing.b
tuple b
Definition: T850_SubClassing.py:144
kinematic_closure_helpers.hh
Header file for kinematic_closure_helpers.cc.
ObjexxFCL::format::A
std::string A(int const w, char const c)
char Format
Definition: format.cc:175
y
def y
Definition: .pymolSurfaceMode.py:203
numeric::kinematic_closure::bridgeObjects
void bridgeObjects(const utility::vector1< utility::vector1< Real > > &atoms, const utility::vector1< Real > &dt, const utility::vector1< Real > &da, const utility::vector1< Real > &db, const utility::vector1< int > &pivots, const utility::vector1< int > &order, utility::vector1< utility::vector1< Real > > &t_ang, utility::vector1< utility::vector1< Real > > &b_ang, utility::vector1< utility::vector1< Real > > &b_len, int &nsol)
Solve the triaxial loop closure problem for a system of atoms.
Definition: bridgeObjects.cc:680
demo.D060_Folding.f
tuple f
Definition: D060_Folding.py:445
numeric::kinematic_closure::dixon_sturm
void dixon_sturm(const utility::vector1< utility::vector1< Real > > &A, const utility::vector1< utility::vector1< Real > > &B, const utility::vector1< utility::vector1< Real > > &C, const utility::vector1< utility::vector1< Real > > &D, const utility::vector1< int > &order, utility::vector1< utility::vector1< Real > > &cos, utility::vector1< utility::vector1< Real > > &sin, utility::vector1< utility::vector1< Real > > &tau, int &nsol)
Definition: dixon.cc:492
numeric::kinematic_closure::rotateX
void rotateX(const utility::vector1< utility::vector1< Real > > &R, const Real &c, const Real &s, utility::vector1< utility::vector1< Real > > &Rx)
Definition: bridgeObjects.cc:610
Generated on Sun Jan 17 2016 15:21:19 for Rosetta by   doxygen 1.8.7