NAMD
fitrms.C
Go to the documentation of this file.
1 /*
2 
3 Code in this file was taken from PyMol v0.90 and user by permissing under
4 the following license agreement contained in the PyMol distribution.
5 Trivial modifications have been made to permit incorporation into VMD.
6 
7 
8 
9 PyMOL Copyright Notice
10 ======================
11 
12 The PyMOL source code is copyrighted, but you can freely use and copy
13 it as long as you don't change or remove any of the copyright notices.
14 
15 ----------------------------------------------------------------------
16 PyMOL is Copyright 1998-2003 by Warren L. DeLano of
17 DeLano Scientific LLC, San Carlos, CA, USA (www.delanoscientific.com).
18 
19  All Rights Reserved
20 
21 Permission to use, copy, modify, distribute, and distribute modified
22 versions of this software and its documentation for any purpose and
23 without fee is hereby granted, provided that the above copyright
24 notice appear in all copies and that both the copyright notice and
25 this permission notice appear in supporting documentation, and that
26 the names of Warren L. DeLano or DeLano Scientific LLC not be used in
27 advertising or publicity pertaining to distribution of the software
28 without specific, written prior permission.
29 
30 WARREN LYFORD DELANO AND DELANO SCIENTIFIC LLC DISCLAIM ALL WARRANTIES
31 WITH REGARD TO THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF
32 MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL WARREN LYFORD DELANO
33 OR DELANO SCIENTIFIC LLC BE LIABLE FOR ANY SPECIAL, INDIRECT OR
34 CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
35 OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE
36 OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE
37 USE OR PERFORMANCE OF THIS SOFTWARE.
38 ----------------------------------------------------------------------
39 
40 Where indicated, portions of the PyMOL system are instead protected
41 under the copyrights of the respective authors. However, all code in
42 the PyMOL system is released as non-restrictive open-source software
43 under the above license or an equivalent license.
44 
45 PyMOL Trademark Notice
46 ======================
47 
48 PyMOL(TM) is a trademark of DeLano Scientific LLC. Derivate software
49 which contains PyMOL source code must be plainly distinguished from
50 the PyMOL package distributed by DeLano Scientific LLC in all publicity,
51 advertising, and documentation.
52 
53 The slogans, "Includes PyMOL(TM).", "Based on PyMOL(TM) technology.",
54 "Contains PyMOL(TM) source code.", and "Built using PyMOL(TM).", may
55 be used in advertising, publicity, and documentation of derivate
56 software provided that the notice, "PyMOL is a trademark of DeLano
57 Scientific LLC.", is included in a footnote or at the end of the document.
58 
59 All other endorsements employing the PyMOL trademark require specific,
60 written prior permission.
61 
62 --Warren L. DeLano (warren@delanoscientific.com)
63 
64 */
65 
66 #include <stdio.h>
67 #include <math.h>
68 
69 #include "common.h"
70 
71 #ifdef R_SMALL
72 #undef R_SMALL
73 #endif
74 #define R_SMALL 0.000000001
75 
76 static void normalize3d(BigReal *v) {
77  BigReal vlen;
78  vlen = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
79  if (vlen > R_SMALL) {
80  v[0] /= vlen;
81  v[1] /= vlen;
82  v[2] /= vlen;
83  } else {
84  v[0] = 0;
85  v[1] = 0;
86  v[2] = 0;
87  }
88 }
89 
90 /*========================================================================*/
91 BigReal MatrixFitRMS(int n, BigReal *v1, BigReal *v2, const BigReal *wt, BigReal *ttt)
92 {
93  /*
94  Subroutine to do the actual RMS fitting of two sets of vector coordinates
95  This routine does not rotate the actual coordinates, but instead returns
96  the RMS fitting value, along with the center-of-mass translation vectors
97  T1 and T2 and the rotation vector M, which rotates the translated
98  coordinates of molecule 2 onto the translated coordinates of molecule 1.
99  */
100 
101  BigReal *vv1,*vv2;
102  BigReal m[3][3],aa[3][3],x[3],xx[3];
103  BigReal sumwt, tol, sig, gam;
104  BigReal sg, bb, cc, err, etmp, tmp;
105  int a, b, c, maxiter, iters, ix, iy, iz;
106  BigReal t1[3],t2[3];
107 
108  /* Initialize arrays. */
109 
110  for(a=0;a<3;a++) {
111  for(b=0;b<3;b++) {
112  m[a][b] = 0.0F;
113  aa[a][b] = 0.0F;
114  }
115  m[a][a] = 1.0F;
116  t1[a]=0.0F;
117  t2[a]=0.0F;
118  }
119 
120  sumwt = 0.0F;
121 #if 0
122  tol = SettingGet(cSetting_fit_tolerance);
123 #else
124  tol = 0.00001F;
125 #endif
126 #if 0
127  maxiter = (int)SettingGet(cSetting_fit_iterations);
128 #else
129  maxiter = 10000;
130 #endif
131 
132  /* Calculate center-of-mass vectors */
133 
134  vv1=v1;
135  vv2=v2;
136 
137  if(wt) {
138  for(c=0;c<n;c++)
139  {
140  for(a=0;a<3;a++) {
141  t1[a] += wt[c]*vv1[a];
142  t2[a] += wt[c]*vv2[a];
143  }
144  if (wt[c]!=0.0F) {
145  sumwt = sumwt + wt[c];
146  } else {
147  sumwt = sumwt + 1.0F; /* WHAT IS THIS? */
148  }
149  vv1+=3;
150  vv2+=3;
151  }
152  } else {
153  for(c=0;c<n;c++)
154  {
155  for(a=0;a<3;a++) {
156  t1[a] += vv1[a];
157  t2[a] += vv2[a];
158  }
159  sumwt+=1.0F;
160  vv1+=3;
161  vv2+=3;
162  }
163  }
164  if(sumwt==0.0F) sumwt = 1.0F;
165  for(a=0;a<3;a++) {
166  t1[a] /= sumwt;
167  t2[a] /= sumwt;
168  }
169  /* Calculate correlation matrix */
170  vv1=v1;
171  vv2=v2;
172  for(c=0;c<n;c++)
173  {
174  if(wt) {
175  for(a=0;a<3;a++) {
176  x[a] = wt[c]*(vv1[a] - t1[a]);
177  xx[a] = wt[c]*(vv2[a] - t2[a]);
178  }
179  } else {
180  for(a=0;a<3;a++) {
181  x[a] = vv1[a] - t1[a];
182  xx[a] = vv2[a] - t2[a];
183  }
184  }
185  for(a=0;a<3;a++)
186  for(b=0;b<3;b++)
187  aa[a][b] = aa[a][b] + xx[a]*x[b];
188  vv1+=3;
189  vv2+=3;
190  }
191  if(n>1) {
192  /* Primary iteration scheme to determine rotation matrix for molecule 2 */
193  iters = 0;
194  while(1) {
195  /* for(a=0;a<3;a++)
196  {
197  for(b=0;b<3;b++)
198  printf("%8.3f ",m[a][b]);
199  printf("\n");
200  }
201  for(a=0;a<3;a++)
202  {
203  for(b=0;b<3;b++)
204  printf("%8.3f ",aa[a][b]);
205  printf("\n");
206  }
207  printf("\n");
208  */
209 
210  /* IX, IY, and IZ rotate 1-2-3, 2-3-1, 3-1-2, etc.*/
211  iz = (iters+1) % 3;
212  iy = (iz+1) % 3;
213  ix = (iy+1) % 3;
214  sig = aa[iz][iy] - aa[iy][iz];
215  gam = aa[iy][iy] + aa[iz][iz];
216 
217  if(iters>=maxiter)
218  {
219 #if 0
220  PRINTFB(FB_Matrix,FB_Details)
221 #else
222  fprintf(stderr,
223 #endif
224  " Matrix: Warning: no convergence (%1.8f<%1.8f after %d iterations).\n",(BigReal)tol,(BigReal)gam,iters
225 #if 0
226  ENDFB;
227 #else
228  );
229 #endif
230  break;
231  }
232 
233  /* Determine size of off-diagonal element. If off-diagonals exceed the
234  diagonal elements * tolerance, perform Jacobi rotation. */
235  tmp = sig*sig + gam*gam;
236  sg = sqrt(tmp);
237  if((sg!=0.0F) &&(fabs(sig)>(tol*fabs(gam)))) {
238  sg = 1.0F / sg;
239  for(a=0;a<3;a++)
240  {
241  bb = gam*aa[iy][a] + sig*aa[iz][a];
242  cc = gam*aa[iz][a] - sig*aa[iy][a];
243  aa[iy][a] = bb*sg;
244  aa[iz][a] = cc*sg;
245 
246  bb = gam*m[iy][a] + sig*m[iz][a];
247  cc = gam*m[iz][a] - sig*m[iy][a];
248  m[iy][a] = bb*sg;
249  m[iz][a] = cc*sg;
250  }
251  } else {
252  break;
253  }
254  iters++;
255  }
256  }
257  /* At this point, we should have a converged rotation matrix (M). Calculate
258  the weighted RMS error. */
259  err = 0.0F;
260  vv1=v1;
261  vv2=v2;
262 
263  normalize3d(m[0]);
264  normalize3d(m[1]);
265  normalize3d(m[2]);
266  for(c=0;c<n;c++) {
267  etmp = 0.0F;
268  for(a=0;a<3;a++) {
269  tmp = m[a][0]*(vv2[0]-t2[0])
270  + m[a][1]*(vv2[1]-t2[1])
271  + m[a][2]*(vv2[2]-t2[2]);
272  tmp = (vv1[a]-t1[a])-tmp;
273  etmp += tmp*tmp;
274  }
275  if(wt)
276  err += wt[c] * etmp;
277  else
278  err += etmp;
279  vv1+=3;
280  vv2+=3;
281  }
282 
283  err=err/sumwt;
284  err=sqrt(err);
285 
286  ttt[0]=(BigReal)m[0][0];
287  ttt[1]=(BigReal)m[0][1];
288  ttt[2]=(BigReal)m[0][2];
289  ttt[3]=(BigReal)-t1[0];
290  ttt[4]=(BigReal)m[1][0];
291  ttt[5]=(BigReal)m[1][1];
292  ttt[6]=(BigReal)m[1][2];
293  ttt[7]=(BigReal)-t1[1];
294  ttt[8]=(BigReal)m[2][0];
295  ttt[9]=(BigReal)m[2][1];
296  ttt[10]=(BigReal)m[2][2];
297  ttt[11]=(BigReal)-t1[2];
298  ttt[12]=(BigReal)t2[0];
299  ttt[13]=(BigReal)t2[1];
300  ttt[14]=(BigReal)t2[2];
301  ttt[15]=1.0F; /* for compatibility with normal 4x4 matrices */
302 
303 #if 0
304  if(fabs(err)<R_SMALL4)
305  err=0.0F;
306 #endif
307 
308  return((BigReal)err);
309 }
310 
311 
#define R_SMALL
Definition: fitrms.C:74
virial xx
static void normalize3d(BigReal *v)
Definition: fitrms.C:76
BigReal MatrixFitRMS(int n, BigReal *v1, BigReal *v2, const BigReal *wt, BigReal *ttt)
Definition: fitrms.C:91
gridSize x
double BigReal
Definition: common.h:114