Measure.C

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