Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

Measure.C

Go to the documentation of this file.
00001 /***************************************************************************
00002  *cr
00003  *cr            (C) Copyright 1995-2019 The Board of Trustees of the
00004  *cr                        University of Illinois
00005  *cr                         All Rights Reserved
00006  *cr
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  * RCS INFORMATION:
00011  *
00012  *      $RCSfile: Measure.C,v $
00013  *      $Author: dgomes $        $Locker:  $             $State: Exp $
00014  *      $Revision: 1.156 $       $Date: 2024/05/16 14:56:56 $
00015  *
00016  ***************************************************************************
00017  * DESCRIPTION:
00018  *   Code to measure atom distances, angles, dihedrals, etc.
00019  ***************************************************************************/
00020 
00021 #include <stdio.h>
00022 #include <stdlib.h>
00023 #include <math.h>
00024 #include "Measure.h"
00025 #include "AtomSel.h"
00026 #include "Matrix4.h"
00027 #include "utilities.h"
00028 #include "fitrms.h"
00029 #include "ResizeArray.h"
00030 #include "SpatialSearch.h"  // for find_within()
00031 #include "MoleculeList.h"
00032 #include "Inform.h"
00033 #include "Timestep.h"
00034 #include "VMDApp.h"
00035 #include "WKFThreads.h"
00036 #include "WKFUtils.h"
00037 
00038 // Standard functions available to everyone
00039 static const char *measure_error_messages[] = {
00040   "no atom selection",                                  // -1
00041   "no atoms in selection",                              // -2
00042   "incorrect number of weights for selection",          // -3
00043   "internal error: NULL pointer given",                 // -4
00044   "bad weight sum, would cause divide by zero",         // -5
00045   "molecule was deleted(?)",                            // -6
00046   "cannot understand weight parameter",                 // -7
00047   "non-number given as parameter",                      // -8
00048   "two selections don't have the same number of atoms", // -9
00049   "internal error: out of range",                       // -10
00050   "no coordinates in selection",                        // -11
00051   "couldn't compute eigenvalue/vectors",                // -12
00052   "unknown Tcl problem",                                // -13
00053   "no atom radii",                                      // -14
00054   "order parameter contains out-of-range atom index",   // -15
00055   "order parameter not supported in old RMS fit mode",  // -16
00056   "specified frames are out of range",                  // -17
00057   "both selections must reference the same molecule",   // -18
00058   "one atom appears twice in list",                     // -19
00059   "molecule contains no frames",                        // -20
00060   "invalid atom id",                                    // -21
00061   "cutoff must be smaller than cell dimension",         // -22
00062   "Zero volmap gridsize"                                // -23
00063 };
00064   
00065 const char *measure_error(int errnum) {
00066   if (errnum >= 0 || errnum < -23) 
00067     return "bad error number";
00068   return measure_error_messages[-errnum - 1];
00069 }
00070 
00071 int measure_move(const AtomSel *sel, float *framepos, const Matrix4 &mat) {
00072   if (!sel)       return MEASURE_ERR_NOSEL;
00073   if (!framepos)  return MEASURE_ERR_NOFRAMEPOS;
00074 
00075   // and apply it to the coordinates
00076   int i;
00077   float *pos = framepos;
00078   if (mat.mat[3]==0 && mat.mat[7]==0 && mat.mat[11]==0 && mat.mat[15]==1) {
00079     // then this is just a rotation followed by a translation.  Do it
00080     // the fast way.
00081     const float dx = mat.mat[12];
00082     const float dy = mat.mat[13];
00083     const float dz = mat.mat[14];
00084     const float *m = mat.mat;
00085     if (sel->selected == sel->num_atoms) {
00086       for (i=0; i<sel->num_atoms; i++) {
00087         float x = pos[0]*m[0] + pos[1]*m[4] + pos[2]*m[8] + dx;
00088         float y = pos[0]*m[1] + pos[1]*m[5] + pos[2]*m[9] + dy;
00089         float z = pos[0]*m[2] + pos[1]*m[6] + pos[2]*m[10] + dz;
00090         pos[0] = x;
00091         pos[1] = y;
00092         pos[2] = z;
00093         pos += 3;
00094       }
00095     } else {
00096       pos += sel->firstsel*3L;
00097       for (i=sel->firstsel; i<=sel->lastsel; i++) {
00098         if (sel->on[i]) {
00099           float x = pos[0]*m[0] + pos[1]*m[4] + pos[2]*m[8] + dx;
00100           float y = pos[0]*m[1] + pos[1]*m[5] + pos[2]*m[9] + dy;
00101           float z = pos[0]*m[2] + pos[1]*m[6] + pos[2]*m[10] + dz;
00102           pos[0] = x;
00103           pos[1] = y;
00104           pos[2] = z;
00105         }
00106         pos += 3;
00107       }
00108     }
00109   } else {
00110     pos += sel->firstsel*3L;
00111     for (i=sel->firstsel; i<=sel->lastsel; i++) {
00112       if (sel->on[i]) {
00113         mat.multpoint3d(pos, pos);
00114       }
00115       pos += 3;
00116     }
00117   }
00118   return MEASURE_NOERR;
00119 }
00120 
00121 // compute the sum of a set of weights which correspond either to
00122 // all atoms, or to selected atoms
00123 int measure_sumweights(const AtomSel *sel, int numweights, 
00124                        const float *weights, float *weightsum) {
00125   if (!sel)       return MEASURE_ERR_NOSEL;
00126   if (!weights)   return MEASURE_ERR_NOWEIGHT;
00127   if (!weightsum) return MEASURE_ERR_GENERAL;
00128 
00129   double sum = 0;
00130 
00131   if (numweights == sel->num_atoms) {
00132     int i;
00133     for (i=sel->firstsel; i<=sel->lastsel; i++) {
00134       if (sel->on[i]) {
00135         sum += weights[i];
00136       }
00137     }
00138   } else if (numweights == sel->selected) {
00139     int i, j;
00140     for (j=0,i=sel->firstsel; i<=sel->lastsel; i++) {
00141       if (sel->on[i]) {
00142         sum += weights[j];
00143         j++;
00144       }
00145     }
00146   } else {
00147     return MEASURE_ERR_BADWEIGHTPARM;  
00148   }
00149 
00150   *weightsum = (float) sum;
00151   return MEASURE_NOERR;
00152 }
00153 
00154 extern int measure_center_perresidue(MoleculeList *mlist, const AtomSel *sel, const float *framepos,
00155                    const float *weight, float *com) {
00156   if (!sel)      return MEASURE_ERR_NOSEL;
00157   if (!framepos) return MEASURE_ERR_NOFRAMEPOS;
00158   if (!weight)   return MEASURE_ERR_NOWEIGHT;
00159   if (!com)      return MEASURE_ERR_NOCOM;
00160 
00161   Molecule *mol = mlist->mol_from_id(sel->molid());
00162   int residue = mol->atom(sel->firstsel)->uniq_resid;
00163   int rescount = 0;
00164   int i, j = 0;
00165   float x=0, y=0, z=0, w=0;
00166   for (i=sel->firstsel; i<=sel->lastsel; i++) {
00167     if (sel->on[i]) {
00168       if (residue != mol->atom(i)->uniq_resid) {
00169         com[3*rescount    ] = x / w;
00170         com[3*rescount + 1] = y / w;
00171         com[3*rescount + 2] = z / w;
00172         residue = mol->atom(i)->uniq_resid;
00173         rescount++;
00174         x = 0;
00175         y = 0;
00176         z = 0;
00177         w = 0;
00178       }
00179       float tw = weight[j];
00180       w += tw;
00181       x += tw * framepos[3*i  ];
00182       y += tw * framepos[3*i+1];
00183       z += tw * framepos[3*i+2];
00184       j++;
00185     }
00186   }
00187   com[3*rescount    ] = x / w;
00188   com[3*rescount + 1] = y / w;
00189   com[3*rescount + 2] = z / w;
00190   rescount++;
00191   return rescount;
00192 }
00193 
00194 // compute the center of mass
00195 // return -5 if total weight == 0, otherwise 0 for success.
00196 int measure_center(const AtomSel *sel, const float *framepos,
00197                    const float *weight, float *com) {
00198   if (!sel)      return MEASURE_ERR_NOSEL;
00199   if (!framepos) return MEASURE_ERR_NOFRAMEPOS;
00200   if (!weight)   return MEASURE_ERR_NOWEIGHT;
00201   if (!com)      return MEASURE_ERR_NOCOM;   
00202 
00203 #if 0 && (__cplusplus >= 201103L)
00204   double x, y, z, w;
00205   w=x=y=z=0;
00206 
00207   int j=0;
00208 
00209 #if 1
00210   sel->for_selected_lambda([&] (int i) {
00211                              float tw = weight[j];
00212                              w += double(tw);
00213 
00214                              long ind = 3L * i;
00215                              x += double(tw * framepos[ind  ]);
00216                              y += double(tw * framepos[ind+1]);
00217                              z += double(tw * framepos[ind+2]);
00218                              j++; 
00219                            });
00220 #else
00221   auto lambda_com = [&] (int i) {
00222     float tw = weight[j];
00223     w += double(tw);
00224 
00225     long ind = 3L * i;
00226     x += double(tw * framepos[ind  ]);
00227     y += double(tw * framepos[ind+1]);
00228     z += double(tw * framepos[ind+2]);
00229     j++; 
00230   };
00231 
00232   sel->for_selected_lambda(lambda_com);
00233 #endif
00234 
00235   if (w == 0) {
00236     return MEASURE_ERR_BADWEIGHTSUM;
00237   }
00238 
00239   x/=w;
00240   y/=w;
00241   z/=w;
00242 
00243   com[0] = float(x);
00244   com[1] = float(y);
00245   com[2] = float(z);
00246 #else
00247   // use double precision floating point in case we have a large selection
00248   int i, j;
00249   double x, y, z, w;
00250   j=0;
00251   w=x=y=z=0;
00252   for (i=sel->firstsel; i<=sel->lastsel; i++) {
00253     if (sel->on[i]) {
00254       float tw = weight[j];
00255       w += double(tw);
00256       x += double(tw * framepos[3L*i  ]);
00257       y += double(tw * framepos[3L*i+1]);
00258       z += double(tw * framepos[3L*i+2]);
00259       j++;
00260     }
00261   }
00262 
00263   if (w == 0) {
00264     return MEASURE_ERR_BADWEIGHTSUM;
00265   }
00266 
00267   x/=w;
00268   y/=w;
00269   z/=w;
00270 
00271   com[0] = float(x);
00272   com[1] = float(y);
00273   com[2] = float(z);
00274 #endif
00275 
00276   return MEASURE_NOERR;
00277 }
00278 
00279 
00280 // compute the axis aligned aligned bounding box for the selected atoms
00281 int measure_minmax(int num, const int *on, const float *framepos, 
00282                    const float *radii, float *min_coord, float *max_coord) {
00283   int i;
00284   float minx, miny, minz;
00285   float maxx, maxy, maxz;
00286 
00287   if (!on)                      return MEASURE_ERR_NOSEL;
00288   if (num == 0)                 return MEASURE_ERR_NOATOMS;
00289   if (!min_coord || !max_coord) return MEASURE_ERR_NOMINMAXCOORDS;
00290 
00291   vec_zero(min_coord);
00292   vec_zero(max_coord);
00293 
00294   // find first selected atom
00295   int firstsel = 0;
00296   int lastsel = -1;
00297   if (analyze_selection_aligned_dispatch(NULL, num, on, &firstsel, &lastsel, NULL))
00298     return MEASURE_NOERR; // no atoms selected
00299 
00300   // initialize minmax coords to the first selected atom
00301   i=firstsel;
00302   if (radii == NULL) {
00303     // calculate bounding box of atom centers
00304     minx = maxx = framepos[i*3L  ];
00305     miny = maxy = framepos[i*3L+1];
00306     minz = maxz = framepos[i*3L+2];
00307   } else {
00308     // calculate bounding box for atoms /w given radii 
00309     minx = framepos[i*3L  ] - radii[i];
00310     maxx = framepos[i*3L  ] + radii[i];
00311     miny = framepos[i*3L+1] - radii[i];
00312     maxy = framepos[i*3L+1] + radii[i];
00313     minz = framepos[i*3L+2] - radii[i];
00314     maxz = framepos[i*3L+2] + radii[i];
00315   }
00316 
00317   // continue looping from there until we finish
00318   if (radii == NULL) {
00319     // calculate bounding box of atom centers
00320 #if 0
00321     minmax_selected_3fv_aligned(framepos, on, atomSel->num_atoms,
00322                                 firstsel, lastsel, fmin, fmax);
00323 #else
00324     for (i++; i<=lastsel; i++) {
00325       if (on[i]) {
00326         long ind = i * 3L;
00327         float tmpx = framepos[ind  ];
00328         if (tmpx < minx) minx = tmpx; 
00329         if (tmpx > maxx) maxx = tmpx;
00330   
00331         float tmpy = framepos[ind+1];
00332         if (tmpy < miny) miny = tmpy; 
00333         if (tmpy > maxy) maxy = tmpy;
00334   
00335         float tmpz = framepos[ind+2];
00336         if (tmpz < minz) minz = tmpz; 
00337         if (tmpz > maxz) maxz = tmpz;
00338       }
00339     }
00340 #endif
00341   } else {
00342     // calculate bounding box for atoms /w given radii 
00343     for (i++; i<=lastsel; i++) {
00344       if (on[i]) {
00345         long ind = i * 3L;
00346         float mintmpx = framepos[ind  ] - radii[i];
00347         float maxtmpx = framepos[ind  ] + radii[i];
00348         if (mintmpx < minx) minx = mintmpx;
00349         if (maxtmpx > maxx) maxx = maxtmpx;
00350   
00351         float mintmpy = framepos[ind+1] - radii[i];
00352         float maxtmpy = framepos[ind+1] + radii[i];
00353         if (mintmpy < miny) miny = mintmpy; 
00354         if (maxtmpy > maxy) maxy = maxtmpy;
00355   
00356         float mintmpz = framepos[ind+2] - radii[i];
00357         float maxtmpz = framepos[ind+2] + radii[i];
00358         if (mintmpz < minz) minz = mintmpz; 
00359         if (maxtmpz > maxz) maxz = maxtmpz;
00360       }
00361     }
00362   }
00363 
00364   // set the final min/max output values
00365   min_coord[0]=minx;
00366   min_coord[1]=miny;
00367   min_coord[2]=minz;
00368   max_coord[0]=maxx;
00369   max_coord[1]=maxy;
00370   max_coord[2]=maxz;
00371 
00372   return MEASURE_NOERR;
00373 }
00374 
00375 /*I'm going to assume that *avpos points to an array the length of 3*sel. The return
00376 value will indicate the ACTUAL number of residues in the selection,
00377 which tells the caller how many values should be in the returned list, or a number
00378 less than zero on error.
00379 */
00380 int measure_avpos_perresidue(const AtomSel *sel, MoleculeList *mlist, 
00381                          int start, int end, int step, float *avpos)  {
00382   //Get the per-atom average position. We'll be accumulating on this array.
00383   int retval = measure_avpos(sel, mlist, start, end, step, avpos);
00384   if (retval != MEASURE_NOERR) {
00385     return retval;
00386   }
00387   Molecule *mol = mlist->mol_from_id(sel->molid());
00388   int residue = mol->atom(sel->firstsel)->uniq_resid;
00389   int rescount = 0;
00390   int ressize = 0;
00391   float accumulate[3] = {0.0,0.0,0.0};
00392   int j = 0, k, i;
00393   for (i = sel->firstsel; i <= sel->lastsel; i++) {
00394     if (sel->on[i]) {
00395       if (residue != mol->atom(i)->uniq_resid) {
00396         for (k = 0; k < 3; k++) {
00397           avpos[3*rescount + k] = accumulate[k] / ressize;
00398           accumulate[k] = 0.0;
00399         }
00400         residue = mol->atom(i)->uniq_resid;
00401         ressize = 0;
00402         rescount++;
00403       }
00404       for (k = 0; k < 3; k++) {
00405         accumulate[k] += avpos[3*j + k];
00406       }
00407       j++;
00408       ressize++;
00409     }
00410   }
00411   for (k = 0; k < 3; k++) {
00412     avpos[3*rescount + k] = accumulate[k] / ressize;
00413   }
00414   rescount++;
00415   return rescount;
00416 }
00417 
00418 // Calculate average position of selected atoms over selected frames
00419 extern int measure_avpos(const AtomSel *sel, MoleculeList *mlist, 
00420                          int start, int end, int step, float *avpos) {
00421   if (!sel)                     return MEASURE_ERR_NOSEL;
00422   if (sel->num_atoms == 0)      return MEASURE_ERR_NOATOMS;
00423 
00424   Molecule *mymol = mlist->mol_from_id(sel->molid());
00425   int maxframes = mymol->numframes();
00426   
00427   // accept value of -1 meaning "all" frames
00428   if (end == -1)
00429     end = maxframes-1;
00430 
00431   if (maxframes == 0 || start < 0 || start > end || 
00432       end >= maxframes || step <= 0)
00433     return MEASURE_ERR_BADFRAMERANGE;
00434 
00435   long i;
00436   for (i=0; i<(3L*sel->selected); i++)
00437     avpos[i] = 0.0f;
00438 
00439   long frame, avcount, j;
00440   for (avcount=0,frame=start; frame<=end; avcount++,frame+=step) {
00441     const float *framepos = (mymol->get_frame(frame))->pos;
00442     for (j=0,i=sel->firstsel; i<=sel->lastsel; i++) {
00443       if (sel->on[i]) {
00444         avpos[j*3L    ] += framepos[i*3L    ];
00445         avpos[j*3L + 1] += framepos[i*3L + 1];
00446         avpos[j*3L + 2] += framepos[i*3L + 2];
00447         j++;
00448       }
00449     } 
00450   }
00451 
00452   float avinv = 1.0f / (float) avcount;
00453   for (j=0; j<(3L*sel->selected); j++) {
00454     avpos[j] *= avinv;
00455   } 
00456 
00457   return MEASURE_NOERR;
00458 }
00459 
00460 
00461 // Calculate dipole moment for selected atoms
00462 extern int measure_dipole(const AtomSel *sel, MoleculeList *mlist, 
00463                           float *dipole, int unitsdebye, int usecenter) {
00464   if (!sel)                     return MEASURE_ERR_NOSEL;
00465   if (sel->num_atoms == 0)      return MEASURE_ERR_NOATOMS;
00466 
00467   Molecule *mymol = mlist->mol_from_id(sel->molid());
00468   double  rvec[3] = {0, 0, 0};
00469   double qrvec[3] = {0, 0, 0};
00470   double mrvec[3] = {0, 0, 0};
00471   double totalq = 0.0;
00472   double totalm = 0.0;
00473   int i;
00474 
00475   // get atom coordinates
00476   const float *framepos = sel->coordinates(mlist);
00477 
00478   // get atom charges
00479   const float *q = mymol->charge();
00480   const float *m = mymol->mass();
00481 
00482   for (i=sel->firstsel; i<=sel->lastsel; i++) {
00483     if (sel->on[i]) {
00484       int ind = i * 3;
00485       rvec[0] += framepos[ind    ];
00486       rvec[1] += framepos[ind + 1];
00487       rvec[2] += framepos[ind + 2];
00488 
00489       qrvec[0] += framepos[ind    ] * q[i];
00490       qrvec[1] += framepos[ind + 1] * q[i];
00491       qrvec[2] += framepos[ind + 2] * q[i];
00492 
00493       mrvec[0] += framepos[ind    ] * m[i];
00494       mrvec[1] += framepos[ind + 1] * m[i];
00495       mrvec[2] += framepos[ind + 2] * m[i];
00496 
00497       totalq += q[i];
00498       totalm += m[i];
00499     }
00500   }
00501 
00502   // fall back to geometrical center when bad or no masses
00503   if (totalm < 0.0001)
00504     usecenter=1;
00505   
00506   switch (usecenter) {
00507     case 1:
00508     {
00509         double rscale = totalq / sel->selected; 
00510         dipole[0] = (float) (qrvec[0] - (rvec[0] * rscale)); 
00511         dipole[1] = (float) (qrvec[1] - (rvec[1] * rscale)); 
00512         dipole[2] = (float) (qrvec[2] - (rvec[2] * rscale)); 
00513         break;
00514     }
00515 
00516     case -1: 
00517     {
00518         double mscale = totalq / totalm; 
00519         dipole[0] = (float) (qrvec[0] - (mrvec[0] * mscale)); 
00520         dipole[1] = (float) (qrvec[1] - (mrvec[1] * mscale)); 
00521         dipole[2] = (float) (qrvec[2] - (mrvec[2] * mscale)); 
00522         break;
00523     }
00524     
00525     case 0: // fallthrough
00526     default: 
00527     {
00528         dipole[0] = (float) qrvec[0];
00529         dipole[1] = (float) qrvec[1];
00530         dipole[2] = (float) qrvec[2];
00531         break;
00532     }
00533   }
00534 
00535   // According to the values in
00536   // http://www.physics.nist.gov/cuu/Constants/index.html
00537   // 1 e*A = 4.80320440079 D
00538   // 1 D = 1E-18 Fr*cm = 0.208194346224 e*A
00539   // 1 e*A = 1.60217653 E-29 C*m
00540   // 1 C*m = 6.24150947961 E+28 e*A
00541   // 1 e*A = 1.88972613458 e*a0
00542   // 1 e*a0 = 0.529177208115 e*A
00543 
00544   if (unitsdebye) {
00545     // 1 e*A = 4.80320440079 D 
00546     // latest CODATA (2006) gives:
00547     //         4.80320425132073913031
00548     dipole[0] *= 4.80320425132f;
00549     dipole[1] *= 4.80320425132f;
00550     dipole[2] *= 4.80320425132f;
00551   }
00552  
00553   return MEASURE_NOERR;
00554 }
00555 
00556 
00557 extern int measure_hbonds(Molecule *mol, AtomSel *sel1, AtomSel *sel2, double cutoff, double maxangle, int *donlist, int *hydlist, int *acclist, int maxsize) {
00558   int hbondcount = 0;
00559   const float *pos = sel1->coordinates(mol->app->moleculeList);
00560   const int *A = sel1->on;
00561   const int *B = sel2 ? sel2->on : sel1->on;
00562   GridSearchPair *pairlist = vmd_gridsearch2(pos, sel1->num_atoms, A, B, (float) cutoff, sel1->num_atoms * 27);
00563   GridSearchPair *p, *tmp;
00564   float donortoH[3], Htoacceptor[3];
00565   
00566   for (p=pairlist; p != NULL; p=tmp) {
00567     MolAtom *a1 = mol->atom(p->ind1); 
00568     MolAtom *a2 = mol->atom(p->ind2); 
00569     
00570     // neither the donor nor acceptor may be hydrogens
00571     if (mol->atom(p->ind1)->atomType == ATOMHYDROGEN ||
00572         mol->atom(p->ind2)->atomType == ATOMHYDROGEN) {
00573       tmp = p->next;
00574       free(p);
00575       continue;
00576     } 
00577     if (!a1->bonded(p->ind2)) {
00578       int b1 = a1->bonds;
00579       int b2 = a2->bonds;
00580       const float *coor1 = pos + 3*p->ind1;
00581       const float *coor2 = pos + 3*p->ind2;
00582       int k;
00583       // first treat sel1 as donor
00584       for (k=0; k<b1; k++) {
00585         const int hindex = a1->bondTo[k];
00586         if (mol->atom(hindex)->atomType == ATOMHYDROGEN) {         
00587           const float *hydrogen = pos + 3*hindex;
00588           vec_sub(donortoH,hydrogen,coor1);
00589           vec_sub(Htoacceptor,coor2,hydrogen);
00590           if (angle(donortoH, Htoacceptor)  < maxangle ) {
00591             if (hbondcount < maxsize) {
00592               donlist[hbondcount] = p->ind1;
00593               acclist[hbondcount] = p->ind2;
00594               hydlist[hbondcount] = hindex;
00595             }
00596             hbondcount++;
00597           }
00598         }
00599       }
00600       // if only one atom selection was given, treat sel2 as a donor as well
00601       if (!sel2) {
00602         for (k=0; k<b2; k++) {
00603           const int hindex = a2->bondTo[k];
00604           if (mol->atom(hindex)->atomType == ATOMHYDROGEN) {
00605             const float *hydrogen = pos + 3*hindex;
00606             vec_sub(donortoH,hydrogen,coor2);
00607             vec_sub(Htoacceptor,coor1,hydrogen);
00608             if (angle(donortoH, Htoacceptor)  < maxangle ) {
00609               if (hbondcount < maxsize) {
00610                 donlist[hbondcount] = p->ind2;
00611                 acclist[hbondcount] = p->ind1;
00612                 hydlist[hbondcount] = hindex;
00613               }
00614               hbondcount++;
00615             }
00616           }
00617         }
00618       } 
00619     }
00620     tmp = p->next;
00621     free(p);
00622   }
00623   return hbondcount;
00624 }
00625 
00626 int measure_rmsf_perresidue(const AtomSel *sel, MoleculeList *mlist, 
00627                         int start, int end, int step, float *rmsf) {
00628   if (!sel)                     return MEASURE_ERR_NOSEL;
00629   if (sel->num_atoms == 0)      return MEASURE_ERR_NOATOMS;
00630 
00631   Molecule *mymol = mlist->mol_from_id(sel->molid());
00632   int maxframes = mymol->numframes();
00633 
00634   // accept value of -1 meaning "all" frames
00635   if (end == -1)
00636     end = maxframes-1;
00637 
00638   if (maxframes == 0 || start < 0 || start > end ||
00639       end >= maxframes || step <= 0)
00640     return MEASURE_ERR_BADFRAMERANGE;
00641 
00642   int i;
00643   for (i=0; i<sel->selected; i++)
00644     rmsf[i] = 0.0f;
00645 
00646   int rc; 
00647   float *avpos = new float[3*sel->selected];
00648   //rc will be the number of residues.
00649   rc = measure_avpos_perresidue(sel, mlist, start, end, step, avpos);
00650   if (rc <= 0) {
00651     delete [] avpos;
00652     return rc;
00653   }
00654   
00655   int frame, avcount;
00656   float *center = new float[3*rc];
00657   float *weights = new float[sel->selected];
00658   for (i = 0; i < sel->selected; i++) {
00659     weights[i] = 1.0;
00660   }
00661   // calculate per-residue variance here. Its a simple calculation, since we have measure_center_perresidue.
00662   for (avcount=0,frame=start; frame<=end; avcount++,frame+=step) {
00663     const float *framepos = (mymol->get_frame(frame))->pos;
00664     measure_center_perresidue(mlist, sel, framepos, weights, center);
00665     for (i = 0; i < rc; i++) {
00666       rmsf[i] += distance2(&avpos[3*i], &center[3*i]);
00667     }
00668   }
00669   delete [] center;
00670   delete [] weights;
00671 
00672   float avinv = 1.0f / (float) avcount;
00673   for (i=0; i<rc; i++) {
00674     rmsf[i] = sqrtf(rmsf[i] * avinv);
00675   }
00676 
00677   delete [] avpos;
00678   return rc;
00679 }
00680 
00681 // Calculate RMS fluctuation of selected atoms over selected frames
00682 extern int measure_rmsf(const AtomSel *sel, MoleculeList *mlist, 
00683                         int start, int end, int step, float *rmsf) {
00684   if (!sel)                     return MEASURE_ERR_NOSEL;
00685   if (sel->num_atoms == 0)      return MEASURE_ERR_NOATOMS;
00686 
00687   Molecule *mymol = mlist->mol_from_id(sel->molid());
00688   int maxframes = mymol->numframes();
00689 
00690   // accept value of -1 meaning "all" frames
00691   if (end == -1)
00692     end = maxframes-1;
00693 
00694   if (maxframes == 0 || start < 0 || start > end ||
00695       end >= maxframes || step <= 0)
00696     return MEASURE_ERR_BADFRAMERANGE;
00697 
00698   int i;
00699   for (i=0; i<sel->selected; i++)
00700     rmsf[i] = 0.0f;
00701 
00702   int rc; 
00703   float *avpos = new float[3L*sel->selected];
00704   rc = measure_avpos(sel, mlist, start, end, step, avpos);
00705 
00706   if (rc != MEASURE_NOERR) {
00707     delete [] avpos;
00708     return rc;
00709   }
00710 
00711   // calculate per-atom variance here
00712   int frame, avcount, j;
00713   for (avcount=0,frame=start; frame<=end; avcount++,frame+=step) {
00714     const float *framepos = (mymol->get_frame(frame))->pos;
00715     for (j=0,i=sel->firstsel; i<=sel->lastsel; i++) {
00716       if (sel->on[i]) {
00717         rmsf[j] += distance2(&avpos[3L*j], &framepos[3L*i]);
00718         j++;
00719       }
00720     }
00721   }
00722 
00723   float avinv = 1.0f / (float) avcount;
00724   for (j=0; j<sel->selected; j++) {
00725     rmsf[j] = sqrtf(rmsf[j] * avinv);
00726   }
00727 
00728   delete [] avpos;
00729   return MEASURE_NOERR;
00730 }
00731 
00732 
00734 //  rgyr := sqrt(sum (mass(n) ( r(n) - r(com) )^2)/sum(mass(n)))
00735 //  The return value, a float, is put in 'float *rgyr'
00736 //  The function return value is 0 if ok, <0 if not
00737 int measure_rgyr(const AtomSel *sel, MoleculeList *mlist, const float *weight, 
00738                  float *rgyr) {
00739   if (!rgyr) return MEASURE_ERR_NORGYR;
00740   if (!sel) return MEASURE_ERR_NOSEL;
00741   if (!mlist) return MEASURE_ERR_GENERAL;
00742 
00743   const float *framepos = sel->coordinates(mlist);
00744 
00745   // compute the center of mass with the current weights
00746   float com[3];
00747   int ret_val = measure_center(sel, framepos, weight, com);
00748   if (ret_val < 0) 
00749     return ret_val;
00750   
00751   // measure center of gyration
00752   int i, j;
00753   float total_w, sum;
00754 
00755   total_w=sum=0;
00756   for (j=0,i=sel->firstsel; i<=sel->lastsel; i++) {
00757     if (sel->on[i]) {
00758       float w = weight[j];
00759       total_w += w;
00760       sum += w * distance2(framepos + 3L*i, com);
00761       j++;
00762     } 
00763   }
00764 
00765   if (total_w == 0) {
00766     return MEASURE_ERR_BADWEIGHTSUM;
00767   }
00768 
00769   // and finalize the computation
00770   *rgyr = sqrtf(sum/total_w);
00771 
00772   return MEASURE_NOERR;
00773 }
00774 
00775 /*I'm going to assume that *rmsd points to an array the length of sel1. The return
00776 value will indicate the ACTUAL number of residues in the selection (specifically sel1),
00777 which tells the caller how many values should be in the returned list, or a number
00778 less than zero on error.
00779 */
00780 int measure_rmsd_perresidue(const AtomSel *sel1, const AtomSel *sel2, MoleculeList *mlist,
00781                  int num,
00782                  float *weight, float *rmsd) {
00783   if (!sel1 || !sel2)   return MEASURE_ERR_NOSEL;
00784   if (sel1->selected < 1 || sel2->selected < 1) return MEASURE_ERR_NOSEL;
00785   if (!weight || !rmsd) return MEASURE_ERR_NOWEIGHT;
00786 
00787   // the number of selected atoms must be the same
00788   if (sel1->selected != sel2->selected) return MEASURE_ERR_MISMATCHEDCNT;
00789 
00790   // need to know how to traverse the list of weights
00791   // there could be 1 weight per atom (sel_flg == 1) or 
00792   // 1 weight per selected atom (sel_flg == 0)
00793   int sel_flg;
00794   if (num == sel1->num_atoms) {
00795     sel_flg = 1; // using all elements
00796   } else {
00797     sel_flg = 0; // using elements from selection
00798   }
00799 
00800   // temporary variables
00801   float tmp_w;
00802   int w_index = 0;              // the term in the weight field to use
00803   int sel1ind = sel1->firstsel; // start from the first selected atom
00804   int sel2ind = sel2->firstsel; // start from the first selected atom
00805   float wsum = 0;
00806   float twsum = 0;
00807   float rmsdsum = 0;
00808   const float *framepos1 = sel1->coordinates(mlist);
00809   const float *framepos2 = sel2->coordinates(mlist);
00810   Molecule *mol = mlist->mol_from_id(sel1->molid());
00811 
00812 
00813   int count = sel1->selected;
00814   int residue = mol->atom(sel1ind)->uniq_resid;
00815   int rescount = 0;
00816   while (count-- > 0) {
00817     // find next 'on' atom in sel1 and sel2
00818     // loop is safe since we already stop the on count > 0 above
00819     while (!sel1->on[sel1ind])
00820       sel1ind++;
00821     while (!sel2->on[sel2ind])
00822       sel2ind++;
00823     if (residue != mol->atom(sel1ind)->uniq_resid) {
00824       rmsd[rescount] = sqrtf(rmsdsum / wsum);
00825       rmsdsum = 0;
00826       twsum += wsum;
00827       wsum = 0;
00828       residue = mol->atom(sel1ind)->uniq_resid;
00829       rescount++;
00830     }
00831     // the weight offset to use depends on how many terms there are
00832     if (sel_flg == 0) {
00833       tmp_w = weight[w_index++];
00834     } else {
00835       tmp_w = weight[sel1ind]; // use the first selection for the weights
00836     }
00837 
00838     // sum the calculated rmsd and weight values
00839     rmsdsum += tmp_w * distance2(framepos1 + 3*sel1ind, framepos2 + 3*sel2ind);
00840     wsum += tmp_w;
00841 
00842     // and advance to the next atom pair
00843     sel1ind++;
00844     sel2ind++;
00845   }
00846   twsum += wsum;
00847   // check weight sum
00848   if (twsum == 0) {
00849     return MEASURE_ERR_BADWEIGHTSUM;
00850   }
00851 
00852   // finish the rmsd calcs
00853   rmsd[rescount++] = sqrtf(rmsdsum / wsum);
00854 
00855   return rescount; // and say rmsd is OK
00856 }
00857 
00859 //  1) if num == sel.selected ; assumes there is one weight per 
00860 //           selected atom
00861 //  2) if num == sel.num_atoms; assumes weight[i] is for atom[i]
00862 //  returns 0 and value in rmsd if good
00863 //   return < 0 if invalid
00864 //  Function is::=  rmsd = 
00865 //    sqrt(sum(weight(n) * sqr(r1(i(n))-r2(i(n))))/sum(weight(n)) / N
00866 int measure_rmsd(const AtomSel *sel1, const AtomSel *sel2,
00867                  int num, const float *framepos1, const float *framepos2,
00868                  float *weight, float *rmsd) {
00869   if (!sel1 || !sel2)   return MEASURE_ERR_NOSEL;
00870   if (sel1->selected < 1 || sel2->selected < 1) return MEASURE_ERR_NOSEL;
00871   if (!weight || !rmsd) return MEASURE_ERR_NOWEIGHT;
00872 
00873   // the number of selected atoms must be the same
00874   if (sel1->selected != sel2->selected) return MEASURE_ERR_MISMATCHEDCNT;
00875 
00876   // need to know how to traverse the list of weights
00877   // there could be 1 weight per atom (sel_flg == 1) or 
00878   // 1 weight per selected atom (sel_flg == 0)
00879   int sel_flg;
00880   if (num == sel1->num_atoms) {
00881     sel_flg = 1; // using all elements
00882   } else {
00883     sel_flg = 0; // using elements from selection
00884   }
00885 
00886   // temporary variables
00887   float tmp_w;
00888   int w_index = 0;              // the term in the weight field to use
00889   int sel1ind = sel1->firstsel; // start from the first selected atom
00890   int sel2ind = sel2->firstsel; // start from the first selected atom
00891   float wsum = 0;
00892   float rmsdsum = 0;
00893 
00894   *rmsd = 10000000; // if we bail out, return a huge value
00895 
00896   // compute the rmsd
00897   int count = sel1->selected;
00898   while (count-- > 0) {
00899     // find next 'on' atom in sel1 and sel2
00900     // loop is safe since we already stop the on count > 0 above
00901     while (!sel1->on[sel1ind])
00902       sel1ind++;
00903     while (!sel2->on[sel2ind])
00904       sel2ind++;
00905 
00906     // the weight offset to use depends on how many terms there are
00907     if (sel_flg == 0) {
00908       tmp_w = weight[w_index++];
00909     } else {
00910       tmp_w = weight[sel1ind]; // use the first selection for the weights
00911     }
00912 
00913     // sum the calculated rmsd and weight values
00914     rmsdsum += tmp_w * distance2(framepos1 + 3L*sel1ind, framepos2 + 3L*sel2ind);
00915     wsum += tmp_w;
00916 
00917     // and advance to the next atom pair
00918     sel1ind++;
00919     sel2ind++;
00920   }
00921 
00922   // check weight sum
00923   if (wsum == 0) {
00924     return MEASURE_ERR_BADWEIGHTSUM;
00925   }
00926 
00927   // finish the rmsd calcs
00928   *rmsd = sqrtf(rmsdsum / wsum);
00929 
00930   return MEASURE_NOERR; // and say rmsd is OK
00931 }
00932 
00933 /* jacobi.C, taken from Numerical Recipes and specialized to 3x3 case */
00934 
00935 #define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\
00936         a[k][l]=h+s*(g-h*tau);
00937 
00938 static int jacobi(float a[4][4], float d[3], float v[3][3])
00939 {
00940   int n=3;
00941   int j,iq,ip,i;
00942   float tresh,theta,tau,t,sm,s,h,g,c,*b,*z;
00943   
00944   b=new float[n];
00945   z=new float[n];
00946   for (ip=0;ip<n;ip++) {
00947     for (iq=0;iq<n;iq++) v[ip][iq]=0.0;
00948     v[ip][ip]=1.0;
00949   }
00950   for (ip=0;ip<n;ip++) {
00951     b[ip]=d[ip]=a[ip][ip];
00952     z[ip]=0.0;
00953   }
00954   for (i=1;i<=50;i++) {
00955     sm=0.0;
00956     for (ip=0;ip<n-1;ip++) {
00957       for (iq=ip+1;iq<n;iq++)
00958         sm += (float) fabs(a[ip][iq]);
00959     }
00960     if (sm == 0.0) {
00961       delete [] z;
00962       delete [] b;
00963       return 0; // Exit normally
00964     }
00965     if (i < 4)
00966       tresh=0.2f*sm/(n*n);
00967     else
00968       tresh=0.0f;
00969     for (ip=0;ip<n-1;ip++) {
00970       for (iq=ip+1;iq<n;iq++) {
00971         g=100.0f * fabsf(a[ip][iq]);
00972         if (i > 4 && (float)(fabs(d[ip])+g) == (float)fabs(d[ip])
00973             && (float)(fabs(d[iq])+g) == (float)fabs(d[iq]))
00974           a[ip][iq]=0.0;
00975         else if (fabs(a[ip][iq]) > tresh) {
00976           h=d[iq]-d[ip];
00977           if ((float)(fabs(h)+g) == (float)fabs(h))
00978             t=(a[ip][iq])/h;
00979           else {
00980             theta=0.5f*h/(a[ip][iq]);
00981             t=1.0f/(fabsf(theta)+sqrtf(1.0f+theta*theta));
00982             if (theta < 0.0f) t = -t;
00983           }
00984           c=1.0f/sqrtf(1+t*t);
00985           s=t*c;
00986           tau=s/(1.0f+c);
00987           h=t*a[ip][iq];
00988           z[ip] -= h;
00989           z[iq] += h;
00990           d[ip] -= h;
00991           d[iq] += h;
00992           a[ip][iq]=0.0;
00993           for (j=0;j<=ip-1;j++) {
00994             ROTATE(a,j,ip,j,iq)
00995               }
00996           for (j=ip+1;j<=iq-1;j++) {
00997             ROTATE(a,ip,j,j,iq)
00998               }
00999           for (j=iq+1;j<n;j++) {
01000             ROTATE(a,ip,j,iq,j)
01001               }
01002           for (j=0;j<n;j++) {
01003             ROTATE(v,j,ip,j,iq)
01004               }
01005         }
01006       }
01007     }
01008     for (ip=0;ip<n;ip++) {
01009       b[ip] += z[ip];
01010       d[ip]=b[ip];
01011       z[ip]=0.0;
01012     }
01013   }
01014   delete [] b;
01015   delete [] z;
01016   return 1; // Failed to converge
01017 }
01018 #undef ROTATE
01019 
01020 static int transvecinv(const double *v, Matrix4 &res) {
01021   double x, y, z;
01022   x=v[0];
01023   y=v[1];
01024   z=v[2];
01025   if (x == 0.0 && y == 0.0) {
01026     if (z == 0.0) {
01027       return -1;
01028     }
01029     if (z > 0)
01030       res.rot(90, 'y');
01031     else
01032       res.rot(-90, 'y');
01033     return 0;
01034   }
01035   double theta = atan2(y,x);
01036   double length = sqrt(x*x + y*y);
01037   double phi = atan2(z,length);
01038   res.rot((float) RADTODEG(phi),  'y');
01039   res.rot((float) RADTODEG(-theta), 'z');
01040   return 0;
01041 }
01042  
01043 static int transvec(const double *v, Matrix4 &res) {
01044   double x, y, z;
01045   x=v[0];
01046   y=v[1];
01047   z=v[2];
01048   if (x == 0.0 && y == 0.0) {
01049     if (z == 0.0) {
01050       return -1;
01051     }
01052     if (z > 0)
01053       res.rot(-90, 'y');
01054     else
01055       res.rot(90, 'y');
01056     return 0;
01057   }
01058   double theta = atan2(y,x);
01059   double length = sqrt(x*x + y*y);
01060   double phi = atan2(z,length);
01061   res.rot((float) RADTODEG(theta), 'z');
01062   res.rot((float) RADTODEG(-phi),  'y');
01063   return 0;
01064 }
01065 
01066 static Matrix4 myfit2(const float *x, const float *y, 
01067                const float *comx, const float *comy) {
01068 
01069   Matrix4 res;
01070   double dx[3], dy[3];
01071   dx[0] = x[0] - comx[0];
01072   dx[1] = x[1] - comx[1];
01073   dx[2] = x[2] - comx[2];
01074   dy[0] = y[0] - comy[0];
01075   dy[1] = y[1] - comy[1];
01076   dy[2] = y[2] - comy[2];
01077   
01078   res.translate(comy[0], comy[1], comy[2]);
01079   transvec(dy, res);
01080   transvecinv(dx, res);
01081   res.translate(-comx[0], -comx[1], -comx[2]);
01082   return res;
01083 }
01084 
01085 static Matrix4 myfit3(const float *x1, const float *x2, 
01086                       const float *y1, const float *y2,
01087                       const float *comx, const float *comy) {
01088    
01089   Matrix4 mx, my, rx1, ry1;
01090   double dx1[3], dy1[3], angle;
01091   float dx2[3], dy2[3], x2t[3], y2t[3];
01092 
01093   for (int i=0; i<3; i++) {
01094     dx1[i] = x1[i] - comx[i];
01095     dx2[i] = x2[i] - comx[i];
01096     dy1[i] = y1[i] - comy[i];
01097     dy2[i] = y2[i] - comy[i];
01098   }
01099 
01100   // At some point, multmatrix went from pre-multiplying, as the code of 
01101   // Matrix.C itself suggests, to post multiplying, which is what the 
01102   // users must have expected. Thus my.multmatrix(mx) is the same as 
01103   // my = my * mx, not mx * my.  This means that you use the matrix 
01104   // conventions of openGL (first matrix in the code is the last 
01105   // matrix applied)
01106   transvecinv(dx1, rx1);
01107   rx1.multpoint3d(dx2, x2t);
01108   angle = atan2(x2t[2], x2t[1]);
01109   mx.rot((float) RADTODEG(angle), 'x'); 
01110   mx.multmatrix(rx1);
01111   mx.translate(-comx[0], -comx[1], -comx[2]);
01112   
01113   transvecinv(dy1, ry1);
01114   ry1.multpoint3d(dy2, y2t);
01115   angle = atan2(y2t[2], y2t[1]);
01116   my.rot((float) RADTODEG(angle), 'x');
01117   my.multmatrix(ry1);
01118   my.translate(-comy[0], -comy[1], -comy[2]);
01119   my.inverse();
01120   my.multmatrix(mx);
01121   return my;
01122 }
01123 
01124 // find the best fit alignment to take the first structure into the second
01125 // Put the result in the matrix 'mat'
01126 //  This algorithm comes from Kabsch, Acta Cryst. (1978) A34, 827-828.
01127 // Need the 2nd weight for the com calculation
01128 int measure_fit(const AtomSel *sel1, const AtomSel *sel2, const float *x, 
01129                 const float *y, const float *weight, 
01130                 const int *order, Matrix4 *mat) {
01131   float comx[3];
01132   float comy[3];
01133   int num = sel1->selected;
01134   int ret_val;
01135   ret_val = measure_center(sel1, x, weight, comx);
01136   if (ret_val < 0) {
01137     return ret_val;
01138   }
01139   ret_val = measure_center(sel2, y, weight, comy);
01140   if (ret_val < 0) {
01141     return ret_val;
01142   }
01143 
01144   // the Kabsch method won't work of the number of atoms is less than 4
01145   // (and won't work in some cases of n > 4; I think it works so long as
01146   // three or more planes are needed to intersect all the data points
01147   switch (sel1->selected) {
01148   case 1: { // simple center of mass alignment
01149     Matrix4 tmp;
01150     tmp.translate(-comx[0], -comx[1], -comx[2]);
01151     tmp.translate(comy[0], comy[1], comy[2]);
01152     memcpy(mat->mat, tmp.mat, 16L*sizeof(float));
01153     return MEASURE_NOERR;
01154   }
01155   case 3:
01156   case 2: { 
01157     // find the first (n-1) points (from each molecule)
01158     int pts[6], count = 0;
01159     int n;
01160     for (n=sel1->firstsel; n<=sel1->lastsel; n++) {
01161       if (sel1->on[n]) {
01162         pts[count++] = n;
01163         if (sel1->selected == 2) {
01164           count++;                   // will put y data in pts[3]
01165           break;
01166         }
01167         if (count == 2) break;
01168       }
01169     }
01170     for (n=sel2->firstsel; n<=sel2->lastsel; n++) {
01171       if (sel2->on[n]) {
01172         pts[count++] = n;
01173         if (sel1->selected == 2) {
01174           count++;
01175           break;
01176         }
01177         if (count == 4) break;
01178       }
01179     }
01180     if (count != 4) {
01181       return MEASURE_ERR_MISMATCHEDCNT;
01182     }
01183     
01184     // reorder the sel2 atoms according to the order parameter
01185     if (order != NULL) {
01186       int i; 
01187       int tmp[6];
01188       memcpy(tmp, pts, sizeof(pts));
01189       for (i=0; i<num; i++) {
01190         pts[i + num] = tmp[num + order[i]]; // order indices are 0-based
01191       } 
01192     }
01193 
01194     if (sel1->selected == 2) {
01195       *mat = myfit2(x+3L*pts[0], y+3L*pts[2], comx, comy);
01196       ret_val = 0;
01197     } else {
01198       *mat = myfit3(x+3L*pts[0], x+3L*pts[1], y+3L*pts[2], y+3L*pts[3], comx, comy);
01199       ret_val = 0;
01200     }  
01201     if (ret_val != 0) {
01202       return MEASURE_ERR_GENERAL;
01203     }
01204 
01205     return 0;
01206   }
01207   default:
01208     break;
01209   }
01210   // at this point I know all the data values are good
01211  
01212 
01213   // use the new RMS fit implementation by default unless told otherwise 
01214   char *opt = getenv("VMDRMSFITMETHOD");
01215   if (!opt || strcmp(opt, "oldvmd")) {
01216     int i, k;
01217     float *v1, *v2;
01218     v1 = new float[3L*num];
01219     v2 = new float[3L*num];
01220     for (k=0, i=sel1->firstsel; i<=sel1->lastsel; i++) {
01221       if (sel1->on[i]) {
01222         long ind = 3L * i;
01223         v1[k++] = x[ind    ];
01224         v1[k++] = x[ind + 1];
01225         v1[k++] = x[ind + 2];
01226       }
01227     }
01228     for (k=0, i=sel2->firstsel; i<=sel2->lastsel; i++) {
01229       if (sel2->on[i]) {
01230         long ind = 3L * i;
01231         v2[k++] = y[ind    ];
01232         v2[k++] = y[ind + 1];
01233         v2[k++] = y[ind + 2];
01234       }
01235     }
01236 
01237     // reorder the sel2 atoms according to the order parameter
01238     if (order != NULL) {
01239       int i; 
01240       float *tmp = new float[3L*num];
01241       memcpy(tmp, v2, 3L*num*sizeof(float));
01242       for (i=0; i<num; i++) {
01243         long ind = 3L * i;
01244         long idx = 3L * order[i]; // order indices are 0-based
01245         v2[ind    ] = tmp[idx    ];
01246         v2[ind + 1] = tmp[idx + 1];
01247         v2[ind + 2] = tmp[idx + 2];
01248       } 
01249       delete[] tmp;
01250     }
01251 
01252     float tmp[16];
01253     // MatrixFitRMS returns RMS distance of fitted molecule.  Would be nice
01254     // to return this information to the user since it would make computing
01255     // the fitted RMSD much faster (no need to get the matrix, apply the
01256     // transformation, recompute RMS).
01257     MatrixFitRMS(num, v1, v2, weight, tmp);
01258 
01259     delete [] v1;
01260     delete [] v2;
01261     //fprintf(stderr, "got err %f\n", err);
01262     // output from FitRMS is a 3x3 matrix, plus a pre-translation stored in
01263     // row 3, and a post-translation stored in column 3.
01264     float pre[3], post[3];
01265     for (i=0; i<3; i++) {
01266       post[i] = tmp[4L*i+3];
01267       tmp[4L*i+3] = 0;
01268     }
01269     for (i=0; i<3; i++) {
01270       pre[i] = tmp[12+i];
01271       tmp[12+i] = 0;
01272     }
01273     Matrix4 result;
01274     result.translate(pre);
01275     result.multmatrix(Matrix4(tmp));
01276     result.translate(post);
01277     memcpy(mat->mat, result.mat, 16L*sizeof(float));
01278     return 0;
01279   }
01280 
01281   // the old RMS fit code doesn't support reordering of sel2 currently
01282   if (order != NULL) {
01283     return MEASURE_ERR_NOTSUP;
01284   }
01285 
01286   // a) compute R = r(i,j) = sum( w(n) * (y(n,i)-comy(i)) * (x(n,j)-comx(j)))
01287   Matrix4 R;
01288   int i,j;
01289   float scale = (float) num * num;
01290   for (i=0; i<3; i++) {
01291     for (j=0; j<3; j++) {
01292       float tmp = 0;
01293       int nx = 0, ny = 0, k = 0; 
01294       while (nx < sel1->num_atoms && ny < sel2->num_atoms) { 
01295         if (!sel1->on[nx]) {
01296           nx++;
01297           continue;
01298         }
01299         if (!sel2->on[ny]) {
01300           ny++;
01301           continue;
01302         }
01303 
01304         // found both, so get data
01305         
01306         tmp += weight[k] * (y[3L*ny+i] - comy[i]) * (x[3L*nx+j] - comx[j]) /
01307           scale;
01308         nx++;
01309         ny++;
01310         k++;
01311       }
01312       R.mat[4L*i+j] = tmp;
01313     }
01314   }
01315 
01316   // b) 1) form R~R
01317   Matrix4 Rt;
01318   for (i=0; i<3; i++) {
01319     for (j=0; j<3; j++) {
01320       Rt.mat[4L*i+j] = R.mat[4L*j+i];
01321     }
01322   }
01323   Matrix4 RtR(R);
01324   RtR.multmatrix(Rt);
01325 
01326   // b) 2) find the eigenvalues and eigenvectors
01327   float evalue[3];
01328   float evector[3][3];
01329   float tmpmat[4][4];
01330   for (i=0; i<4; i++)
01331     for (j=0; j<4; j++)
01332       tmpmat[i][j]=RtR.mat[4L*i+j];
01333 
01334   if(jacobi(tmpmat,evalue,evector) != 0) return MEASURE_ERR_NONZEROJACOBI;
01335   // transposition the evector matrix to put the vectors in rows
01336   float vectmp;
01337   vectmp=evector[0][1]; evector[0][1]=evector[1][0]; evector[1][0]=vectmp;
01338   vectmp=evector[0][2]; evector[0][2]=evector[2][0]; evector[2][0]=vectmp;
01339   vectmp=evector[2][1]; evector[2][1]=evector[1][2]; evector[1][2]=vectmp;
01340 
01341 
01342   // b) 4) sort so that the eigenvalues are from largest to smallest
01343   //      (or rather so a[0] is eigenvector with largest eigenvalue, ...)
01344   float *a[3];
01345   a[0] = evector[0];
01346   a[1] = evector[1];
01347   a[2] = evector[2];
01348 #define SWAP(qq,ww) {                                           \
01349     float v; float *v1;                                         \
01350     v = evalue[qq]; evalue[qq] = evalue[ww]; evalue[ww] = v;    \
01351     v1 = a[qq]; a[qq] = a[ww]; a[ww] = v1;                      \
01352 }
01353   if (evalue[0] < evalue[1]) {
01354     SWAP(0, 1);
01355   }
01356   if (evalue[0] < evalue[2]) {
01357     SWAP(0, 2);
01358   }
01359   if (evalue[1] < evalue[2]) {
01360     SWAP(1, 2);
01361   }
01362 
01363   
01364   // c) determine b(i) = R*a(i)
01365   float b[3][3];
01366 
01367   Rt.multpoint3d(a[0], b[0]);
01368   vec_normalize(b[0]);
01369 
01370   Rt.multpoint3d(a[1], b[1]);
01371   vec_normalize(b[1]);
01372 
01373   Rt.multpoint3d(a[2], b[2]);
01374   vec_normalize(b[2]);
01375 
01376   // d) compute U = u(i,j) = sum(b(k,i) * a(k,j))
01377   Matrix4 U;
01378   for (i=0; i<3; i++) {
01379     for (j=0; j<3; j++) {
01380       float *tmp = &(U.mat[4L*j+i]);
01381       *tmp = 0;
01382       for (int k=0; k<3; k++) {
01383         *tmp += b[k][i] * a[k][j];
01384       }
01385     }
01386   }
01387 
01388   // Check the determinant of U.  If it's negative, we need to
01389   // flip the sign of the last row.
01390   float *um = U.mat;
01391   float det = 
01392     um[0]*(um[4+1]*um[8+2] - um[4+2]*um[8+1]) -
01393     um[1]*(um[4+0]*um[8+2] - um[4+2]*um[8+0]) +
01394     um[2]*(um[4+0]*um[8+1] - um[4+1]*um[8+0]);
01395   if (det < 0) {
01396     for (int q=0; q<3; q++) um[8+q] = -um[8+q];
01397   }
01398 
01399   // e) apply the offset for com
01400   Matrix4 tx;
01401   tx.translate(-comx[0], -comx[1], -comx[2]);
01402   Matrix4 ty;
01403   ty.translate(comy[0], comy[1], comy[2]);
01404   //  U.multmatrix(com);
01405   ty.multmatrix(U);
01406   ty.multmatrix(tx);
01407   memcpy(mat->mat, ty.mat, 16L*sizeof(float));
01408   return MEASURE_NOERR;
01409 }
01410 
01411 // For different values of the random seed, the computed SASA's of brH.pdb 
01412 // converge to within 1% of each other when the number of points is about
01413 // 500.  We therefore use 500 as the default number.
01414 #define NPTS 500 
01415 
01416 extern int measure_sasa_perresidue(const AtomSel *sel, const float *framepos,
01417     const float *radius, float srad, float *sasa, 
01418     ResizeArray<float> *sasapts, const AtomSel *restrictsel,
01419     const int *nsamples, int *rescount,MoleculeList *mlist) {
01420 
01421   // check arguments
01422   if (!sel) return MEASURE_ERR_NOSEL;
01423   if (!sel->selected) {
01424     *sasa = 0;
01425     return MEASURE_NOERR;
01426   }
01427   if (!framepos) return MEASURE_ERR_NOFRAMEPOS;
01428   if (!radius)   return MEASURE_ERR_NORADII;
01429   if (restrictsel && restrictsel->num_atoms != sel->num_atoms)
01430     return MEASURE_ERR_MISMATCHEDCNT;
01431 
01432   Molecule *mol = mlist->mol_from_id(sel->molid());
01433   int residue = mol->atom(sel->firstsel)->uniq_resid;
01434   if (restrictsel)
01435     residue = mol->atom(restrictsel->firstsel)->uniq_resid;
01436   *rescount = 0;
01437 
01438   // If there is no restrictsel, set it to be the input selection. Otherwise, we would compute 
01439   //the SASA for each residue in isolation, which is almost certainly not what the user intended.
01440   if (!restrictsel)
01441     restrictsel = sel;
01442 
01443   int i;
01444   int npts = nsamples ? *nsamples : NPTS;
01445   float maxrad = -1;
01446 
01447 #if 0
01448   // Query min/max atom radii for the entire molecule
01449   mol->get_radii_minmax(minrad, maxrad);
01450 #endif
01451 
01452   // find biggest atom radius 
01453   for (i=0; i<sel->num_atoms; i++) {
01454     float rad = radius[i];
01455     if (maxrad < rad) maxrad = rad;
01456   }
01457 
01458   // Find atoms within maxrad of each other.  
01459   // build a list of pairs for each atom
01460   ResizeArray<int> *pairlist = new ResizeArray<int>[sel->num_atoms];
01461   {
01462     GridSearchPair *pairs;
01463     pairs = vmd_gridsearch1(framepos, sel->num_atoms, sel->on, 
01464                             2.0f * (maxrad + srad), 0, sel->num_atoms * 1000);
01465 
01466     GridSearchPair *p, *tmp; 
01467     for (p = pairs; p != NULL; p = tmp) {
01468       int ind1=p->ind1;
01469       int ind2=p->ind2;
01470       pairlist[ind1].append(ind2);
01471       pairlist[ind2].append(ind1);
01472       tmp = p->next;
01473       free(p);
01474     }
01475   }
01476 
01477   static const float RAND_MAX_INV = 1.0f/float(VMD_RAND_MAX);
01478   // Seed the random number generator before each calculation.  This gives
01479   // reproducible results and still allows a more accurate answer to be
01480   // obtained by increasing the samples size.  I don't know if this is a
01481   // "good" seed value or not, I just picked something random-looking.
01482   vmd_srandom(38572111);
01483 
01484   // All the spheres use the same random points.  
01485   float *spherepts = new float[3L*npts];
01486   for (i=0; i<npts; i++) {
01487     float u1 = (float) vmd_random();
01488     float u2 = (float) vmd_random();
01489     float z = 2.0f*u1*RAND_MAX_INV -1.0f;
01490     float phi = (float) (2.0f*VMD_PI*u2*RAND_MAX_INV);
01491     float R = sqrtf(1.0f-z*z);
01492     float sphi, cphi;
01493     sincosf(phi, &sphi, &cphi);
01494     spherepts[3L*i  ] = R*cphi;
01495     spherepts[3L*i+1] = R*sphi;
01496     spherepts[3L*i+2] = z;
01497   }
01498 
01499   const float prefac = (float) (4 * VMD_PI / npts);
01500   float totarea = 0.0f;
01501   // compute area for each atom based on its pairlist
01502   for (i=sel->firstsel; i<=sel->lastsel; i++) {
01503     if (sel->on[i]) {
01504       // only atoms in restrictsel contribute
01505       if (restrictsel && !restrictsel->on[i]) continue;
01506       if (residue != mol->atom(i)->uniq_resid) {
01507         //We are in a new residue, so we need to start modifying the value again.
01508         sasa[*rescount] = totarea;
01509         totarea = 0.0f;
01510         (*rescount)++;
01511         residue = mol->atom(i)->uniq_resid;
01512       }
01513       
01514       const float *loc = framepos+3L*i;
01515       float rad = radius[i]+srad;
01516       float surfpos[3];
01517       int surfpts = npts;
01518       const ResizeArray<int> &nbrs = pairlist[i];
01519       for (int j=0; j<npts; j++) {
01520         surfpos[0] = loc[0] + rad*spherepts[3L*j  ];
01521         surfpos[1] = loc[1] + rad*spherepts[3L*j+1];
01522         surfpos[2] = loc[2] + rad*spherepts[3L*j+2];
01523         int on = 1;
01524         for (int k=0; k<nbrs.num(); k++) {
01525           int ind = nbrs[k];
01526           const float *nbrloc = framepos+3L*ind;
01527           float radsq = radius[ind]+srad; radsq *= radsq;
01528           float dx = surfpos[0]-nbrloc[0];
01529           float dy = surfpos[1]-nbrloc[1];
01530           float dz = surfpos[2]-nbrloc[2];
01531           if (dx*dx + dy*dy + dz*dz <= radsq) {
01532             on = 0;
01533             break;
01534           }
01535         }
01536         if (on) {
01537           if (sasapts) {
01538             sasapts->append3(&surfpos[0]);
01539           }
01540         } else {
01541           surfpts--;
01542         }
01543       }
01544       float atomarea = prefac * rad * rad * surfpts;
01545       totarea += atomarea;
01546     }
01547   }
01548 
01549   delete [] pairlist;
01550   delete [] spherepts;
01551   sasa[*rescount] = totarea;
01552   (*rescount)++;
01553   return 0;
01554 }
01555 
01556 extern int measure_sasa(const AtomSel *sel, const float *framepos,
01557     const float *radius, float srad, float *sasa, 
01558     ResizeArray<float> *sasapts, const AtomSel *restrictsel,
01559     const int *nsamples) {
01560 
01561   // check arguments
01562   if (!sel) return MEASURE_ERR_NOSEL;
01563   if (!sel->selected) {
01564     *sasa = 0;
01565     return MEASURE_NOERR;
01566   }
01567   if (!framepos) return MEASURE_ERR_NOFRAMEPOS;
01568   if (!radius)   return MEASURE_ERR_NORADII;
01569   if (restrictsel && restrictsel->num_atoms != sel->num_atoms)
01570     return MEASURE_ERR_MISMATCHEDCNT;
01571 
01572   int i;
01573   int npts = nsamples ? *nsamples : NPTS;
01574   float maxrad = -1;
01575 
01576 #if 0
01577   // Query min/max atom radii for the entire molecule
01578   mol->get_radii_minmax(minrad, maxrad);
01579 #endif
01580 
01581   // find biggest atom radius 
01582   for (i=0; i<sel->num_atoms; i++) {
01583     float rad = radius[i];
01584     if (maxrad < rad) maxrad = rad;
01585   }
01586 
01587   // Find atoms within maxrad of each other.  
01588   // build a list of pairs for each atom
01589   ResizeArray<int> *pairlist = new ResizeArray<int>[sel->num_atoms];
01590   {
01591     GridSearchPair *pairs;
01592     pairs = vmd_gridsearch1(framepos, sel->num_atoms, sel->on, 
01593                             2.0f * (maxrad + srad), 0, sel->num_atoms * 1000);
01594 
01595     GridSearchPair *p, *tmp; 
01596     for (p = pairs; p != NULL; p = tmp) {
01597       int ind1=p->ind1;
01598       int ind2=p->ind2;
01599       pairlist[ind1].append(ind2);
01600       pairlist[ind2].append(ind1);
01601       tmp = p->next;
01602       free(p);
01603     }
01604   }
01605 
01606   static const float RAND_MAX_INV = 1.0f/float(VMD_RAND_MAX);
01607   // Seed the random number generator before each calculation.  This gives
01608   // reproducible results and still allows a more accurate answer to be
01609   // obtained by increasing the samples size.  I don't know if this is a
01610   // "good" seed value or not, I just picked something random-looking.
01611   vmd_srandom(38572111);
01612 
01613   // All the spheres use the same random points.  
01614   float *spherepts = new float[3L*npts];
01615   for (i=0; i<npts; i++) {
01616     float u1 = (float) vmd_random();
01617     float u2 = (float) vmd_random();
01618     float z = 2.0f*u1*RAND_MAX_INV -1.0f;
01619     float phi = (float) (2.0f*VMD_PI*u2*RAND_MAX_INV);
01620     float R = sqrtf(1.0f-z*z);
01621     float sphi, cphi;
01622     sincosf(phi, &sphi, &cphi);
01623     spherepts[3L*i  ] = R*cphi;
01624     spherepts[3L*i+1] = R*sphi;
01625     spherepts[3L*i+2] = z;
01626   }
01627 
01628   const float prefac = (float) (4 * VMD_PI / npts);
01629   float totarea = 0.0f;
01630   // compute area for each atom based on its pairlist
01631   for (i=sel->firstsel; i<=sel->lastsel; i++) {
01632     if (sel->on[i]) {
01633       // only atoms in restrictsel contribute
01634       if (restrictsel && !restrictsel->on[i]) continue;
01635       const float *loc = framepos+3L*i;
01636       float rad = radius[i]+srad;
01637       float surfpos[3];
01638       int surfpts = npts;
01639       const ResizeArray<int> &nbrs = pairlist[i];
01640       for (int j=0; j<npts; j++) {
01641         surfpos[0] = loc[0] + rad*spherepts[3L*j  ];
01642         surfpos[1] = loc[1] + rad*spherepts[3L*j+1];
01643         surfpos[2] = loc[2] + rad*spherepts[3L*j+2];
01644         int on = 1;
01645         for (int k=0; k<nbrs.num(); k++) {
01646           int ind = nbrs[k];
01647           const float *nbrloc = framepos+3L*ind;
01648           float radsq = radius[ind]+srad; radsq *= radsq;
01649           float dx = surfpos[0]-nbrloc[0];
01650           float dy = surfpos[1]-nbrloc[1];
01651           float dz = surfpos[2]-nbrloc[2];
01652           if (dx*dx + dy*dy + dz*dz <= radsq) {
01653             on = 0;
01654             break;
01655           }
01656         }
01657         if (on) {
01658           if (sasapts) {
01659             sasapts->append3(&surfpos[0]);
01660           }
01661         } else {
01662           surfpts--;
01663         }
01664       }
01665       float atomarea = prefac * rad * rad * surfpts;
01666       totarea += atomarea;
01667     }
01668   }
01669 
01670   delete [] pairlist;
01671   delete [] spherepts;
01672   *sasa = totarea;
01673   return 0;
01674 }
01675 
01676 
01677 
01678 #if 1
01679 
01680 // #define DEBUGSASATHR 1
01681 
01682 typedef struct {
01683   MoleculeList *mlist;
01684   const AtomSel **sellist; 
01685   int numsels;
01686   float srad; 
01687   float *sasalist;
01688   int nsamples;
01689   float *spherepts;
01690 } sasathreadparms;
01691 
01692 // For different values of the random seed, the computed SASA's of brH.pdb 
01693 // converge to within 1% of each other when the number of points is about
01694 // 500.  We therefore use 500 as the default number.
01695 static void * measure_sasa_thread(void *voidparms) {
01696   int threadid;
01697   sasathreadparms *parms = NULL;
01698   wkf_threadlaunch_getdata(voidparms, (void **) &parms);
01699   wkf_threadlaunch_getid(voidparms, &threadid, NULL);
01700 #if defined(DEBUGSASATHR)
01701 printf("sasathread[%d] running...\n", threadid);
01702 #endif
01703 
01704   /*
01705    * copy in per-thread parameters
01706    */
01707   MoleculeList *mlist = parms->mlist;
01708   const AtomSel **sellist = parms->sellist;
01709 //  int numsels = parms->numsels;
01710   float srad = parms->srad;
01711   float *sasalist = parms->sasalist;
01712   const int npts = parms->nsamples;
01713   float *spherepts = parms->spherepts;
01714 
01715   int i, selidx;
01716   wkf_tasktile_t tile;
01717   while (wkf_threadlaunch_next_tile(voidparms, 1, &tile) != WKF_SCHED_DONE) {
01718 #if defined(DEBUGSASATHR)
01719 printf("sasathread[%d] running idx %d to %d\n", threadid, tile.start, tile.end);
01720 #endif
01721     for (selidx=tile.start; selidx<tile.end; selidx++) {
01722       const AtomSel *sel = sellist[selidx];
01723       Molecule *mol = mlist->mol_from_id(sel->molid());
01724 
01725       if (!sel->selected) {
01726         sasalist[selidx] = 0;
01727         continue;
01728       }
01729 
01730       const float *framepos = sel->coordinates(mlist);
01731       if (!framepos) {
01732 #if defined(DEBUGSASATHR)
01733 printf("measure_sasalist: failed to get coords!!!\n");
01734 #endif
01735         return NULL; // MEASURE_ERR_NOFRAMEPOS;
01736       }
01737 
01738       const float *radius = mol->extraflt.data("radius");
01739       if (!radius) {
01740 #if defined(DEBUGSASATHR)
01741 printf("measure_sasalist: failed to get radii!!!\n");
01742 #endif
01743         return NULL; // MEASURE_ERR_NORADII;
01744       }
01745 
01746 
01747       float minrad=-1, maxrad=-1;
01748 #if 1
01749       // Query min/max atom radii for the entire molecule
01750       mol->get_radii_minmax(minrad, maxrad);
01751 #else
01752       // find biggest atom radius 
01753       for (i=0; i<sel->num_atoms; i++) {
01754         float rad = radius[i];
01755         if (maxrad < rad) maxrad = rad;
01756       }
01757 #endif
01758 
01759       // Find atoms within maxrad of each other.  
01760       // build a list of pairs for each atom
01761       ResizeArray<int> *pairlist = new ResizeArray<int>[sel->num_atoms];
01762       {
01763         GridSearchPair *pairs;
01764         pairs = vmd_gridsearch1(framepos, sel->num_atoms, sel->on, 
01765                                 2.0f * (maxrad + srad), 0, 
01766                                 sel->num_atoms * 1000);
01767 
01768         GridSearchPair *p, *tmp; 
01769         for (p = pairs; p != NULL; p = tmp) {
01770           int ind1=p->ind1;
01771           int ind2=p->ind2;
01772           pairlist[ind1].append(ind2);
01773           pairlist[ind2].append(ind1);
01774           tmp = p->next;
01775           free(p);
01776         }
01777       }
01778 
01779       const float prefac = (float) (4 * VMD_PI / npts);
01780       float totarea = 0.0f;
01781       // compute area for each atom based on its pairlist
01782       for (i=sel->firstsel; i<=sel->lastsel; i++) {
01783         if (sel->on[i]) {
01784           const float *loc = framepos+3L*i;
01785           float rad = radius[i]+srad;
01786           float surfpos[3];
01787           int surfpts = npts;
01788           const ResizeArray<int> &nbrs = pairlist[i];
01789           for (int j=0; j<npts; j++) {
01790             surfpos[0] = loc[0] + rad*spherepts[3L*j  ];
01791             surfpos[1] = loc[1] + rad*spherepts[3L*j+1];
01792             surfpos[2] = loc[2] + rad*spherepts[3L*j+2];
01793             int on = 1;
01794             for (int k=0; k<nbrs.num(); k++) {
01795               int ind = nbrs[k];
01796               const float *nbrloc = framepos+3L*ind;
01797               float radsq = radius[ind]+srad; radsq *= radsq;
01798               float dx = surfpos[0]-nbrloc[0];
01799               float dy = surfpos[1]-nbrloc[1];
01800               float dz = surfpos[2]-nbrloc[2];
01801               if (dx*dx + dy*dy + dz*dz <= radsq) {
01802                 on = 0;
01803                 break;
01804               }
01805             }
01806             if (!on) {
01807               surfpts--;
01808             }
01809           }
01810           float atomarea = prefac * rad * rad * surfpts;
01811           totarea += atomarea;
01812         }
01813       }
01814 
01815       delete [] pairlist;
01816       sasalist[selidx] = totarea;
01817     }
01818   }
01819 
01820   return NULL;
01821 }
01822 
01823 #if 1
01824 
01825 // For different values of the random seed, the computed SASA's of brH.pdb 
01826 // converge to within 1% of each other when the number of points is about
01827 // 500.  We therefore use 500 as the default number.
01828 extern int measure_sasalist(MoleculeList *mlist,
01829                             const AtomSel **sellist, int numsels,
01830                             float srad, float *sasalist, const int *nsamples) {
01831 
01832   // check arguments
01833   if (!sellist) return MEASURE_ERR_NOSEL;
01834 
01835   int i, rc;
01836   int npts = nsamples ? *nsamples : NPTS;
01837 
01838 #if defined(VMDTHREADS)
01839   int numprocs = wkf_thread_numprocessors();
01840 #else
01841   int numprocs = 1;
01842 #endif
01843 
01844 #if defined(DEBUGSASATHR)
01845 printf("sasaprocs: %d\n", numprocs);
01846 #endif
01847 
01848   static const float RAND_MAX_INV = 1.0f/float(VMD_RAND_MAX);
01849 
01850   // Seed the random number generator before each calculation.  This gives
01851   // reproducible results and still allows a more accurate answer to be
01852   // obtained by increasing the samples size.  I don't know if this is a
01853   // "good" seed value or not, I just picked something random-looking.
01854   vmd_srandom(38572111);
01855 
01856   // All the spheres use the same random points.  
01857   float *spherepts = new float[3L*npts];
01858   for (i=0; i<npts; i++) {
01859     float u1 = (float) vmd_random();
01860     float u2 = (float) vmd_random();
01861     float z = 2.0f*u1*RAND_MAX_INV -1.0f;
01862     float phi = (float) (2.0f*VMD_PI*u2*RAND_MAX_INV);
01863     float R = sqrtf(1.0f-z*z);
01864     float sphi, cphi;
01865     sincosf(phi, &sphi, &cphi);
01866     spherepts[3L*i  ] = R*cphi;
01867     spherepts[3L*i+1] = R*sphi;
01868     spherepts[3L*i+2] = z;
01869   }
01870 
01871   sasathreadparms parms;
01872   parms.mlist = mlist;
01873   parms.sellist = sellist;
01874   parms.numsels = numsels;
01875   parms.srad = srad;
01876   parms.sasalist = sasalist;
01877   parms.nsamples = npts;
01878   parms.spherepts = spherepts;
01879 
01880 
01881   /* spawn child threads to do the work */
01882   wkf_tasktile_t tile;
01883   tile.start=0;
01884   tile.end=numsels;
01885   rc = wkf_threadlaunch(numprocs, &parms, measure_sasa_thread, &tile);
01886 
01887   delete [] spherepts;
01888 
01889   return rc;
01890 }
01891 
01892 
01893 #else
01894 
01895 // For different values of the random seed, the computed SASA's of brH.pdb 
01896 // converge to within 1% of each other when the number of points is about
01897 // 500.  We therefore use 500 as the default number.
01898 extern int measure_sasalist(MoleculeList *mlist,
01899                             const AtomSel **sellist, int numsels,
01900                             float srad, float *sasalist, const int *nsamples) {
01901 
01902   // check arguments
01903   if (!sellist) return MEASURE_ERR_NOSEL;
01904 
01905   int i;
01906   int npts = nsamples ? *nsamples : NPTS;
01907 
01908   int selidx;
01909   for (selidx=0; selidx<numsels; selidx++) {
01910     const AtomSel *sel = sellist[selidx];
01911     Molecule *mol = mlist->mol_from_id(sel->molid());
01912 
01913     if (!sel->selected) {
01914       sasalist[selidx] = 0;
01915       continue;
01916     }
01917 
01918     const float *framepos = sel->coordinates(mlist);
01919     if (!framepos) {
01920 #if defined(DEBUGSASATHR)
01921 printf("measure_sasalist: failed to get coords!!!\n");
01922 #endif
01923       return MEASURE_ERR_NOFRAMEPOS;
01924     }
01925 
01926     const float *radius = mol->extraflt.data("radius");
01927     if (!radius) {
01928 #if defined(DEBUGSASATHR)
01929 printf("measure_sasalist: failed to get radii!!!\n");
01930 #endif
01931       return MEASURE_ERR_NORADII;
01932     }
01933 
01934     float minrad=-1, maxrad=-1;
01935 #if 1
01936     // Query min/max atom radii for the entire molecule
01937     mol->get_radii_minmax(minrad, maxrad);
01938 #else
01939     // find biggest atom radius 
01940     for (i=0; i<sel->num_atoms; i++) {
01941       float rad = radius[i];
01942       if (maxrad < rad) maxrad = rad;
01943     }
01944 #endif
01945 
01946     // Find atoms within maxrad of each other.  
01947     // build a list of pairs for each atom
01948     ResizeArray<int> *pairlist = new ResizeArray<int>[sel->num_atoms];
01949     {
01950       GridSearchPair *pairs;
01951       pairs = vmd_gridsearch1(framepos, sel->num_atoms, sel->on, 
01952                               2.0f * (maxrad + srad), 0, sel->num_atoms * 1000);
01953 
01954       GridSearchPair *p, *tmp; 
01955       for (p = pairs; p != NULL; p = tmp) {
01956         int ind1=p->ind1;
01957         int ind2=p->ind2;
01958         pairlist[ind1].append(ind2);
01959         pairlist[ind2].append(ind1);
01960         tmp = p->next;
01961         free(p);
01962       }
01963     }
01964 
01965     static const float RAND_MAX_INV = 1.0f/VMD_RAND_MAX;
01966     // Seed the random number generator before each calculation.  This gives
01967     // reproducible results and still allows a more accurate answer to be
01968     // obtained by increasing the samples size.  I don't know if this is a
01969     // "good" seed value or not, I just picked something random-looking.
01970     vmd_srandom(38572111);
01971 
01972     // All the spheres use the same random points.  
01973     float *spherepts = new float[3L*npts];
01974     for (i=0; i<npts; i++) {
01975       float u1 = (float) vmd_random();
01976       float u2 = (float) vmd_random();
01977       float z = 2.0f*u1*RAND_MAX_INV -1.0f;
01978       float phi = (float) (2.0f*VMD_PI*u2*RAND_MAX_INV);
01979       float R = sqrtf(1.0f-z*z);
01980       float sphi, cphi;
01981       sincosf(phi, &sphi, &cphi);
01982       spherepts[3L*i  ] = R*cphi;
01983       spherepts[3L*i+1] = R*sphi;
01984       spherepts[3L*i+2] = z;
01985     }
01986 
01987     const float prefac = (float) (4 * VMD_PI / npts);
01988     float totarea = 0.0f;
01989     // compute area for each atom based on its pairlist
01990     for (i=sel->firstsel; i<=sel->lastsel; i++) {
01991       if (sel->on[i]) {
01992         const float *loc = framepos+3L*i;
01993         float rad = radius[i]+srad;
01994         float surfpos[3];
01995         int surfpts = npts;
01996         const ResizeArray<int> &nbrs = pairlist[i];
01997         for (int j=0; j<npts; j++) {
01998           surfpos[0] = loc[0] + rad*spherepts[3L*j  ];
01999           surfpos[1] = loc[1] + rad*spherepts[3L*j+1];
02000           surfpos[2] = loc[2] + rad*spherepts[3L*j+2];
02001           int on = 1;
02002           for (int k=0; k<nbrs.num(); k++) {
02003             int ind = nbrs[k];
02004             const float *nbrloc = framepos+3L*ind;
02005             float radsq = radius[ind]+srad; radsq *= radsq;
02006             float dx = surfpos[0]-nbrloc[0];
02007             float dy = surfpos[1]-nbrloc[1];
02008             float dz = surfpos[2]-nbrloc[2];
02009             if (dx*dx + dy*dy + dz*dz <= radsq) {
02010               on = 0;
02011               break;
02012             }
02013           }
02014           if (!on) {
02015             surfpts--;
02016           }
02017         }
02018         float atomarea = prefac * rad * rad * surfpts;
02019         totarea += atomarea;
02020       }
02021     }
02022 
02023     delete [] pairlist;
02024     delete [] spherepts;
02025     sasalist[selidx] = totarea;
02026   }
02027 
02028   return 0;
02029 }
02030 
02031 #endif
02032 #endif
02033 
02034 
02035 
02036 //
02037 // Calculate g(r)
02038 //
02039 
02040 // find the minimum image distance for one coordinate component 
02041 // and square the result (orthogonal cells only).
02042 static float fix_pbc_n_sqr(float delta, const float boxby2) {
02043   while (delta >= boxby2) { delta -= 2.0f * boxby2; }
02044   while (delta < -boxby2) { delta += 2.0f * boxby2; }
02045   return delta * delta;
02046 }
02047 
02048 // calculate the minimum distance between two positions with respect 
02049 // to the periodic cell (orthogonal cells only).
02050 static float min_dist_with_pbc(const float *a, const float *b, 
02051                                 const float *boxby2) {
02052   float distsqr;
02053   distsqr  = fix_pbc_n_sqr(a[0] - b[0], boxby2[0]);
02054   distsqr += fix_pbc_n_sqr(a[1] - b[1], boxby2[1]);
02055   distsqr += fix_pbc_n_sqr(a[2] - b[2], boxby2[2]);
02056   return sqrtf(distsqr);
02057 }
02058 
02063 static inline double spherical_cap(const double &radius, const double &boxby2) {
02064   return (VMD_PI / 3.0 * (radius - boxby2) * (radius - boxby2)
02065           * ( 2.0 * radius + boxby2));
02066 }
02067 
02068 
02069 typedef struct {
02070   int threadid;
02071   int threadcount;
02072   int count_o_start;
02073   int count_o_end;
02074   const float *olist;
02075   int count_i;
02076   const float *ilist;
02077   int count_h;
02078   int *hlist;
02079   float delta;
02080   const float *boxby2;
02081   wkfmsgtimer *msgtp;
02082   int curframe;
02083   int maxframe;
02084 } gofrparms_t;
02085     
02086 // calculate the non-normalized pair-distribution function
02087 // for two lists of atom coordinates and store the resulting
02088 // histogram in the hlist array. orthogonal cell version.
02089 //
02090 // NOTE: this is just the workhorse function. special issues related
02091 // to atoms present in both lists have to be dealt with in the uplevel
02092 // functions, that then also can do various post-processing steps.
02093 extern "C" void * measure_gofr_orth(void *voidparms) {
02094   // handle per-thread parameters
02095   gofrparms_t *parms = (gofrparms_t *) voidparms;
02096   const int count_o_start = parms->count_o_start;
02097   const int count_o_end   = parms->count_o_end;
02098   const int count_i       = parms->count_i;
02099   const int count_h       = parms->count_h;
02100   const float *olist      = parms->olist;
02101   const float *ilist      = parms->ilist;
02102   const float *boxby2     = parms->boxby2;
02103   wkfmsgtimer *msgtp         = parms->msgtp;
02104   const int curframe      = parms->curframe;
02105   const int maxframe      = parms->maxframe;
02106   int *hlist              = parms->hlist;
02107   // other local variables
02108   int i, j, idx;
02109   float dist;
02110   const float deltascale = 1.0f / parms->delta;
02111   int msgcnt=0;
02112 
02113   // loop over the chunk of pairs that was associated to this thread.
02114   for (i=count_o_start; i<count_o_end; ++i) {
02115 
02116     // print progress messages only for thread(s) that have
02117     // a timer defined (usually only tid==0).
02118     if (msgtp && wkf_msg_timer_timeout(msgtp)) {
02119       char tmpbuf[1024];
02120       if (msgcnt==0) 
02121         sprintf(tmpbuf, "timestep %d of %d", curframe, maxframe);
02122       else
02123         sprintf(tmpbuf, "timestep %d of %d: (%6.2f %% complete)", curframe, maxframe, 
02124                 (100.0f * i) / (float) (count_o_end - count_o_start + 1));
02125       msgInfo << "vmd_measure_gofr_orth: " << tmpbuf << sendmsg;
02126       msgcnt++;
02127       // XXX we should update the display here...
02128     }
02129 
02130     for (j=0; j<count_i; ++j) {
02131       // calculate distance and add to histogram
02132       dist = min_dist_with_pbc(&olist[i*3L], &ilist[j*3L], boxby2);
02133       idx = (int) (dist * deltascale);
02134       if ((idx >= 0) && (idx < count_h)) 
02135         ++hlist[idx];
02136     }
02137   }
02138 
02139   return MEASURE_NOERR;
02140 }
02141 
02142 // main entry point for 'measure gofr'.
02143 // tasks:
02144 // - sanity check on arguments.
02145 // - select proper algorithm for PBC treatment.
02146 int measure_gofr(AtomSel *sel1, AtomSel *sel2, MoleculeList *mlist,
02147                  const int count_h, double *gofr, double *numint, double *histog,
02148                  const float delta, int first, int last, int step, int *framecntr,
02149                  int usepbc, int selupdate) {
02150   int i, j, frame;
02151   float a, b, c, alpha, beta, gamma;
02152   int isortho=0;    // orthogonal unit cell not assumed by default.
02153   int duplicates=0; // counter for duplicate atoms in both selections.
02154 
02155   // initialize a/b/c/alpha/beta/gamma to arbitrary defaults to please the compiler.
02156   a=b=c=9999999.0;
02157   alpha=beta=gamma=90.0;
02158 
02159   // reset counter for total, skipped, and _orth processed frames.
02160   framecntr[0]=framecntr[1]=framecntr[2]=0;
02161 
02162   // First round of sanity checks.
02163   // neither list can be undefined
02164   if (!sel1 || !sel2 ) {
02165     return MEASURE_ERR_NOSEL;
02166   }
02167 
02168   // make sure that both selections are from the same molecule
02169   // so that we know that PBC unit cell info is the same for both
02170   if (sel2->molid() != sel1->molid()) {
02171     return MEASURE_ERR_MISMATCHEDMOLS;
02172   }
02173 
02174   Molecule *mymol = mlist->mol_from_id(sel1->molid());
02175   int maxframe = mymol->numframes() - 1;
02176   int nframes = 0;
02177 
02178   if (last == -1)
02179     last = maxframe;
02180 
02181   if ((last < first) || (last < 0) || (step <=0) || (first < -1)
02182       || (maxframe < 0) || (last > maxframe)) {
02183       msgErr << "measure gofr: bad frame range given." 
02184              << " max. allowed frame#: " << maxframe << sendmsg;
02185     return MEASURE_ERR_BADFRAMERANGE;
02186   }
02187 
02188   // test for non-orthogonal PBC cells, zero volume cells, etc.
02189   if (usepbc) { 
02190     for (isortho=1, nframes=0, frame=first; frame <=last; ++nframes, frame += step) {
02191       const Timestep *ts;
02192 
02193       if (first == -1) {
02194         // use current frame only. don't loop.
02195         ts = sel1->timestep(mlist);
02196         frame=last;
02197       } else {
02198         ts = mymol->get_frame(frame);
02199       }
02200       // get periodic cell information for current frame
02201       a = ts->a_length;
02202       b = ts->b_length;
02203       c = ts->c_length;
02204       alpha = ts->alpha;
02205       beta = ts->beta;
02206       gamma = ts->gamma;
02207 
02208       // check validity of PBC cell side lengths
02209       if (fabsf(a*b*c) < 0.0001) {
02210         msgErr << "measure gofr: unit cell volume is zero." << sendmsg;
02211         return MEASURE_ERR_GENERAL;
02212       }
02213 
02214       // check PBC unit cell shape to select proper low level algorithm.
02215       if ((alpha != 90.0) || (beta != 90.0) || (gamma != 90.0)) {
02216         isortho=0;
02217       }
02218     }
02219   } else {
02220     // initialize a/b/c/alpha/beta/gamma to arbitrary defaults
02221     isortho=1;
02222     a=b=c=9999999.0;
02223     alpha=beta=gamma=90.0;
02224   }
02225 
02226   // until we can handle non-orthogonal periodic cells, this is fatal
02227   if (!isortho) {
02228     msgErr << "measure gofr: only orthorhombic cells are supported (for now)." << sendmsg;
02229     return MEASURE_ERR_GENERAL;
02230   }
02231 
02232   // clear the result arrays
02233   for (i=0; i<count_h; ++i) {
02234     gofr[i] = numint[i] = histog[i] = 0.0;
02235   }
02236 
02237   // pre-allocate coordinate buffers of the max size we'll
02238   // ever need, so we don't have to reallocate if/when atom
02239   // selections are updated on-the-fly
02240   float *sel1coords = new float[3L*sel1->num_atoms];
02241   float *sel2coords = new float[3L*sel2->num_atoms];
02242 
02243   // setup status message timer
02244   wkfmsgtimer *msgt = wkf_msg_timer_create(5);
02245 
02246   // threading setup.
02247   wkf_thread_t *threads;
02248   gofrparms_t  *parms;
02249 #if defined(VMDTHREADS)
02250   int numprocs = wkf_thread_numprocessors();
02251 #else
02252   int numprocs = 1;
02253 #endif
02254 
02255   threads = new wkf_thread_t[numprocs];
02256   memset(threads, 0, numprocs * sizeof(wkf_thread_t));
02257 
02258   // allocate and (partially) initialize array of per-thread parameters
02259   parms = new gofrparms_t[numprocs];
02260   for (i=0; i<numprocs; ++i) {
02261     parms[i].threadid = i;
02262     parms[i].threadcount = numprocs;
02263     parms[i].delta = (float) delta;
02264     parms[i].msgtp = NULL;
02265     parms[i].count_h = count_h;
02266     parms[i].hlist = new int[count_h];
02267   }
02268 
02269   msgInfo << "measure gofr: using multi-threaded implementation with " 
02270           << numprocs << ((numprocs > 1) ? " CPUs" : " CPU") << sendmsg;
02271 
02272   for (nframes=0,frame=first; frame <=last; ++nframes, frame += step) {
02273     const Timestep *ts1, *ts2;
02274 
02275     if (frame  == -1) {
02276       // use current frame only. don't loop.
02277       ts1 = sel1->timestep(mlist);
02278       ts2 = sel2->timestep(mlist);
02279       frame=last;
02280     } else {
02281       sel1->which_frame = frame;
02282       sel2->which_frame = frame;
02283       ts1 = ts2 = mymol->get_frame(frame); // requires sels from same mol
02284     }
02285 
02286     if (usepbc) {
02287       // get periodic cell information for current frame
02288       a     = ts1->a_length;
02289       b     = ts1->b_length;
02290       c     = ts1->c_length;
02291       alpha = ts1->alpha;
02292       beta  = ts1->beta;
02293       gamma = ts1->gamma;
02294     }
02295 
02296     // compute half periodic cell size
02297     float boxby2[3];
02298     boxby2[0] = 0.5f * a;
02299     boxby2[1] = 0.5f * b;
02300     boxby2[2] = 0.5f * c;
02301 
02302     // update the selections if the user desires it
02303     if (selupdate) {
02304       if (sel1->change(NULL, mymol) != AtomSel::PARSE_SUCCESS)
02305         msgErr << "measure gofr: failed to evaluate atom selection update";
02306       if (sel2->change(NULL, mymol) != AtomSel::PARSE_SUCCESS)
02307         msgErr << "measure gofr: failed to evaluate atom selection update";
02308     }
02309 
02310     // check for duplicate atoms in the two lists, as these will have
02311     // to be subtracted back out of the first histogram slot
02312     if (sel2->molid() == sel1->molid()) {
02313       int i;
02314       for (i=0, duplicates=0; i<sel1->num_atoms; ++i) {
02315         if (sel1->on[i] && sel2->on[i])
02316           ++duplicates;
02317       }
02318     }
02319 
02320     // copy selected atoms to the two coordinate lists
02321     // requires that selections come from the same molecule
02322     const float *framepos = ts1->pos;
02323     for (i=0, j=0; i<sel1->num_atoms; ++i) {
02324       if (sel1->on[i]) {
02325         long a = i*3L;
02326         sel1coords[j    ] = framepos[a    ];
02327         sel1coords[j + 1] = framepos[a + 1];
02328         sel1coords[j + 2] = framepos[a + 2];
02329         j+=3;
02330       }
02331     }
02332     framepos = ts2->pos;
02333     for (i=0, j=0; i<sel2->num_atoms; ++i) {
02334       if (sel2->on[i]) {
02335         long a = i*3L;
02336         sel2coords[j    ] = framepos[a    ];
02337         sel2coords[j + 1] = framepos[a + 1];
02338         sel2coords[j + 2] = framepos[a + 2];
02339         j+=3;
02340       }
02341     }
02342 
02343     // clear the histogram for this frame
02344     // and set up argument structure for threaded execution.
02345     int maxframe = (int) ((last - first + 1) / ((float) step));
02346     for (i=0; i<numprocs; ++i) {
02347       memset(parms[i].hlist, 0, count_h * sizeof(int));
02348       parms[i].boxby2 = boxby2;
02349       parms[i].curframe = frame;
02350       parms[i].maxframe = maxframe;
02351     }
02352     parms[0].msgtp = msgt;
02353     
02354     if (isortho && sel1->selected && sel2->selected) {
02355       int count_o = sel1->selected;
02356       int count_i = sel2->selected;
02357       const float *olist = sel1coords;
02358       const float *ilist = sel2coords;
02359       // make sure the outer loop is the longer one to have 
02360       // better threading efficiency and cache utilization.
02361       if (count_o < count_i) {
02362         count_o = sel2->selected;
02363         count_i = sel1->selected;
02364         olist = sel2coords;
02365         ilist = sel1coords;
02366       }
02367 
02368       // distribute outer loop across threads in fixed size chunks.
02369       // this should work very well for small numbers of threads.
02370       // thrdelta is the chunk size and we need it to be at least 1 
02371       // _and_ numprocs*thrdelta >= count_o.
02372       int thrdelta = (count_o + (numprocs-1)) / numprocs;
02373       int o_min = 0;
02374       int o_max = thrdelta;
02375       for (i=0; i<numprocs; ++i, o_min += thrdelta, o_max += thrdelta) {
02376         if (o_max >  count_o)  o_max = count_o; // curb loop to max
02377         if (o_min >= count_o)  o_max = - 1;     // no work for this thread. too little data.
02378         parms[i].count_o_start = o_min;
02379         parms[i].count_o_end   = o_max;
02380         parms[i].count_i       = count_i;
02381         parms[i].olist         = olist;
02382         parms[i].ilist         = ilist;
02383       }
02384       
02385       // do the gofr calculation for orthogonal boxes.
02386       // XXX. non-orthogonal box not supported yet. detected and handled above.
02387 #if defined(VMDTHREADS)
02388       for (i=0; i<numprocs; ++i) {
02389         wkf_thread_create(&threads[i], measure_gofr_orth, &parms[i]);
02390       }
02391       for (i=0; i<numprocs; ++i) {
02392         wkf_thread_join(threads[i], NULL);
02393       } 
02394 #else
02395       measure_gofr_orth((void *) &parms[0]);
02396 #endif
02397       ++framecntr[2]; // frame processed with _orth algorithm
02398     } else {
02399       ++framecntr[1]; // frame skipped
02400     }
02401     ++framecntr[0];   // total frames.
02402 
02403     // correct the first histogram slot for the number of atoms that are 
02404     // present in both lists. they'll end up in the first histogram bin. 
02405     // we subtract only from the first thread histogram which is always defined.
02406     parms[0].hlist[0] -= duplicates;
02407 
02408     // in case of going 'into the edges', we should cut 
02409     // off the part that is not properly normalized to 
02410     // not confuse people that don't know about this.
02411     int h_max=count_h;
02412     float smallside=a;
02413     if (isortho && usepbc) {
02414       if(b < smallside) {
02415         smallside=b;
02416       }
02417       if(c < smallside) {
02418         smallside=c;
02419       }
02420       h_max=(int) (sqrtf(0.5f)*smallside/delta) +1;
02421       if (h_max > count_h) {
02422         h_max=count_h;
02423       }
02424     }
02425     // compute normalization function.
02426     double all=0.0;
02427     double pair_dens = 0.0;
02428     if (sel1->selected && sel2->selected) {
02429       if (usepbc) {
02430         pair_dens = a * b * c / ((double)sel1->selected * (double)sel2->selected - (double)duplicates);
02431       } else { // assume a particle volume of 30 \AA^3 (~ 1 water).
02432         pair_dens = 30.0 * (double)sel1->selected / 
02433           ((double)sel1->selected * (double)sel2->selected - (double)duplicates);
02434       }
02435     }
02436     
02437     // XXX for orthogonal boxes, we can reduce this to rmax < sqrt(0.5)*smallest side
02438     for (i=0; i<h_max; ++i) {
02439       // radius of inner and outer sphere that form the spherical slice
02440       double r_in  = delta * (double)i;
02441       double r_out = delta * (double)(i+1);
02442       double slice_vol = 4.0 / 3.0 * VMD_PI
02443         * ((r_out * r_out * r_out) - (r_in * r_in * r_in));
02444 
02445       if (isortho && usepbc) {
02446         // add correction for 0.5*box < r <= sqrt(0.5)*box
02447         if (r_out > boxby2[0]) {
02448           slice_vol -= 2.0 * spherical_cap(r_out, boxby2[0]);
02449         }
02450         if (r_out > boxby2[1]) {
02451           slice_vol -= 2.0 * spherical_cap(r_out, boxby2[1]);
02452         }
02453         if (r_out > boxby2[2]) {
02454           slice_vol -= 2.0 * spherical_cap(r_out, boxby2[2]);
02455         }
02456         if (r_in > boxby2[0]) {
02457           slice_vol += 2.0 * spherical_cap(r_in, boxby2[0]);
02458         }
02459         if (r_in > boxby2[1]) {
02460           slice_vol += 2.0 * spherical_cap(r_in, boxby2[1]);
02461         }
02462         if (r_in > boxby2[2]) {
02463           slice_vol += 2.0 * spherical_cap(r_in, boxby2[2]);
02464         }
02465       }
02466 
02467       double normf = pair_dens / slice_vol;
02468       double histv = 0.0;
02469       for (j=0; j<numprocs; ++j) {
02470         histv += (double) parms[j].hlist[i];
02471       }
02472       gofr[i] += normf * histv;
02473       all     += histv;
02474       if (sel1->selected) {
02475         numint[i] += all / (double)(sel1->selected);
02476       }
02477       histog[i] += histv;
02478     }
02479   }
02480   delete [] sel1coords;
02481   delete [] sel2coords;
02482 
02483   double norm = 1.0 / (double) nframes;
02484   for (i=0; i<count_h; ++i) {
02485     gofr[i]   *= norm;
02486     numint[i] *= norm;
02487     histog[i] *= norm;
02488   }
02489   wkf_msg_timer_destroy(msgt);
02490 
02491   // release thread-related storage
02492   for (i=0; i<numprocs; ++i) {
02493     delete [] parms[i].hlist;
02494   }
02495   delete [] threads;
02496   delete [] parms;
02497 
02498   return MEASURE_NOERR;
02499 }
02500 
02501 
02502 int measure_geom(MoleculeList *mlist, int *molid, int *atmid, ResizeArray<float> *gValues,
02503                  int frame, int first, int last, int defmolid, int geomtype) {
02504   int i, ret_val;
02505   for(i=0; i < geomtype; i++) {
02506     // make sure an atom is not repeated in this list
02507     if(i > 0 && molid[i-1]==molid[i] && atmid[i-1]==atmid[i]) {
02508       printf("measure_geom: %i/%i %i/%i\n", molid[i-1],atmid[i-1],molid[i],atmid[i]);
02509       return MEASURE_ERR_REPEATEDATOM;
02510     }
02511   }
02512 
02513   float value;
02514   int max_ts, orig_ts;
02515 
02516   // use the default molecule to determine which frames to cycle through
02517   Molecule *mol = mlist->mol_from_id(defmolid);
02518   if( !mol )
02519     return MEASURE_ERR_NOMOLECULE;
02520   
02521   // get current frame number and make sure there are frames
02522   if((orig_ts = mol->frame()) < 0)
02523     return MEASURE_ERR_NOFRAMES;
02524   
02525   // get the max frame number and determine frame range
02526   max_ts = mol->numframes()-1;
02527   if (frame<0) {
02528     if (first<0 && last<0) first = last = orig_ts; 
02529     if (last<0 || last>max_ts) last = max_ts;
02530     if (first<0) first = 0;
02531   } else {
02532     if (frame>max_ts) frame = max_ts;
02533     first = last = frame; 
02534   }
02535   
02536   // go through all the frames, calculating values
02537   for(i=first; i <= last; i++) {
02538     mol->override_current_frame(i);
02539     switch (geomtype) {
02540     case MEASURE_BOND:
02541       if ((ret_val=calculate_bond(mlist, molid, atmid, &value))<0)
02542         return ret_val;
02543       gValues->append(value);
02544       break;
02545     case MEASURE_ANGLE:
02546       if ((ret_val=calculate_angle(mlist, molid, atmid, &value))<0)
02547         return ret_val;
02548       gValues->append(value);
02549       break;
02550     case MEASURE_DIHED:
02551       if ((ret_val=calculate_dihed(mlist, molid, atmid, &value))<0)
02552         return ret_val;
02553       gValues->append(value);
02554       break;
02555     }
02556   }
02557   
02558   // reset the current frame
02559   mol->override_current_frame(orig_ts);
02560   
02561   return MEASURE_NOERR;
02562 }
02563   
02564   
02565 // calculate the value of this geometry, and return it
02566 int calculate_bond(MoleculeList *mlist, int *molid, int *atmid, float *value) {
02567 
02568   // get coords to calculate distance 
02569   int ret_val;
02570   float pos1[3], pos2[3];
02571   if ((ret_val=normal_atom_coord(mlist->mol_from_id(molid[0]), atmid[0], pos1))<0)
02572     return ret_val;
02573   if ((ret_val=normal_atom_coord(mlist->mol_from_id(molid[1]), atmid[1], pos2))<0)
02574     return ret_val;
02575   
02576   vec_sub(pos2, pos2, pos1);
02577   *value = norm(pos2);
02578 
02579   return MEASURE_NOERR;
02580 }
02581 
02582 // calculate the value of this geometry, and return it
02583 int calculate_angle(MoleculeList *mlist, int *molid, int *atmid, float *value) {
02584 
02585   // get coords to calculate distance 
02586   int ret_val;
02587   float pos1[3], pos2[3], pos3[3], r1[3], r2[3];
02588   if((ret_val=normal_atom_coord(mlist->mol_from_id(molid[0]), atmid[0], pos1))<0)
02589     return ret_val;
02590   if((ret_val=normal_atom_coord(mlist->mol_from_id(molid[1]), atmid[1], pos2))<0)
02591     return ret_val;
02592   if((ret_val=normal_atom_coord(mlist->mol_from_id(molid[2]), atmid[2], pos3))<0)
02593     return ret_val;
02594 
02595   vec_sub(r1, pos1, pos2);
02596   vec_sub(r2, pos3, pos2);
02597   *value = angle(r1, r2);
02598 
02599   return MEASURE_NOERR;
02600 }
02601 
02602 // calculate the value of this geometry, and return it
02603 int calculate_dihed(MoleculeList *mlist, int *molid, int *atmid, float *value) {
02604 
02605   // get coords to calculate distance 
02606   int ret_val;
02607   float pos1[3], pos2[3], pos3[3], pos4[3]; 
02608   if((ret_val=normal_atom_coord(mlist->mol_from_id(molid[0]), atmid[0], pos1))<0)
02609     return ret_val;
02610   if((ret_val=normal_atom_coord(mlist->mol_from_id(molid[1]), atmid[1], pos2))<0)
02611     return ret_val;
02612   if((ret_val=normal_atom_coord(mlist->mol_from_id(molid[2]), atmid[2], pos3))<0)
02613     return ret_val;
02614   if((ret_val=normal_atom_coord(mlist->mol_from_id(molid[3]), atmid[3], pos4))<0)
02615     return ret_val;
02616 
02617   *value = dihedral(pos1, pos2, pos3, pos4);
02618 
02619   return MEASURE_NOERR;
02620 }
02621 
02622 
02623 // for the given Molecule, find the UNTRANSFORMED coords for the given atom
02624 // return Molecule pointer if successful, NULL otherwise.
02625 int normal_atom_coord(Molecule *mol, int a, float *pos) {
02626   Timestep *now;
02627 
02628   int cell[3];
02629   memset(cell, 0, 3L*sizeof(int));
02630 
02631   // get the molecule pointer, and get the coords for the current timestep
02632   int ret_val = check_mol(mol, a);
02633   if (ret_val<0) 
02634     return ret_val;
02635 
02636   if ((now = mol->current())) {
02637     memcpy((void *)pos, (void *)(now->pos + 3L*a), 3L*sizeof(float));
02638     
02639     // Apply periodic image transformation before returning
02640     Matrix4 mat;
02641     now->get_transform_from_cell(cell, mat);
02642     mat.multpoint3d(pos, pos);
02643     
02644     return MEASURE_NOERR;
02645   }
02646   
02647   // if here, error (i.e. molecule contains no frames)
02648   return MEASURE_ERR_NOFRAMES;
02649 }
02650 
02651 
02652 // check whether the given molecule & atom index is OK
02653 // if OK, return Molecule pointer; otherwise, return NULL
02654 int check_mol(Molecule *mol, int a) {
02655 
02656   if (!mol)
02657     return MEASURE_ERR_NOMOLECULE;
02658   if (a < 0 || a >= mol->nAtoms)
02659     return MEASURE_ERR_BADATOMID;
02660   
02661   return MEASURE_NOERR;
02662 }
02663 
02664 
02665 int measure_energy(MoleculeList *mlist, int *molid, int *atmid, int natoms, ResizeArray<float> *gValues,
02666                  int frame, int first, int last, int defmolid, double *params, int geomtype) {
02667   int i, ret_val;
02668   for(i=0; i < natoms; i++) {
02669     // make sure an atom is not repeated in this list
02670     if(i > 0 && molid[i-1]==molid[i] && atmid[i-1]==atmid[i]) {
02671       printf("measure_energy: %i/%i %i/%i\n", molid[i-1],atmid[i-1],molid[i],atmid[i]);
02672       return MEASURE_ERR_REPEATEDATOM;
02673     }
02674   }
02675 
02676   float value;
02677   int max_ts, orig_ts;
02678 
02679   // use the default molecule to determine which frames to cycle through
02680   Molecule *mol = mlist->mol_from_id(defmolid);
02681   if( !mol )
02682     return MEASURE_ERR_NOMOLECULE;
02683   
02684   // get current frame number and make sure there are frames
02685   if((orig_ts = mol->frame()) < 0)
02686     return MEASURE_ERR_NOFRAMES;
02687   
02688   // get the max frame number and determine frame range
02689   max_ts = mol->numframes()-1;
02690   if (frame==-1) {
02691     if (first<0 && last<0) first = last = orig_ts; 
02692     if (last<0 || last>max_ts) last = max_ts;
02693     if (first<0) first = 0;
02694   } else {
02695     if (frame>max_ts || frame==-2) frame = max_ts;
02696     first = last = frame; 
02697   }
02698   
02699   // go through all the frames, calculating values
02700   for(i=first; i <= last; i++) {
02701     mol->override_current_frame(i);
02702     switch (geomtype) {
02703     case MEASURE_BOND:
02704       if ((ret_val=compute_bond_energy(mlist, molid, atmid, &value, (float) params[0], (float) params[1]))<0)
02705         return ret_val;
02706       gValues->append(value);
02707       break;
02708     case MEASURE_ANGLE:
02709       if ((ret_val=compute_angle_energy(mlist, molid, atmid, &value, (float) params[0], (float) params[1], (float) params[2], (float) params[3]))<0)
02710         return ret_val;
02711       gValues->append(value);
02712       break;
02713     case MEASURE_DIHED:
02714       if ((ret_val=compute_dihed_energy(mlist, molid, atmid, &value, (float) params[0], int(params[1]), (float) params[2]))<0)
02715         return ret_val;
02716       gValues->append(value);
02717       break;
02718     case MEASURE_IMPRP:
02719       if ((ret_val=compute_imprp_energy(mlist, molid, atmid, &value, (float) params[0], (float) params[1]))<0)
02720         return ret_val;
02721       gValues->append(value);
02722       break;
02723     case MEASURE_VDW:
02724       if ((ret_val=compute_vdw_energy(mlist, molid, atmid, &value, (float) params[0], (float) params[1], (float) params[2], (float) params[3], (float) params[4], (float) params[5]))<0)
02725         return ret_val;
02726       gValues->append(value);
02727       break;
02728     case MEASURE_ELECT:
02729       if ((ret_val=compute_elect_energy(mlist, molid, atmid, &value, (float) params[0], (float) params[1], (bool) params[2], (bool) params[3], (float) params[4]))<0)
02730         return ret_val;
02731       gValues->append(value);
02732       break;
02733     }
02734   }
02735   
02736   // reset the current frame
02737   mol->override_current_frame(orig_ts);
02738   
02739   return MEASURE_NOERR;
02740 }
02741   
02742 // calculate the energy of this geometry
02743 int compute_bond_energy(MoleculeList *mlist, int *molid, int *atmid, float *energy, float k, float x0) {
02744   int ret_val;
02745   float dist;
02746 
02747   // Get the coordinates
02748   if ((ret_val=calculate_bond(mlist, molid, atmid, &dist))<0)
02749         return ret_val;
02750   float x = dist-x0;
02751   *energy = k*x*x;
02752 
02753   return MEASURE_NOERR;
02754 }
02755 
02756 // calculate the energy of this geometry
02757 int compute_angle_energy(MoleculeList *mlist, int *molid, int *atmid, float *energy,
02758                          float k, float x0, float kub, float s0) {
02759   int ret_val;
02760   float value;
02761 
02762   // Get the coordinates
02763   if ((ret_val=calculate_angle(mlist, molid, atmid, &value))<0)
02764         return ret_val;
02765   float x = (float) DEGTORAD((value-x0));
02766   float s = 0.0f;
02767 
02768   if (kub>0.0f) {
02769     int twoatoms[2];
02770     twoatoms[0] = atmid[0];
02771     twoatoms[1] = atmid[2];
02772     if ((ret_val=calculate_bond(mlist, molid, twoatoms, &value))<0)
02773       return ret_val;
02774     s = value-s0;
02775   }
02776 
02777   *energy = k*x*x + kub*s*s;
02778 
02779   return MEASURE_NOERR;
02780 }
02781 
02782 // calculate the energy of this geometry
02783 int compute_dihed_energy(MoleculeList *mlist, int *molid, int *atmid, float *energy,
02784                          float k, int n, float delta) {
02785   int ret_val;
02786   float value;
02787 
02788   // Get the coordinates
02789   if ((ret_val=calculate_dihed(mlist, molid, atmid, &value))<0)
02790         return ret_val;
02791   *energy = k*(1+cosf((float) (DEGTORAD((n*value-delta)))));
02792 
02793   return MEASURE_NOERR;
02794 }
02795 
02796 // calculate the energy of this geometry
02797 int compute_imprp_energy(MoleculeList *mlist, int *molid, int *atmid, float *energy,
02798                          float k, float x0) {
02799   int ret_val;
02800   float value;
02801 
02802   // Get the coordinates
02803   if ((ret_val=calculate_dihed(mlist, molid, atmid, &value))<0)
02804         return ret_val;
02805   float x = (float) (DEGTORAD((value-x0)));
02806   *energy = k*x*x;
02807 
02808   return MEASURE_NOERR;
02809 }
02810 
02811 // Calculate the VDW energy for specified pair of atoms
02812 // VDW energy:                                               
02813 // Evdw = eps * ((Rmin/dist)**12 - 2*(Rmin/dist)**6)         
02814 // eps = sqrt(eps1*eps2),  Rmin = Rmin1+Rmin2                
02815 int compute_vdw_energy(MoleculeList *mlist, int *molid, int *atmid, float *energy, float eps1, float rmin1,
02816                          float eps2, float rmin2, float cutoff, float switchdist) {
02817   int ret_val;
02818   float dist;
02819 
02820   // Get the coordinates
02821   if ((ret_val=calculate_bond(mlist, molid, atmid, &dist))<0)
02822     return ret_val;
02823 
02824   float sw=1.0;
02825   if (switchdist>0.0 && cutoff>0.0) {
02826     if (dist>=cutoff) {
02827       sw = 0.0;
02828     } else if (dist>=switchdist) {
02829       // This is the CHARMM switching function
02830       float dist2 = dist*dist;
02831       float cut2 = cutoff*cutoff;
02832       float switch2 = switchdist*switchdist;
02833       float s = cut2-dist2;
02834       float range = cut2-switch2;
02835       sw = s*s*(cut2+2*dist2-3*switch2)/(range*range*range);
02836     }
02837   }
02838 
02839   float term6 = (float) powf((rmin1+rmin2)/dist,6);
02840   *energy = sqrtf(eps1*eps2)*(term6*term6 - 2.0f*term6)*sw;
02841 
02842   return MEASURE_NOERR;
02843 }
02844 
02845 int compute_elect_energy(MoleculeList *mlist, int *molid, int *atmid, float *energy, float q1, float q2,
02846                          bool flag1, bool flag2, float cutoff) {
02847   int ret_val;
02848   float dist;
02849 
02850   // Get the coordinates
02851   if ((ret_val=calculate_bond(mlist, molid, atmid, &dist))<0)
02852     return ret_val;
02853 
02854   // Get atom charges
02855   if (!flag1) q1 = mlist->mol_from_id(molid[0])->charge()[atmid[0]];
02856   if (!flag2) q2 = mlist->mol_from_id(molid[0])->charge()[atmid[1]];
02857 
02858   if (cutoff>0.0) {
02859     if (dist<cutoff) {
02860       float efac = 1.0f-dist*dist/(cutoff*cutoff);
02861       *energy = 332.0636f*q1*q2/dist*efac*efac;
02862     } else {
02863       *energy = 0.0f;
02864     }
02865   } else {
02866     *energy = 332.0636f*q1*q2/dist;
02867   }
02868 
02869   return MEASURE_NOERR;
02870 }
02871  
02872 
02873 // Compute the center of mass for a given selection.
02874 // The result is put in rcom which has to have a size of at least 3.
02875 static void center_of_mass(AtomSel *sel, MoleculeList *mlist, float *rcom) {
02876   int i;
02877   float m = 0, mtot = 0;
02878   Molecule *mol = mlist->mol_from_id(sel->molid());
02879 
02880   // get atom masses
02881   const float *mass = mol->mass();
02882 
02883   // get atom coordinates
02884   const float *pos = sel->coordinates(mlist);
02885 
02886   memset(rcom, 0, 3L*sizeof(float));
02887 
02888   // center of mass
02889   for (i=sel->firstsel; i<=sel->lastsel; i++) {
02890     if (sel->on[i]) {
02891       long ind = i * 3L;
02892 
02893       m = mass[i];
02894 
02895       rcom[0] += m*pos[ind    ];
02896       rcom[1] += m*pos[ind + 1];
02897       rcom[2] += m*pos[ind + 2];
02898 
02899       // total mass
02900       mtot += m;
02901     }
02902   }
02903 
02904   rcom[0] /= mtot;
02905   rcom[1] /= mtot;
02906   rcom[2] /= mtot;
02907 }
02908 
02909 
02910 // Calculate principle axes and moments of inertia for selected atoms.
02911 // The corresponding eigenvalues are also returned, they can be used
02912 // to see if two axes are equivalent. The center of mass will be put
02913 // in parameter rcom.
02914 // The user can provide his own set of coordinates in coor. If this
02915 // parameter is NULL then the coordinates from the selection are used.
02916 extern int measure_inertia(AtomSel *sel, MoleculeList *mlist, const float *coor, float rcom[3],
02917                            float priaxes[3][3], float itensor[4][4], float evalue[3]) {
02918   if (!sel)                     return MEASURE_ERR_NOSEL;
02919   if (sel->num_atoms == 0)      return MEASURE_ERR_NOATOMS;
02920 
02921   Molecule *mol = mlist->mol_from_id(sel->molid());
02922 
02923   float x, y, z, m;
02924   float Ixx=0, Iyy=0, Izz=0, Ixy=0, Ixz=0, Iyz=0;
02925   int i,j=0;
02926 
02927   // need to put 3x3 inertia tensor into 4x4 matrix for jacobi eigensolver
02928   // itensor = {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 1}};
02929   memset(itensor, 0, 16L*sizeof(float));
02930   itensor[3][3] = 1.0;
02931 
02932   // compute center of mass
02933   center_of_mass(sel, mlist, rcom);
02934 
02935   // get atom coordinates
02936   const float *pos = sel->coordinates(mlist);
02937 
02938   // get atom masses
02939   const float *mass = mol->mass();
02940 
02941 
02942   // moments of inertia tensor
02943   for (i=sel->firstsel; i<=sel->lastsel; i++) {
02944     if (sel->on[i]) {
02945       // position relative to COM
02946       if (coor) {
02947         // use user provided coordinates
02948         x = coor[j*3L    ] - rcom[0];
02949         y = coor[j*3L + 1] - rcom[1];
02950         z = coor[j*3L + 2] - rcom[2];
02951         j++;
02952       } else {
02953         // use coordinates from selection
02954         x = pos[i*3L    ] - rcom[0];
02955         y = pos[i*3L + 1] - rcom[1];
02956         z = pos[i*3L + 2] - rcom[2];
02957       }
02958 
02959       m = mass[i];
02960 
02961       Ixx += m*(y*y+z*z);
02962       Iyy += m*(x*x+z*z);
02963       Izz += m*(x*x+y*y);
02964       Ixy -= m*x*y;
02965       Ixz -= m*x*z;
02966       Iyz -= m*y*z;
02967     }
02968   }
02969 
02970   itensor[0][0] = Ixx;
02971   itensor[1][1] = Iyy;
02972   itensor[2][2] = Izz;
02973   itensor[0][1] = Ixy;
02974   itensor[1][0] = Ixy;
02975   itensor[0][2] = Ixz;
02976   itensor[2][0] = Ixz;
02977   itensor[1][2] = Iyz;
02978   itensor[2][1] = Iyz;
02979 
02980   // Find the eigenvalues and eigenvectors of moments of inertia tensor.
02981   // The eigenvectors correspond to the principle axes of inertia.
02982   float evector[3][3];
02983   if (jacobi(itensor,evalue,evector) != 0) return MEASURE_ERR_NONZEROJACOBI;
02984 
02985   // transpose the evector matrix to put the vectors in rows
02986   float vectmp;
02987   vectmp=evector[0][1]; evector[0][1]=evector[1][0]; evector[1][0]=vectmp;
02988   vectmp=evector[0][2]; evector[0][2]=evector[2][0]; evector[2][0]=vectmp;
02989   vectmp=evector[2][1]; evector[2][1]=evector[1][2]; evector[1][2]=vectmp;
02990 
02991 
02992   // sort so that the eigenvalues are from largest to smallest
02993   // (or rather so a[0] is eigenvector with largest eigenvalue, ...)
02994   float *a[3];
02995   a[0] = evector[0];
02996   a[1] = evector[1];
02997   a[2] = evector[2];
02998   // The code for SWAP is copied from measure_fit().
02999   // It swaps rows in the eigenvector matrix.
03000 #define SWAP(qq,ww) {                                           \
03001     float v; float *v1;                                         \
03002     v = evalue[qq]; evalue[qq] = evalue[ww]; evalue[ww] = v;    \
03003     v1 = a[qq]; a[qq] = a[ww]; a[ww] = v1;                      \
03004 }
03005   if (evalue[0] < evalue[1]) {
03006     SWAP(0, 1);
03007   }
03008   if (evalue[0] < evalue[2]) {
03009     SWAP(0, 2);
03010   }
03011   if (evalue[1] < evalue[2]) {
03012     SWAP(1, 2);
03013   }
03014 
03015 #if 0
03016   // If the 2nd and 3rd eigenvalues are identical and not close to zero
03017   // then the corresponding axes are not unique. 
03018   if (evalue[1]/evalue[0]>0.1 && fabs(evalue[1]-evalue[2])/evalue[0]<0.05
03019       && fabs(evalue[0]-evalue[1])/evalue[0]>0.05) {
03020     msgInfo << "Principal axes of inertia 2 and 3 are not unique!" << sendmsg;
03021   }
03022 #endif
03023 
03024   for (i=0; i<3; i++) {
03025     for (j=0; j<3; j++) 
03026       priaxes[i][j] = a[i][j];
03027   }
03028 
03029   return MEASURE_NOERR;
03030 }
03031 

Generated on Fri Oct 4 02:44:13 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002