NAMD
Macros | Functions
fitrms.C File Reference
#include <stdio.h>
#include <math.h>
#include "common.h"

Go to the source code of this file.

Macros

#define R_SMALL   0.000000001
 

Functions

static void normalize3d (BigReal *v)
 
BigReal MatrixFitRMS (int n, BigReal *v1, BigReal *v2, const BigReal *wt, BigReal *ttt)
 

Macro Definition Documentation

#define R_SMALL   0.000000001

Definition at line 74 of file fitrms.C.

Referenced by normalize3d().

Function Documentation

BigReal MatrixFitRMS ( int  n,
BigReal v1,
BigReal v2,
const BigReal wt,
BigReal ttt 
)

Definition at line 91 of file fitrms.C.

References normalize3d(), x, and xx.

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 }
virial xx
static void normalize3d(BigReal *v)
Definition: fitrms.C:76
gridSize x
double BigReal
Definition: common.h:114
static void normalize3d ( BigReal v)
static

Definition at line 76 of file fitrms.C.

References R_SMALL.

Referenced by MatrixFitRMS().

76  {
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 }
#define R_SMALL
Definition: fitrms.C:74
double BigReal
Definition: common.h:114