fitrms.C

Go to the documentation of this file.
00001 /*
00002 
00003 Code in this file was taken from PyMol v0.90 and user by permissing under
00004 the following license agreement contained in the PyMol distribution.  
00005 Trivial modifications have been made to permit incorporation into VMD.
00006 
00007 
00008 
00009 PyMOL Copyright Notice
00010 ======================
00011 
00012 The PyMOL source code is copyrighted, but you can freely use and copy
00013 it as long as you don't change or remove any of the copyright notices.
00014 
00015 ----------------------------------------------------------------------
00016 PyMOL is Copyright 1998-2003 by Warren L. DeLano of 
00017 DeLano Scientific LLC, San Carlos, CA, USA (www.delanoscientific.com).
00018 
00019                         All Rights Reserved
00020 
00021 Permission to use, copy, modify, distribute, and distribute modified 
00022 versions of this software and its documentation for any purpose and 
00023 without fee is hereby granted, provided that the above copyright 
00024 notice appear in all copies and that both the copyright notice and 
00025 this permission notice appear in supporting documentation, and that 
00026 the names of Warren L. DeLano or DeLano Scientific LLC not be used in 
00027 advertising or publicity pertaining to distribution of the software 
00028 without specific, written prior permission.
00029 
00030 WARREN LYFORD DELANO AND DELANO SCIENTIFIC LLC DISCLAIM ALL WARRANTIES 
00031 WITH REGARD TO THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF
00032 MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL WARREN LYFORD DELANO
00033 OR DELANO SCIENTIFIC LLC BE LIABLE FOR ANY SPECIAL, INDIRECT OR 
00034 CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
00035 OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE
00036 OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE 
00037 USE OR PERFORMANCE OF THIS SOFTWARE.
00038 ----------------------------------------------------------------------
00039 
00040 Where indicated, portions of the PyMOL system are instead protected
00041 under the copyrights of the respective authors.  However, all code in
00042 the PyMOL system is released as non-restrictive open-source software
00043 under the above license or an equivalent license.  
00044 
00045 PyMOL Trademark Notice
00046 ======================
00047 
00048 PyMOL(TM) is a trademark of DeLano Scientific LLC.  Derivate software
00049 which contains PyMOL source code must be plainly distinguished from
00050 the PyMOL package distributed by DeLano Scientific LLC in all publicity,
00051 advertising, and documentation.
00052 
00053 The slogans, "Includes PyMOL(TM).", "Based on PyMOL(TM) technology.",
00054 "Contains PyMOL(TM) source code.", and "Built using PyMOL(TM).", may
00055 be used in advertising, publicity, and documentation of derivate
00056 software provided that the notice, "PyMOL is a trademark of DeLano
00057 Scientific LLC.", is included in a footnote or at the end of the document.
00058 
00059 All other endorsements employing the PyMOL trademark require specific,
00060 written prior permission.
00061 
00062 --Warren L. DeLano (warren@delanoscientific.com)
00063 
00064 */
00065 
00066 #include <stdio.h>
00067 #include <math.h>
00068 
00069 #include "common.h"
00070 
00071 #ifdef R_SMALL
00072 #undef R_SMALL
00073 #endif
00074 #define R_SMALL 0.000000001
00075 
00076 static void normalize3d(BigReal *v) {
00077   BigReal vlen;
00078   vlen = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
00079   if (vlen > R_SMALL) {
00080     v[0] /= vlen;
00081     v[1] /= vlen;
00082     v[2] /= vlen;
00083   } else {
00084     v[0] = 0;
00085     v[1] = 0;
00086     v[2] = 0;
00087   }
00088 }
00089 
00090 /*========================================================================*/
00091 BigReal MatrixFitRMS(int n, BigReal *v1, BigReal *v2, const BigReal *wt, BigReal *ttt)
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 }
00310 
00311 

Generated on Thu Sep 21 01:17:12 2017 for NAMD by  doxygen 1.4.7