00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
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"
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
00039 static const char *measure_error_messages[] = {
00040 "no atom selection",
00041 "no atoms in selection",
00042 "incorrect number of weights for selection",
00043 "internal error: NULL pointer given",
00044 "bad weight sum, would cause divide by zero",
00045 "molecule was deleted(?)",
00046 "cannot understand weight parameter",
00047 "non-number given as parameter",
00048 "two selections don't have the same number of atoms",
00049 "internal error: out of range",
00050 "no coordinates in selection",
00051 "couldn't compute eigenvalue/vectors",
00052 "unknown Tcl problem",
00053 "no atom radii",
00054 "order parameter contains out-of-range atom index",
00055 "order parameter not supported in old RMS fit mode",
00056 "specified frames are out of range",
00057 "both selections must reference the same molecule",
00058 "one atom appears twice in list",
00059 "molecule contains no frames",
00060 "invalid atom id",
00061 "cutoff must be smaller than cell dimension",
00062 "Zero volmap gridsize"
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
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
00080
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
00122
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
00195
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
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
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
00295 int firstsel = 0;
00296 int lastsel = -1;
00297 if (analyze_selection_aligned_dispatch(NULL, num, on, &firstsel, &lastsel, NULL))
00298 return MEASURE_NOERR;
00299
00300
00301 i=firstsel;
00302 if (radii == NULL) {
00303
00304 minx = maxx = framepos[i*3L ];
00305 miny = maxy = framepos[i*3L+1];
00306 minz = maxz = framepos[i*3L+2];
00307 } else {
00308
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
00318 if (radii == NULL) {
00319
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
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
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
00376
00377
00378
00379
00380 int measure_avpos_perresidue(const AtomSel *sel, MoleculeList *mlist,
00381 int start, int end, int step, float *avpos) {
00382
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
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
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
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
00476 const float *framepos = sel->coordinates(mlist);
00477
00478
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
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:
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
00536
00537
00538
00539
00540
00541
00542
00543
00544 if (unitsdebye) {
00545
00546
00547
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
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
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
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
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
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
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], ¢er[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
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
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
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
00735
00736
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
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
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
00770 *rgyr = sqrtf(sum/total_w);
00771
00772 return MEASURE_NOERR;
00773 }
00774
00775
00776
00777
00778
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
00788 if (sel1->selected != sel2->selected) return MEASURE_ERR_MISMATCHEDCNT;
00789
00790
00791
00792
00793 int sel_flg;
00794 if (num == sel1->num_atoms) {
00795 sel_flg = 1;
00796 } else {
00797 sel_flg = 0;
00798 }
00799
00800
00801 float tmp_w;
00802 int w_index = 0;
00803 int sel1ind = sel1->firstsel;
00804 int sel2ind = sel2->firstsel;
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
00818
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
00832 if (sel_flg == 0) {
00833 tmp_w = weight[w_index++];
00834 } else {
00835 tmp_w = weight[sel1ind];
00836 }
00837
00838
00839 rmsdsum += tmp_w * distance2(framepos1 + 3*sel1ind, framepos2 + 3*sel2ind);
00840 wsum += tmp_w;
00841
00842
00843 sel1ind++;
00844 sel2ind++;
00845 }
00846 twsum += wsum;
00847
00848 if (twsum == 0) {
00849 return MEASURE_ERR_BADWEIGHTSUM;
00850 }
00851
00852
00853 rmsd[rescount++] = sqrtf(rmsdsum / wsum);
00854
00855 return rescount;
00856 }
00857
00859
00860
00861
00862
00863
00864
00865
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
00874 if (sel1->selected != sel2->selected) return MEASURE_ERR_MISMATCHEDCNT;
00875
00876
00877
00878
00879 int sel_flg;
00880 if (num == sel1->num_atoms) {
00881 sel_flg = 1;
00882 } else {
00883 sel_flg = 0;
00884 }
00885
00886
00887 float tmp_w;
00888 int w_index = 0;
00889 int sel1ind = sel1->firstsel;
00890 int sel2ind = sel2->firstsel;
00891 float wsum = 0;
00892 float rmsdsum = 0;
00893
00894 *rmsd = 10000000;
00895
00896
00897 int count = sel1->selected;
00898 while (count-- > 0) {
00899
00900
00901 while (!sel1->on[sel1ind])
00902 sel1ind++;
00903 while (!sel2->on[sel2ind])
00904 sel2ind++;
00905
00906
00907 if (sel_flg == 0) {
00908 tmp_w = weight[w_index++];
00909 } else {
00910 tmp_w = weight[sel1ind];
00911 }
00912
00913
00914 rmsdsum += tmp_w * distance2(framepos1 + 3L*sel1ind, framepos2 + 3L*sel2ind);
00915 wsum += tmp_w;
00916
00917
00918 sel1ind++;
00919 sel2ind++;
00920 }
00921
00922
00923 if (wsum == 0) {
00924 return MEASURE_ERR_BADWEIGHTSUM;
00925 }
00926
00927
00928 *rmsd = sqrtf(rmsdsum / wsum);
00929
00930 return MEASURE_NOERR;
00931 }
00932
00933
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;
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;
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
01101
01102
01103
01104
01105
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
01125
01126
01127
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
01145
01146
01147 switch (sel1->selected) {
01148 case 1: {
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
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++;
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
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]];
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
01211
01212
01213
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
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];
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
01254
01255
01256
01257 MatrixFitRMS(num, v1, v2, weight, tmp);
01258
01259 delete [] v1;
01260 delete [] v2;
01261
01262
01263
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
01282 if (order != NULL) {
01283 return MEASURE_ERR_NOTSUP;
01284 }
01285
01286
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
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
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
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
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
01343
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
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
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
01389
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
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
01405 ty.multmatrix(U);
01406 ty.multmatrix(tx);
01407 memcpy(mat->mat, ty.mat, 16L*sizeof(float));
01408 return MEASURE_NOERR;
01409 }
01410
01411
01412
01413
01414 #define NPTS 500
01415 extern int measure_sasa(const AtomSel *sel, const float *framepos,
01416 const float *radius, float srad, float *sasa,
01417 ResizeArray<float> *sasapts, const AtomSel *restrictsel,
01418 const int *nsamples) {
01419
01420
01421 if (!sel) return MEASURE_ERR_NOSEL;
01422 if (!sel->selected) {
01423 *sasa = 0;
01424 return MEASURE_NOERR;
01425 }
01426 if (!framepos) return MEASURE_ERR_NOFRAMEPOS;
01427 if (!radius) return MEASURE_ERR_NORADII;
01428 if (restrictsel && restrictsel->num_atoms != sel->num_atoms)
01429 return MEASURE_ERR_MISMATCHEDCNT;
01430
01431 int i;
01432 int npts = nsamples ? *nsamples : NPTS;
01433 float maxrad = -1;
01434
01435 #if 0
01436
01437 mol->get_radii_minmax(minrad, maxrad);
01438 #endif
01439
01440
01441 for (i=0; i<sel->num_atoms; i++) {
01442 float rad = radius[i];
01443 if (maxrad < rad) maxrad = rad;
01444 }
01445
01446
01447
01448 ResizeArray<int> *pairlist = new ResizeArray<int>[sel->num_atoms];
01449 {
01450 GridSearchPair *pairs;
01451 pairs = vmd_gridsearch1(framepos, sel->num_atoms, sel->on,
01452 2.0f * (maxrad + srad), 0, sel->num_atoms * 1000);
01453
01454 GridSearchPair *p, *tmp;
01455 for (p = pairs; p != NULL; p = tmp) {
01456 int ind1=p->ind1;
01457 int ind2=p->ind2;
01458 pairlist[ind1].append(ind2);
01459 pairlist[ind2].append(ind1);
01460 tmp = p->next;
01461 free(p);
01462 }
01463 }
01464
01465 static const float RAND_MAX_INV = 1.0f/float(VMD_RAND_MAX);
01466
01467
01468
01469
01470 vmd_srandom(38572111);
01471
01472
01473 float *spherepts = new float[3L*npts];
01474 for (i=0; i<npts; i++) {
01475 float u1 = (float) vmd_random();
01476 float u2 = (float) vmd_random();
01477 float z = 2.0f*u1*RAND_MAX_INV -1.0f;
01478 float phi = (float) (2.0f*VMD_PI*u2*RAND_MAX_INV);
01479 float R = sqrtf(1.0f-z*z);
01480 float sphi, cphi;
01481 sincosf(phi, &sphi, &cphi);
01482 spherepts[3L*i ] = R*cphi;
01483 spherepts[3L*i+1] = R*sphi;
01484 spherepts[3L*i+2] = z;
01485 }
01486
01487 const float prefac = (float) (4 * VMD_PI / npts);
01488 float totarea = 0.0f;
01489
01490 for (i=sel->firstsel; i<=sel->lastsel; i++) {
01491 if (sel->on[i]) {
01492
01493 if (restrictsel && !restrictsel->on[i]) continue;
01494 const float *loc = framepos+3L*i;
01495 float rad = radius[i]+srad;
01496 float surfpos[3];
01497 int surfpts = npts;
01498 const ResizeArray<int> &nbrs = pairlist[i];
01499 for (int j=0; j<npts; j++) {
01500 surfpos[0] = loc[0] + rad*spherepts[3L*j ];
01501 surfpos[1] = loc[1] + rad*spherepts[3L*j+1];
01502 surfpos[2] = loc[2] + rad*spherepts[3L*j+2];
01503 int on = 1;
01504 for (int k=0; k<nbrs.num(); k++) {
01505 int ind = nbrs[k];
01506 const float *nbrloc = framepos+3L*ind;
01507 float radsq = radius[ind]+srad; radsq *= radsq;
01508 float dx = surfpos[0]-nbrloc[0];
01509 float dy = surfpos[1]-nbrloc[1];
01510 float dz = surfpos[2]-nbrloc[2];
01511 if (dx*dx + dy*dy + dz*dz <= radsq) {
01512 on = 0;
01513 break;
01514 }
01515 }
01516 if (on) {
01517 if (sasapts) {
01518 sasapts->append3(&surfpos[0]);
01519 }
01520 } else {
01521 surfpts--;
01522 }
01523 }
01524 float atomarea = prefac * rad * rad * surfpts;
01525 totarea += atomarea;
01526 }
01527 }
01528
01529 delete [] pairlist;
01530 delete [] spherepts;
01531 *sasa = totarea;
01532 return 0;
01533 }
01534
01535
01536
01537 #if 1
01538
01539
01540
01541 typedef struct {
01542 MoleculeList *mlist;
01543 const AtomSel **sellist;
01544 int numsels;
01545 float srad;
01546 float *sasalist;
01547 int nsamples;
01548 float *spherepts;
01549 } sasathreadparms;
01550
01551
01552
01553
01554 static void * measure_sasa_thread(void *voidparms) {
01555 int threadid;
01556 sasathreadparms *parms = NULL;
01557 wkf_threadlaunch_getdata(voidparms, (void **) &parms);
01558 wkf_threadlaunch_getid(voidparms, &threadid, NULL);
01559 #if defined(DEBUGSASATHR)
01560 printf("sasathread[%d] running...\n", threadid);
01561 #endif
01562
01563
01564
01565
01566 MoleculeList *mlist = parms->mlist;
01567 const AtomSel **sellist = parms->sellist;
01568
01569 float srad = parms->srad;
01570 float *sasalist = parms->sasalist;
01571 const int npts = parms->nsamples;
01572 float *spherepts = parms->spherepts;
01573
01574 int i, selidx;
01575 wkf_tasktile_t tile;
01576 while (wkf_threadlaunch_next_tile(voidparms, 1, &tile) != WKF_SCHED_DONE) {
01577 #if defined(DEBUGSASATHR)
01578 printf("sasathread[%d] running idx %d to %d\n", threadid, tile.start, tile.end);
01579 #endif
01580 for (selidx=tile.start; selidx<tile.end; selidx++) {
01581 const AtomSel *sel = sellist[selidx];
01582 Molecule *mol = mlist->mol_from_id(sel->molid());
01583
01584 if (!sel->selected) {
01585 sasalist[selidx] = 0;
01586 continue;
01587 }
01588
01589 const float *framepos = sel->coordinates(mlist);
01590 if (!framepos) {
01591 #if defined(DEBUGSASATHR)
01592 printf("measure_sasalist: failed to get coords!!!\n");
01593 #endif
01594 return NULL;
01595 }
01596
01597 const float *radius = mol->extraflt.data("radius");
01598 if (!radius) {
01599 #if defined(DEBUGSASATHR)
01600 printf("measure_sasalist: failed to get radii!!!\n");
01601 #endif
01602 return NULL;
01603 }
01604
01605
01606 float minrad=-1, maxrad=-1;
01607 #if 1
01608
01609 mol->get_radii_minmax(minrad, maxrad);
01610 #else
01611
01612 for (i=0; i<sel->num_atoms; i++) {
01613 float rad = radius[i];
01614 if (maxrad < rad) maxrad = rad;
01615 }
01616 #endif
01617
01618
01619
01620 ResizeArray<int> *pairlist = new ResizeArray<int>[sel->num_atoms];
01621 {
01622 GridSearchPair *pairs;
01623 pairs = vmd_gridsearch1(framepos, sel->num_atoms, sel->on,
01624 2.0f * (maxrad + srad), 0,
01625 sel->num_atoms * 1000);
01626
01627 GridSearchPair *p, *tmp;
01628 for (p = pairs; p != NULL; p = tmp) {
01629 int ind1=p->ind1;
01630 int ind2=p->ind2;
01631 pairlist[ind1].append(ind2);
01632 pairlist[ind2].append(ind1);
01633 tmp = p->next;
01634 free(p);
01635 }
01636 }
01637
01638 const float prefac = (float) (4 * VMD_PI / npts);
01639 float totarea = 0.0f;
01640
01641 for (i=sel->firstsel; i<=sel->lastsel; i++) {
01642 if (sel->on[i]) {
01643 const float *loc = framepos+3L*i;
01644 float rad = radius[i]+srad;
01645 float surfpos[3];
01646 int surfpts = npts;
01647 const ResizeArray<int> &nbrs = pairlist[i];
01648 for (int j=0; j<npts; j++) {
01649 surfpos[0] = loc[0] + rad*spherepts[3L*j ];
01650 surfpos[1] = loc[1] + rad*spherepts[3L*j+1];
01651 surfpos[2] = loc[2] + rad*spherepts[3L*j+2];
01652 int on = 1;
01653 for (int k=0; k<nbrs.num(); k++) {
01654 int ind = nbrs[k];
01655 const float *nbrloc = framepos+3L*ind;
01656 float radsq = radius[ind]+srad; radsq *= radsq;
01657 float dx = surfpos[0]-nbrloc[0];
01658 float dy = surfpos[1]-nbrloc[1];
01659 float dz = surfpos[2]-nbrloc[2];
01660 if (dx*dx + dy*dy + dz*dz <= radsq) {
01661 on = 0;
01662 break;
01663 }
01664 }
01665 if (!on) {
01666 surfpts--;
01667 }
01668 }
01669 float atomarea = prefac * rad * rad * surfpts;
01670 totarea += atomarea;
01671 }
01672 }
01673
01674 delete [] pairlist;
01675 sasalist[selidx] = totarea;
01676 }
01677 }
01678
01679 return NULL;
01680 }
01681
01682 #if 1
01683
01684
01685
01686
01687 extern int measure_sasalist(MoleculeList *mlist,
01688 const AtomSel **sellist, int numsels,
01689 float srad, float *sasalist, const int *nsamples) {
01690
01691
01692 if (!sellist) return MEASURE_ERR_NOSEL;
01693
01694 int i, rc;
01695 int npts = nsamples ? *nsamples : NPTS;
01696
01697 #if defined(VMDTHREADS)
01698 int numprocs = wkf_thread_numprocessors();
01699 #else
01700 int numprocs = 1;
01701 #endif
01702
01703 #if defined(DEBUGSASATHR)
01704 printf("sasaprocs: %d\n", numprocs);
01705 #endif
01706
01707 static const float RAND_MAX_INV = 1.0f/float(VMD_RAND_MAX);
01708
01709
01710
01711
01712
01713 vmd_srandom(38572111);
01714
01715
01716 float *spherepts = new float[3L*npts];
01717 for (i=0; i<npts; i++) {
01718 float u1 = (float) vmd_random();
01719 float u2 = (float) vmd_random();
01720 float z = 2.0f*u1*RAND_MAX_INV -1.0f;
01721 float phi = (float) (2.0f*VMD_PI*u2*RAND_MAX_INV);
01722 float R = sqrtf(1.0f-z*z);
01723 float sphi, cphi;
01724 sincosf(phi, &sphi, &cphi);
01725 spherepts[3L*i ] = R*cphi;
01726 spherepts[3L*i+1] = R*sphi;
01727 spherepts[3L*i+2] = z;
01728 }
01729
01730 sasathreadparms parms;
01731 parms.mlist = mlist;
01732 parms.sellist = sellist;
01733 parms.numsels = numsels;
01734 parms.srad = srad;
01735 parms.sasalist = sasalist;
01736 parms.nsamples = npts;
01737 parms.spherepts = spherepts;
01738
01739
01740
01741 wkf_tasktile_t tile;
01742 tile.start=0;
01743 tile.end=numsels;
01744 rc = wkf_threadlaunch(numprocs, &parms, measure_sasa_thread, &tile);
01745
01746 delete [] spherepts;
01747
01748 return rc;
01749 }
01750
01751
01752 #else
01753
01754
01755
01756
01757 extern int measure_sasalist(MoleculeList *mlist,
01758 const AtomSel **sellist, int numsels,
01759 float srad, float *sasalist, const int *nsamples) {
01760
01761
01762 if (!sellist) return MEASURE_ERR_NOSEL;
01763
01764 int i;
01765 int npts = nsamples ? *nsamples : NPTS;
01766
01767 int selidx;
01768 for (selidx=0; selidx<numsels; selidx++) {
01769 const AtomSel *sel = sellist[selidx];
01770 Molecule *mol = mlist->mol_from_id(sel->molid());
01771
01772 if (!sel->selected) {
01773 sasalist[selidx] = 0;
01774 continue;
01775 }
01776
01777 const float *framepos = sel->coordinates(mlist);
01778 if (!framepos) {
01779 #if defined(DEBUGSASATHR)
01780 printf("measure_sasalist: failed to get coords!!!\n");
01781 #endif
01782 return MEASURE_ERR_NOFRAMEPOS;
01783 }
01784
01785 const float *radius = mol->extraflt.data("radius");
01786 if (!radius) {
01787 #if defined(DEBUGSASATHR)
01788 printf("measure_sasalist: failed to get radii!!!\n");
01789 #endif
01790 return MEASURE_ERR_NORADII;
01791 }
01792
01793 float minrad=-1, maxrad=-1;
01794 #if 1
01795
01796 mol->get_radii_minmax(minrad, maxrad);
01797 #else
01798
01799 for (i=0; i<sel->num_atoms; i++) {
01800 float rad = radius[i];
01801 if (maxrad < rad) maxrad = rad;
01802 }
01803 #endif
01804
01805
01806
01807 ResizeArray<int> *pairlist = new ResizeArray<int>[sel->num_atoms];
01808 {
01809 GridSearchPair *pairs;
01810 pairs = vmd_gridsearch1(framepos, sel->num_atoms, sel->on,
01811 2.0f * (maxrad + srad), 0, sel->num_atoms * 1000);
01812
01813 GridSearchPair *p, *tmp;
01814 for (p = pairs; p != NULL; p = tmp) {
01815 int ind1=p->ind1;
01816 int ind2=p->ind2;
01817 pairlist[ind1].append(ind2);
01818 pairlist[ind2].append(ind1);
01819 tmp = p->next;
01820 free(p);
01821 }
01822 }
01823
01824 static const float RAND_MAX_INV = 1.0f/VMD_RAND_MAX;
01825
01826
01827
01828
01829 vmd_srandom(38572111);
01830
01831
01832 float *spherepts = new float[3L*npts];
01833 for (i=0; i<npts; i++) {
01834 float u1 = (float) vmd_random();
01835 float u2 = (float) vmd_random();
01836 float z = 2.0f*u1*RAND_MAX_INV -1.0f;
01837 float phi = (float) (2.0f*VMD_PI*u2*RAND_MAX_INV);
01838 float R = sqrtf(1.0f-z*z);
01839 float sphi, cphi;
01840 sincosf(phi, &sphi, &cphi);
01841 spherepts[3L*i ] = R*cphi;
01842 spherepts[3L*i+1] = R*sphi;
01843 spherepts[3L*i+2] = z;
01844 }
01845
01846 const float prefac = (float) (4 * VMD_PI / npts);
01847 float totarea = 0.0f;
01848
01849 for (i=sel->firstsel; i<=sel->lastsel; i++) {
01850 if (sel->on[i]) {
01851 const float *loc = framepos+3L*i;
01852 float rad = radius[i]+srad;
01853 float surfpos[3];
01854 int surfpts = npts;
01855 const ResizeArray<int> &nbrs = pairlist[i];
01856 for (int j=0; j<npts; j++) {
01857 surfpos[0] = loc[0] + rad*spherepts[3L*j ];
01858 surfpos[1] = loc[1] + rad*spherepts[3L*j+1];
01859 surfpos[2] = loc[2] + rad*spherepts[3L*j+2];
01860 int on = 1;
01861 for (int k=0; k<nbrs.num(); k++) {
01862 int ind = nbrs[k];
01863 const float *nbrloc = framepos+3L*ind;
01864 float radsq = radius[ind]+srad; radsq *= radsq;
01865 float dx = surfpos[0]-nbrloc[0];
01866 float dy = surfpos[1]-nbrloc[1];
01867 float dz = surfpos[2]-nbrloc[2];
01868 if (dx*dx + dy*dy + dz*dz <= radsq) {
01869 on = 0;
01870 break;
01871 }
01872 }
01873 if (!on) {
01874 surfpts--;
01875 }
01876 }
01877 float atomarea = prefac * rad * rad * surfpts;
01878 totarea += atomarea;
01879 }
01880 }
01881
01882 delete [] pairlist;
01883 delete [] spherepts;
01884 sasalist[selidx] = totarea;
01885 }
01886
01887 return 0;
01888 }
01889
01890 #endif
01891 #endif
01892
01893
01894
01895
01896
01897
01898
01899
01900
01901 static float fix_pbc_n_sqr(float delta, const float boxby2) {
01902 while (delta >= boxby2) { delta -= 2.0f * boxby2; }
01903 while (delta < -boxby2) { delta += 2.0f * boxby2; }
01904 return delta * delta;
01905 }
01906
01907
01908
01909 static float min_dist_with_pbc(const float *a, const float *b,
01910 const float *boxby2) {
01911 float distsqr;
01912 distsqr = fix_pbc_n_sqr(a[0] - b[0], boxby2[0]);
01913 distsqr += fix_pbc_n_sqr(a[1] - b[1], boxby2[1]);
01914 distsqr += fix_pbc_n_sqr(a[2] - b[2], boxby2[2]);
01915 return sqrtf(distsqr);
01916 }
01917
01922 static inline double spherical_cap(const double &radius, const double &boxby2) {
01923 return (VMD_PI / 3.0 * (radius - boxby2) * (radius - boxby2)
01924 * ( 2.0 * radius + boxby2));
01925 }
01926
01927
01928 typedef struct {
01929 int threadid;
01930 int threadcount;
01931 int count_o_start;
01932 int count_o_end;
01933 const float *olist;
01934 int count_i;
01935 const float *ilist;
01936 int count_h;
01937 int *hlist;
01938 float delta;
01939 const float *boxby2;
01940 wkfmsgtimer *msgtp;
01941 int curframe;
01942 int maxframe;
01943 } gofrparms_t;
01944
01945
01946
01947
01948
01949
01950
01951
01952 extern "C" void * measure_gofr_orth(void *voidparms) {
01953
01954 gofrparms_t *parms = (gofrparms_t *) voidparms;
01955 const int count_o_start = parms->count_o_start;
01956 const int count_o_end = parms->count_o_end;
01957 const int count_i = parms->count_i;
01958 const int count_h = parms->count_h;
01959 const float *olist = parms->olist;
01960 const float *ilist = parms->ilist;
01961 const float *boxby2 = parms->boxby2;
01962 wkfmsgtimer *msgtp = parms->msgtp;
01963 const int curframe = parms->curframe;
01964 const int maxframe = parms->maxframe;
01965 int *hlist = parms->hlist;
01966
01967 int i, j, idx;
01968 float dist;
01969 const float deltascale = 1.0f / parms->delta;
01970 int msgcnt=0;
01971
01972
01973 for (i=count_o_start; i<count_o_end; ++i) {
01974
01975
01976
01977 if (msgtp && wkf_msg_timer_timeout(msgtp)) {
01978 char tmpbuf[1024];
01979 if (msgcnt==0)
01980 sprintf(tmpbuf, "timestep %d of %d", curframe, maxframe);
01981 else
01982 sprintf(tmpbuf, "timestep %d of %d: (%6.2f %% complete)", curframe, maxframe,
01983 (100.0f * i) / (float) (count_o_end - count_o_start + 1));
01984 msgInfo << "vmd_measure_gofr_orth: " << tmpbuf << sendmsg;
01985 msgcnt++;
01986
01987 }
01988
01989 for (j=0; j<count_i; ++j) {
01990
01991 dist = min_dist_with_pbc(&olist[i*3L], &ilist[j*3L], boxby2);
01992 idx = (int) (dist * deltascale);
01993 if ((idx >= 0) && (idx < count_h))
01994 ++hlist[idx];
01995 }
01996 }
01997
01998 return MEASURE_NOERR;
01999 }
02000
02001
02002
02003
02004
02005 int measure_gofr(AtomSel *sel1, AtomSel *sel2, MoleculeList *mlist,
02006 const int count_h, double *gofr, double *numint, double *histog,
02007 const float delta, int first, int last, int step, int *framecntr,
02008 int usepbc, int selupdate) {
02009 int i, j, frame;
02010 float a, b, c, alpha, beta, gamma;
02011 int isortho=0;
02012 int duplicates=0;
02013
02014
02015 a=b=c=9999999.0;
02016 alpha=beta=gamma=90.0;
02017
02018
02019 framecntr[0]=framecntr[1]=framecntr[2]=0;
02020
02021
02022
02023 if (!sel1 || !sel2 ) {
02024 return MEASURE_ERR_NOSEL;
02025 }
02026
02027
02028
02029 if (sel2->molid() != sel1->molid()) {
02030 return MEASURE_ERR_MISMATCHEDMOLS;
02031 }
02032
02033 Molecule *mymol = mlist->mol_from_id(sel1->molid());
02034 int maxframe = mymol->numframes() - 1;
02035 int nframes = 0;
02036
02037 if (last == -1)
02038 last = maxframe;
02039
02040 if ((last < first) || (last < 0) || (step <=0) || (first < -1)
02041 || (maxframe < 0) || (last > maxframe)) {
02042 msgErr << "measure gofr: bad frame range given."
02043 << " max. allowed frame#: " << maxframe << sendmsg;
02044 return MEASURE_ERR_BADFRAMERANGE;
02045 }
02046
02047
02048 if (usepbc) {
02049 for (isortho=1, nframes=0, frame=first; frame <=last; ++nframes, frame += step) {
02050 const Timestep *ts;
02051
02052 if (first == -1) {
02053
02054 ts = sel1->timestep(mlist);
02055 frame=last;
02056 } else {
02057 ts = mymol->get_frame(frame);
02058 }
02059
02060 a = ts->a_length;
02061 b = ts->b_length;
02062 c = ts->c_length;
02063 alpha = ts->alpha;
02064 beta = ts->beta;
02065 gamma = ts->gamma;
02066
02067
02068 if (fabsf(a*b*c) < 0.0001) {
02069 msgErr << "measure gofr: unit cell volume is zero." << sendmsg;
02070 return MEASURE_ERR_GENERAL;
02071 }
02072
02073
02074 if ((alpha != 90.0) || (beta != 90.0) || (gamma != 90.0)) {
02075 isortho=0;
02076 }
02077 }
02078 } else {
02079
02080 isortho=1;
02081 a=b=c=9999999.0;
02082 alpha=beta=gamma=90.0;
02083 }
02084
02085
02086 if (!isortho) {
02087 msgErr << "measure gofr: only orthorhombic cells are supported (for now)." << sendmsg;
02088 return MEASURE_ERR_GENERAL;
02089 }
02090
02091
02092 for (i=0; i<count_h; ++i) {
02093 gofr[i] = numint[i] = histog[i] = 0.0;
02094 }
02095
02096
02097
02098
02099 float *sel1coords = new float[3L*sel1->num_atoms];
02100 float *sel2coords = new float[3L*sel2->num_atoms];
02101
02102
02103 wkfmsgtimer *msgt = wkf_msg_timer_create(5);
02104
02105
02106 wkf_thread_t *threads;
02107 gofrparms_t *parms;
02108 #if defined(VMDTHREADS)
02109 int numprocs = wkf_thread_numprocessors();
02110 #else
02111 int numprocs = 1;
02112 #endif
02113
02114 threads = new wkf_thread_t[numprocs];
02115 memset(threads, 0, numprocs * sizeof(wkf_thread_t));
02116
02117
02118 parms = new gofrparms_t[numprocs];
02119 for (i=0; i<numprocs; ++i) {
02120 parms[i].threadid = i;
02121 parms[i].threadcount = numprocs;
02122 parms[i].delta = (float) delta;
02123 parms[i].msgtp = NULL;
02124 parms[i].count_h = count_h;
02125 parms[i].hlist = new int[count_h];
02126 }
02127
02128 msgInfo << "measure gofr: using multi-threaded implementation with "
02129 << numprocs << ((numprocs > 1) ? " CPUs" : " CPU") << sendmsg;
02130
02131 for (nframes=0,frame=first; frame <=last; ++nframes, frame += step) {
02132 const Timestep *ts1, *ts2;
02133
02134 if (frame == -1) {
02135
02136 ts1 = sel1->timestep(mlist);
02137 ts2 = sel2->timestep(mlist);
02138 frame=last;
02139 } else {
02140 sel1->which_frame = frame;
02141 sel2->which_frame = frame;
02142 ts1 = ts2 = mymol->get_frame(frame);
02143 }
02144
02145 if (usepbc) {
02146
02147 a = ts1->a_length;
02148 b = ts1->b_length;
02149 c = ts1->c_length;
02150 alpha = ts1->alpha;
02151 beta = ts1->beta;
02152 gamma = ts1->gamma;
02153 }
02154
02155
02156 float boxby2[3];
02157 boxby2[0] = 0.5f * a;
02158 boxby2[1] = 0.5f * b;
02159 boxby2[2] = 0.5f * c;
02160
02161
02162 if (selupdate) {
02163 if (sel1->change(NULL, mymol) != AtomSel::PARSE_SUCCESS)
02164 msgErr << "measure gofr: failed to evaluate atom selection update";
02165 if (sel2->change(NULL, mymol) != AtomSel::PARSE_SUCCESS)
02166 msgErr << "measure gofr: failed to evaluate atom selection update";
02167 }
02168
02169
02170
02171 if (sel2->molid() == sel1->molid()) {
02172 int i;
02173 for (i=0, duplicates=0; i<sel1->num_atoms; ++i) {
02174 if (sel1->on[i] && sel2->on[i])
02175 ++duplicates;
02176 }
02177 }
02178
02179
02180
02181 const float *framepos = ts1->pos;
02182 for (i=0, j=0; i<sel1->num_atoms; ++i) {
02183 if (sel1->on[i]) {
02184 long a = i*3L;
02185 sel1coords[j ] = framepos[a ];
02186 sel1coords[j + 1] = framepos[a + 1];
02187 sel1coords[j + 2] = framepos[a + 2];
02188 j+=3;
02189 }
02190 }
02191 framepos = ts2->pos;
02192 for (i=0, j=0; i<sel2->num_atoms; ++i) {
02193 if (sel2->on[i]) {
02194 long a = i*3L;
02195 sel2coords[j ] = framepos[a ];
02196 sel2coords[j + 1] = framepos[a + 1];
02197 sel2coords[j + 2] = framepos[a + 2];
02198 j+=3;
02199 }
02200 }
02201
02202
02203
02204 int maxframe = (int) ((last - first + 1) / ((float) step));
02205 for (i=0; i<numprocs; ++i) {
02206 memset(parms[i].hlist, 0, count_h * sizeof(int));
02207 parms[i].boxby2 = boxby2;
02208 parms[i].curframe = frame;
02209 parms[i].maxframe = maxframe;
02210 }
02211 parms[0].msgtp = msgt;
02212
02213 if (isortho && sel1->selected && sel2->selected) {
02214 int count_o = sel1->selected;
02215 int count_i = sel2->selected;
02216 const float *olist = sel1coords;
02217 const float *ilist = sel2coords;
02218
02219
02220 if (count_o < count_i) {
02221 count_o = sel2->selected;
02222 count_i = sel1->selected;
02223 olist = sel2coords;
02224 ilist = sel1coords;
02225 }
02226
02227
02228
02229
02230
02231 int thrdelta = (count_o + (numprocs-1)) / numprocs;
02232 int o_min = 0;
02233 int o_max = thrdelta;
02234 for (i=0; i<numprocs; ++i, o_min += thrdelta, o_max += thrdelta) {
02235 if (o_max > count_o) o_max = count_o;
02236 if (o_min >= count_o) o_max = - 1;
02237 parms[i].count_o_start = o_min;
02238 parms[i].count_o_end = o_max;
02239 parms[i].count_i = count_i;
02240 parms[i].olist = olist;
02241 parms[i].ilist = ilist;
02242 }
02243
02244
02245
02246 #if defined(VMDTHREADS)
02247 for (i=0; i<numprocs; ++i) {
02248 wkf_thread_create(&threads[i], measure_gofr_orth, &parms[i]);
02249 }
02250 for (i=0; i<numprocs; ++i) {
02251 wkf_thread_join(threads[i], NULL);
02252 }
02253 #else
02254 measure_gofr_orth((void *) &parms[0]);
02255 #endif
02256 ++framecntr[2];
02257 } else {
02258 ++framecntr[1];
02259 }
02260 ++framecntr[0];
02261
02262
02263
02264
02265 parms[0].hlist[0] -= duplicates;
02266
02267
02268
02269
02270 int h_max=count_h;
02271 float smallside=a;
02272 if (isortho && usepbc) {
02273 if(b < smallside) {
02274 smallside=b;
02275 }
02276 if(c < smallside) {
02277 smallside=c;
02278 }
02279 h_max=(int) (sqrtf(0.5f)*smallside/delta) +1;
02280 if (h_max > count_h) {
02281 h_max=count_h;
02282 }
02283 }
02284
02285 double all=0.0;
02286 double pair_dens = 0.0;
02287 if (sel1->selected && sel2->selected) {
02288 if (usepbc) {
02289 pair_dens = a * b * c / ((double)sel1->selected * (double)sel2->selected - (double)duplicates);
02290 } else {
02291 pair_dens = 30.0 * (double)sel1->selected /
02292 ((double)sel1->selected * (double)sel2->selected - (double)duplicates);
02293 }
02294 }
02295
02296
02297 for (i=0; i<h_max; ++i) {
02298
02299 double r_in = delta * (double)i;
02300 double r_out = delta * (double)(i+1);
02301 double slice_vol = 4.0 / 3.0 * VMD_PI
02302 * ((r_out * r_out * r_out) - (r_in * r_in * r_in));
02303
02304 if (isortho && usepbc) {
02305
02306 if (r_out > boxby2[0]) {
02307 slice_vol -= 2.0 * spherical_cap(r_out, boxby2[0]);
02308 }
02309 if (r_out > boxby2[1]) {
02310 slice_vol -= 2.0 * spherical_cap(r_out, boxby2[1]);
02311 }
02312 if (r_out > boxby2[2]) {
02313 slice_vol -= 2.0 * spherical_cap(r_out, boxby2[2]);
02314 }
02315 if (r_in > boxby2[0]) {
02316 slice_vol += 2.0 * spherical_cap(r_in, boxby2[0]);
02317 }
02318 if (r_in > boxby2[1]) {
02319 slice_vol += 2.0 * spherical_cap(r_in, boxby2[1]);
02320 }
02321 if (r_in > boxby2[2]) {
02322 slice_vol += 2.0 * spherical_cap(r_in, boxby2[2]);
02323 }
02324 }
02325
02326 double normf = pair_dens / slice_vol;
02327 double histv = 0.0;
02328 for (j=0; j<numprocs; ++j) {
02329 histv += (double) parms[j].hlist[i];
02330 }
02331 gofr[i] += normf * histv;
02332 all += histv;
02333 if (sel1->selected) {
02334 numint[i] += all / (double)(sel1->selected);
02335 }
02336 histog[i] += histv;
02337 }
02338 }
02339 delete [] sel1coords;
02340 delete [] sel2coords;
02341
02342 double norm = 1.0 / (double) nframes;
02343 for (i=0; i<count_h; ++i) {
02344 gofr[i] *= norm;
02345 numint[i] *= norm;
02346 histog[i] *= norm;
02347 }
02348 wkf_msg_timer_destroy(msgt);
02349
02350
02351 for (i=0; i<numprocs; ++i) {
02352 delete [] parms[i].hlist;
02353 }
02354 delete [] threads;
02355 delete [] parms;
02356
02357 return MEASURE_NOERR;
02358 }
02359
02360
02361 int measure_geom(MoleculeList *mlist, int *molid, int *atmid, ResizeArray<float> *gValues,
02362 int frame, int first, int last, int defmolid, int geomtype) {
02363 int i, ret_val;
02364 for(i=0; i < geomtype; i++) {
02365
02366 if(i > 0 && molid[i-1]==molid[i] && atmid[i-1]==atmid[i]) {
02367 printf("measure_geom: %i/%i %i/%i\n", molid[i-1],atmid[i-1],molid[i],atmid[i]);
02368 return MEASURE_ERR_REPEATEDATOM;
02369 }
02370 }
02371
02372 float value;
02373 int max_ts, orig_ts;
02374
02375
02376 Molecule *mol = mlist->mol_from_id(defmolid);
02377 if( !mol )
02378 return MEASURE_ERR_NOMOLECULE;
02379
02380
02381 if((orig_ts = mol->frame()) < 0)
02382 return MEASURE_ERR_NOFRAMES;
02383
02384
02385 max_ts = mol->numframes()-1;
02386 if (frame<0) {
02387 if (first<0 && last<0) first = last = orig_ts;
02388 if (last<0 || last>max_ts) last = max_ts;
02389 if (first<0) first = 0;
02390 } else {
02391 if (frame>max_ts) frame = max_ts;
02392 first = last = frame;
02393 }
02394
02395
02396 for(i=first; i <= last; i++) {
02397 mol->override_current_frame(i);
02398 switch (geomtype) {
02399 case MEASURE_BOND:
02400 if ((ret_val=calculate_bond(mlist, molid, atmid, &value))<0)
02401 return ret_val;
02402 gValues->append(value);
02403 break;
02404 case MEASURE_ANGLE:
02405 if ((ret_val=calculate_angle(mlist, molid, atmid, &value))<0)
02406 return ret_val;
02407 gValues->append(value);
02408 break;
02409 case MEASURE_DIHED:
02410 if ((ret_val=calculate_dihed(mlist, molid, atmid, &value))<0)
02411 return ret_val;
02412 gValues->append(value);
02413 break;
02414 }
02415 }
02416
02417
02418 mol->override_current_frame(orig_ts);
02419
02420 return MEASURE_NOERR;
02421 }
02422
02423
02424
02425 int calculate_bond(MoleculeList *mlist, int *molid, int *atmid, float *value) {
02426
02427
02428 int ret_val;
02429 float pos1[3], pos2[3];
02430 if ((ret_val=normal_atom_coord(mlist->mol_from_id(molid[0]), atmid[0], pos1))<0)
02431 return ret_val;
02432 if ((ret_val=normal_atom_coord(mlist->mol_from_id(molid[1]), atmid[1], pos2))<0)
02433 return ret_val;
02434
02435 vec_sub(pos2, pos2, pos1);
02436 *value = norm(pos2);
02437
02438 return MEASURE_NOERR;
02439 }
02440
02441
02442 int calculate_angle(MoleculeList *mlist, int *molid, int *atmid, float *value) {
02443
02444
02445 int ret_val;
02446 float pos1[3], pos2[3], pos3[3], r1[3], r2[3];
02447 if((ret_val=normal_atom_coord(mlist->mol_from_id(molid[0]), atmid[0], pos1))<0)
02448 return ret_val;
02449 if((ret_val=normal_atom_coord(mlist->mol_from_id(molid[1]), atmid[1], pos2))<0)
02450 return ret_val;
02451 if((ret_val=normal_atom_coord(mlist->mol_from_id(molid[2]), atmid[2], pos3))<0)
02452 return ret_val;
02453
02454 vec_sub(r1, pos1, pos2);
02455 vec_sub(r2, pos3, pos2);
02456 *value = angle(r1, r2);
02457
02458 return MEASURE_NOERR;
02459 }
02460
02461
02462 int calculate_dihed(MoleculeList *mlist, int *molid, int *atmid, float *value) {
02463
02464
02465 int ret_val;
02466 float pos1[3], pos2[3], pos3[3], pos4[3];
02467 if((ret_val=normal_atom_coord(mlist->mol_from_id(molid[0]), atmid[0], pos1))<0)
02468 return ret_val;
02469 if((ret_val=normal_atom_coord(mlist->mol_from_id(molid[1]), atmid[1], pos2))<0)
02470 return ret_val;
02471 if((ret_val=normal_atom_coord(mlist->mol_from_id(molid[2]), atmid[2], pos3))<0)
02472 return ret_val;
02473 if((ret_val=normal_atom_coord(mlist->mol_from_id(molid[3]), atmid[3], pos4))<0)
02474 return ret_val;
02475
02476 *value = dihedral(pos1, pos2, pos3, pos4);
02477
02478 return MEASURE_NOERR;
02479 }
02480
02481
02482
02483
02484 int normal_atom_coord(Molecule *mol, int a, float *pos) {
02485 Timestep *now;
02486
02487 int cell[3];
02488 memset(cell, 0, 3L*sizeof(int));
02489
02490
02491 int ret_val = check_mol(mol, a);
02492 if (ret_val<0)
02493 return ret_val;
02494
02495 if ((now = mol->current())) {
02496 memcpy((void *)pos, (void *)(now->pos + 3L*a), 3L*sizeof(float));
02497
02498
02499 Matrix4 mat;
02500 now->get_transform_from_cell(cell, mat);
02501 mat.multpoint3d(pos, pos);
02502
02503 return MEASURE_NOERR;
02504 }
02505
02506
02507 return MEASURE_ERR_NOFRAMES;
02508 }
02509
02510
02511
02512
02513 int check_mol(Molecule *mol, int a) {
02514
02515 if (!mol)
02516 return MEASURE_ERR_NOMOLECULE;
02517 if (a < 0 || a >= mol->nAtoms)
02518 return MEASURE_ERR_BADATOMID;
02519
02520 return MEASURE_NOERR;
02521 }
02522
02523
02524 int measure_energy(MoleculeList *mlist, int *molid, int *atmid, int natoms, ResizeArray<float> *gValues,
02525 int frame, int first, int last, int defmolid, double *params, int geomtype) {
02526 int i, ret_val;
02527 for(i=0; i < natoms; i++) {
02528
02529 if(i > 0 && molid[i-1]==molid[i] && atmid[i-1]==atmid[i]) {
02530 printf("measure_energy: %i/%i %i/%i\n", molid[i-1],atmid[i-1],molid[i],atmid[i]);
02531 return MEASURE_ERR_REPEATEDATOM;
02532 }
02533 }
02534
02535 float value;
02536 int max_ts, orig_ts;
02537
02538
02539 Molecule *mol = mlist->mol_from_id(defmolid);
02540 if( !mol )
02541 return MEASURE_ERR_NOMOLECULE;
02542
02543
02544 if((orig_ts = mol->frame()) < 0)
02545 return MEASURE_ERR_NOFRAMES;
02546
02547
02548 max_ts = mol->numframes()-1;
02549 if (frame==-1) {
02550 if (first<0 && last<0) first = last = orig_ts;
02551 if (last<0 || last>max_ts) last = max_ts;
02552 if (first<0) first = 0;
02553 } else {
02554 if (frame>max_ts || frame==-2) frame = max_ts;
02555 first = last = frame;
02556 }
02557
02558
02559 for(i=first; i <= last; i++) {
02560 mol->override_current_frame(i);
02561 switch (geomtype) {
02562 case MEASURE_BOND:
02563 if ((ret_val=compute_bond_energy(mlist, molid, atmid, &value, (float) params[0], (float) params[1]))<0)
02564 return ret_val;
02565 gValues->append(value);
02566 break;
02567 case MEASURE_ANGLE:
02568 if ((ret_val=compute_angle_energy(mlist, molid, atmid, &value, (float) params[0], (float) params[1], (float) params[2], (float) params[3]))<0)
02569 return ret_val;
02570 gValues->append(value);
02571 break;
02572 case MEASURE_DIHED:
02573 if ((ret_val=compute_dihed_energy(mlist, molid, atmid, &value, (float) params[0], int(params[1]), (float) params[2]))<0)
02574 return ret_val;
02575 gValues->append(value);
02576 break;
02577 case MEASURE_IMPRP:
02578 if ((ret_val=compute_imprp_energy(mlist, molid, atmid, &value, (float) params[0], (float) params[1]))<0)
02579 return ret_val;
02580 gValues->append(value);
02581 break;
02582 case MEASURE_VDW:
02583 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)
02584 return ret_val;
02585 gValues->append(value);
02586 break;
02587 case MEASURE_ELECT:
02588 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)
02589 return ret_val;
02590 gValues->append(value);
02591 break;
02592 }
02593 }
02594
02595
02596 mol->override_current_frame(orig_ts);
02597
02598 return MEASURE_NOERR;
02599 }
02600
02601
02602 int compute_bond_energy(MoleculeList *mlist, int *molid, int *atmid, float *energy, float k, float x0) {
02603 int ret_val;
02604 float dist;
02605
02606
02607 if ((ret_val=calculate_bond(mlist, molid, atmid, &dist))<0)
02608 return ret_val;
02609 float x = dist-x0;
02610 *energy = k*x*x;
02611
02612 return MEASURE_NOERR;
02613 }
02614
02615
02616 int compute_angle_energy(MoleculeList *mlist, int *molid, int *atmid, float *energy,
02617 float k, float x0, float kub, float s0) {
02618 int ret_val;
02619 float value;
02620
02621
02622 if ((ret_val=calculate_angle(mlist, molid, atmid, &value))<0)
02623 return ret_val;
02624 float x = (float) DEGTORAD((value-x0));
02625 float s = 0.0f;
02626
02627 if (kub>0.0f) {
02628 int twoatoms[2];
02629 twoatoms[0] = atmid[0];
02630 twoatoms[1] = atmid[2];
02631 if ((ret_val=calculate_bond(mlist, molid, twoatoms, &value))<0)
02632 return ret_val;
02633 s = value-s0;
02634 }
02635
02636 *energy = k*x*x + kub*s*s;
02637
02638 return MEASURE_NOERR;
02639 }
02640
02641
02642 int compute_dihed_energy(MoleculeList *mlist, int *molid, int *atmid, float *energy,
02643 float k, int n, float delta) {
02644 int ret_val;
02645 float value;
02646
02647
02648 if ((ret_val=calculate_dihed(mlist, molid, atmid, &value))<0)
02649 return ret_val;
02650 *energy = k*(1+cosf((float) (DEGTORAD((n*value-delta)))));
02651
02652 return MEASURE_NOERR;
02653 }
02654
02655
02656 int compute_imprp_energy(MoleculeList *mlist, int *molid, int *atmid, float *energy,
02657 float k, float x0) {
02658 int ret_val;
02659 float value;
02660
02661
02662 if ((ret_val=calculate_dihed(mlist, molid, atmid, &value))<0)
02663 return ret_val;
02664 float x = (float) (DEGTORAD((value-x0)));
02665 *energy = k*x*x;
02666
02667 return MEASURE_NOERR;
02668 }
02669
02670
02671
02672
02673
02674 int compute_vdw_energy(MoleculeList *mlist, int *molid, int *atmid, float *energy, float eps1, float rmin1,
02675 float eps2, float rmin2, float cutoff, float switchdist) {
02676 int ret_val;
02677 float dist;
02678
02679
02680 if ((ret_val=calculate_bond(mlist, molid, atmid, &dist))<0)
02681 return ret_val;
02682
02683 float sw=1.0;
02684 if (switchdist>0.0 && cutoff>0.0) {
02685 if (dist>=cutoff) {
02686 sw = 0.0;
02687 } else if (dist>=switchdist) {
02688
02689 float dist2 = dist*dist;
02690 float cut2 = cutoff*cutoff;
02691 float switch2 = switchdist*switchdist;
02692 float s = cut2-dist2;
02693 float range = cut2-switch2;
02694 sw = s*s*(cut2+2*dist2-3*switch2)/(range*range*range);
02695 }
02696 }
02697
02698 float term6 = (float) powf((rmin1+rmin2)/dist,6);
02699 *energy = sqrtf(eps1*eps2)*(term6*term6 - 2.0f*term6)*sw;
02700
02701 return MEASURE_NOERR;
02702 }
02703
02704 int compute_elect_energy(MoleculeList *mlist, int *molid, int *atmid, float *energy, float q1, float q2,
02705 bool flag1, bool flag2, float cutoff) {
02706 int ret_val;
02707 float dist;
02708
02709
02710 if ((ret_val=calculate_bond(mlist, molid, atmid, &dist))<0)
02711 return ret_val;
02712
02713
02714 if (!flag1) q1 = mlist->mol_from_id(molid[0])->charge()[atmid[0]];
02715 if (!flag2) q2 = mlist->mol_from_id(molid[0])->charge()[atmid[1]];
02716
02717 if (cutoff>0.0) {
02718 if (dist<cutoff) {
02719 float efac = 1.0f-dist*dist/(cutoff*cutoff);
02720 *energy = 332.0636f*q1*q2/dist*efac*efac;
02721 } else {
02722 *energy = 0.0f;
02723 }
02724 } else {
02725 *energy = 332.0636f*q1*q2/dist;
02726 }
02727
02728 return MEASURE_NOERR;
02729 }
02730
02731
02732
02733
02734 static void center_of_mass(AtomSel *sel, MoleculeList *mlist, float *rcom) {
02735 int i;
02736 float m = 0, mtot = 0;
02737 Molecule *mol = mlist->mol_from_id(sel->molid());
02738
02739
02740 const float *mass = mol->mass();
02741
02742
02743 const float *pos = sel->coordinates(mlist);
02744
02745 memset(rcom, 0, 3L*sizeof(float));
02746
02747
02748 for (i=sel->firstsel; i<=sel->lastsel; i++) {
02749 if (sel->on[i]) {
02750 long ind = i * 3L;
02751
02752 m = mass[i];
02753
02754 rcom[0] += m*pos[ind ];
02755 rcom[1] += m*pos[ind + 1];
02756 rcom[2] += m*pos[ind + 2];
02757
02758
02759 mtot += m;
02760 }
02761 }
02762
02763 rcom[0] /= mtot;
02764 rcom[1] /= mtot;
02765 rcom[2] /= mtot;
02766 }
02767
02768
02769
02770
02771
02772
02773
02774
02775 extern int measure_inertia(AtomSel *sel, MoleculeList *mlist, const float *coor, float rcom[3],
02776 float priaxes[3][3], float itensor[4][4], float evalue[3]) {
02777 if (!sel) return MEASURE_ERR_NOSEL;
02778 if (sel->num_atoms == 0) return MEASURE_ERR_NOATOMS;
02779
02780 Molecule *mol = mlist->mol_from_id(sel->molid());
02781
02782 float x, y, z, m;
02783 float Ixx=0, Iyy=0, Izz=0, Ixy=0, Ixz=0, Iyz=0;
02784 int i,j=0;
02785
02786
02787
02788 memset(itensor, 0, 16L*sizeof(float));
02789 itensor[3][3] = 1.0;
02790
02791
02792 center_of_mass(sel, mlist, rcom);
02793
02794
02795 const float *pos = sel->coordinates(mlist);
02796
02797
02798 const float *mass = mol->mass();
02799
02800
02801
02802 for (i=sel->firstsel; i<=sel->lastsel; i++) {
02803 if (sel->on[i]) {
02804
02805 if (coor) {
02806
02807 x = coor[j*3L ] - rcom[0];
02808 y = coor[j*3L + 1] - rcom[1];
02809 z = coor[j*3L + 2] - rcom[2];
02810 j++;
02811 } else {
02812
02813 x = pos[i*3L ] - rcom[0];
02814 y = pos[i*3L + 1] - rcom[1];
02815 z = pos[i*3L + 2] - rcom[2];
02816 }
02817
02818 m = mass[i];
02819
02820 Ixx += m*(y*y+z*z);
02821 Iyy += m*(x*x+z*z);
02822 Izz += m*(x*x+y*y);
02823 Ixy -= m*x*y;
02824 Ixz -= m*x*z;
02825 Iyz -= m*y*z;
02826 }
02827 }
02828
02829 itensor[0][0] = Ixx;
02830 itensor[1][1] = Iyy;
02831 itensor[2][2] = Izz;
02832 itensor[0][1] = Ixy;
02833 itensor[1][0] = Ixy;
02834 itensor[0][2] = Ixz;
02835 itensor[2][0] = Ixz;
02836 itensor[1][2] = Iyz;
02837 itensor[2][1] = Iyz;
02838
02839
02840
02841 float evector[3][3];
02842 if (jacobi(itensor,evalue,evector) != 0) return MEASURE_ERR_NONZEROJACOBI;
02843
02844
02845 float vectmp;
02846 vectmp=evector[0][1]; evector[0][1]=evector[1][0]; evector[1][0]=vectmp;
02847 vectmp=evector[0][2]; evector[0][2]=evector[2][0]; evector[2][0]=vectmp;
02848 vectmp=evector[2][1]; evector[2][1]=evector[1][2]; evector[1][2]=vectmp;
02849
02850
02851
02852
02853 float *a[3];
02854 a[0] = evector[0];
02855 a[1] = evector[1];
02856 a[2] = evector[2];
02857
02858
02859 #define SWAP(qq,ww) { \
02860 float v; float *v1; \
02861 v = evalue[qq]; evalue[qq] = evalue[ww]; evalue[ww] = v; \
02862 v1 = a[qq]; a[qq] = a[ww]; a[ww] = v1; \
02863 }
02864 if (evalue[0] < evalue[1]) {
02865 SWAP(0, 1);
02866 }
02867 if (evalue[0] < evalue[2]) {
02868 SWAP(0, 2);
02869 }
02870 if (evalue[1] < evalue[2]) {
02871 SWAP(1, 2);
02872 }
02873
02874 #if 0
02875
02876
02877 if (evalue[1]/evalue[0]>0.1 && fabs(evalue[1]-evalue[2])/evalue[0]<0.05
02878 && fabs(evalue[0]-evalue[1])/evalue[0]>0.05) {
02879 msgInfo << "Principal axes of inertia 2 and 3 are not unique!" << sendmsg;
02880 }
02881 #endif
02882
02883 for (i=0; i<3; i++) {
02884 for (j=0; j<3; j++)
02885 priaxes[i][j] = a[i][j];
02886 }
02887
02888 return MEASURE_NOERR;
02889 }
02890