fitrms.h File Reference

#include "common.h"

Go to the source code of this file.

Functions

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


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(), and x.

00092 {
00093   /*
00094         Subroutine to do the actual RMS fitting of two sets of vector coordinates
00095         This routine does not rotate the actual coordinates, but instead returns 
00096         the RMS fitting value, along with the center-of-mass translation vectors 
00097         T1 and T2 and the rotation vector M, which rotates the translated 
00098         coordinates of molecule 2 onto the translated coordinates of molecule 1.
00099   */
00100 
00101   BigReal *vv1,*vv2;
00102   BigReal m[3][3],aa[3][3],x[3],xx[3];
00103   BigReal sumwt, tol, sig, gam;
00104   BigReal sg, bb, cc, err, etmp, tmp;
00105   int a, b, c, maxiter, iters, ix, iy, iz;
00106   BigReal t1[3],t2[3];
00107 
00108   /* Initialize arrays. */
00109 
00110   for(a=0;a<3;a++) {
00111         for(b=0;b<3;b++) {
00112           m[a][b] = 0.0F;
00113           aa[a][b] = 0.0F;
00114         }
00115         m[a][a] = 1.0F;
00116         t1[a]=0.0F;
00117         t2[a]=0.0F;
00118   }
00119 
00120   sumwt = 0.0F;
00121 #if 0
00122   tol = SettingGet(cSetting_fit_tolerance);
00123 #else
00124   tol = 0.00001F;
00125 #endif
00126 #if 0
00127   maxiter = (int)SettingGet(cSetting_fit_iterations);
00128 #else
00129   maxiter = 10000;
00130 #endif
00131 
00132   /* Calculate center-of-mass vectors */
00133 
00134   vv1=v1;
00135   vv2=v2;
00136 
00137   if(wt) {
00138         for(c=0;c<n;c++)
00139           {
00140                 for(a=0;a<3;a++) {
00141                   t1[a] += wt[c]*vv1[a];
00142                   t2[a] += wt[c]*vv2[a];
00143                 }
00144                 if (wt[c]!=0.0F) {
00145                   sumwt = sumwt + wt[c];
00146                 } else {
00147                   sumwt = sumwt + 1.0F; /* WHAT IS THIS? */
00148                 }
00149                 vv1+=3;
00150                 vv2+=3;
00151           }
00152   } else {
00153         for(c=0;c<n;c++)
00154           {
00155                 for(a=0;a<3;a++) {
00156                   t1[a] += vv1[a];
00157                   t2[a] += vv2[a];
00158                 }
00159                 sumwt+=1.0F;
00160                 vv1+=3;
00161                 vv2+=3;
00162           }
00163   }
00164   if(sumwt==0.0F) sumwt = 1.0F;
00165   for(a=0;a<3;a++) {
00166         t1[a] /= sumwt;
00167         t2[a] /= sumwt;
00168   }
00169   /* Calculate correlation matrix */
00170   vv1=v1;
00171   vv2=v2;
00172   for(c=0;c<n;c++)
00173         {
00174           if(wt) {
00175                 for(a=0;a<3;a++) {
00176                   x[a] = wt[c]*(vv1[a] - t1[a]);
00177                   xx[a] = wt[c]*(vv2[a] - t2[a]);
00178                 }
00179           } else {
00180                 for(a=0;a<3;a++) {
00181                   x[a] = vv1[a] - t1[a];
00182                   xx[a] = vv2[a] - t2[a];
00183                 }
00184           }
00185           for(a=0;a<3;a++)
00186                 for(b=0;b<3;b++)
00187                   aa[a][b] = aa[a][b] + xx[a]*x[b];
00188           vv1+=3;
00189           vv2+=3;
00190         }
00191   if(n>1) {
00192     /* Primary iteration scheme to determine rotation matrix for molecule 2 */
00193     iters = 0;
00194     while(1) {
00195       /*        for(a=0;a<3;a++)
00196          {
00197          for(b=0;b<3;b++) 
00198          printf("%8.3f ",m[a][b]);
00199          printf("\n");
00200          }
00201          for(a=0;a<3;a++)
00202          {
00203          for(b=0;b<3;b++) 
00204          printf("%8.3f ",aa[a][b]);
00205          printf("\n");
00206          }
00207          printf("\n");
00208       */
00209       
00210       /* IX, IY, and IZ rotate 1-2-3, 2-3-1, 3-1-2, etc.*/
00211       iz = (iters+1) % 3;
00212       iy = (iz+1) % 3;
00213       ix = (iy+1) % 3;
00214       sig = aa[iz][iy] - aa[iy][iz];
00215       gam = aa[iy][iy] + aa[iz][iz];
00216 
00217       if(iters>=maxiter) 
00218         {
00219 #if 0
00220           PRINTFB(FB_Matrix,FB_Details)
00221 #else
00222             fprintf(stderr,
00223 #endif
00224             " Matrix: Warning: no convergence (%1.8f<%1.8f after %d iterations).\n",(BigReal)tol,(BigReal)gam,iters
00225 #if 0
00226             ENDFB;
00227 #else
00228                 );
00229 #endif
00230           break;
00231         }
00232 
00233       /* Determine size of off-diagonal element.  If off-diagonals exceed the
00234          diagonal elements * tolerance, perform Jacobi rotation. */
00235       tmp = sig*sig + gam*gam;
00236       sg = sqrt(tmp);
00237       if((sg!=0.0F) &&(fabs(sig)>(tol*fabs(gam)))) {
00238         sg = 1.0F / sg;
00239         for(a=0;a<3;a++)
00240           {
00241             bb = gam*aa[iy][a] + sig*aa[iz][a];
00242             cc = gam*aa[iz][a] - sig*aa[iy][a];
00243             aa[iy][a] = bb*sg;
00244             aa[iz][a] = cc*sg;
00245             
00246             bb = gam*m[iy][a] + sig*m[iz][a];
00247             cc = gam*m[iz][a] - sig*m[iy][a];
00248             m[iy][a] = bb*sg;
00249             m[iz][a] = cc*sg;
00250           }
00251       } else {
00252         break;
00253       }
00254       iters++;
00255     }
00256   }
00257   /* At this point, we should have a converged rotation matrix (M).  Calculate
00258          the weighted RMS error. */
00259   err = 0.0F;
00260   vv1=v1;
00261   vv2=v2;
00262 
00263   normalize3d(m[0]);
00264   normalize3d(m[1]);
00265   normalize3d(m[2]);
00266   for(c=0;c<n;c++) {
00267         etmp = 0.0F;
00268         for(a=0;a<3;a++) {
00269           tmp = m[a][0]*(vv2[0]-t2[0])
00270                 + m[a][1]*(vv2[1]-t2[1])
00271                 + m[a][2]*(vv2[2]-t2[2]);
00272           tmp = (vv1[a]-t1[a])-tmp;
00273           etmp += tmp*tmp;
00274         }
00275         if(wt)
00276           err += wt[c] * etmp;
00277         else 
00278           err += etmp;
00279         vv1+=3;
00280         vv2+=3;
00281   }
00282 
00283   err=err/sumwt;
00284   err=sqrt(err);
00285 
00286   ttt[0]=(BigReal)m[0][0];
00287   ttt[1]=(BigReal)m[0][1];
00288   ttt[2]=(BigReal)m[0][2];
00289   ttt[3]=(BigReal)-t1[0];
00290   ttt[4]=(BigReal)m[1][0];
00291   ttt[5]=(BigReal)m[1][1];
00292   ttt[6]=(BigReal)m[1][2];
00293   ttt[7]=(BigReal)-t1[1];
00294   ttt[8]=(BigReal)m[2][0];
00295   ttt[9]=(BigReal)m[2][1];
00296   ttt[10]=(BigReal)m[2][2];
00297   ttt[11]=(BigReal)-t1[2];
00298   ttt[12]=(BigReal)t2[0];
00299   ttt[13]=(BigReal)t2[1];
00300   ttt[14]=(BigReal)t2[2];
00301   ttt[15]=1.0F; /* for compatibility with normal 4x4 matrices */
00302 
00303 #if 0
00304   if(fabs(err)<R_SMALL4)
00305     err=0.0F;
00306 #endif
00307 
00308   return((BigReal)err);
00309 }


Generated on Sat Sep 23 01:17:16 2017 for NAMD by  doxygen 1.4.7