#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) |
|
||||||||||||||||||||||||
|
Definition at line 91 of file fitrms.C. References BigReal, and normalize3d(). 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 }
|
1.3.9.1