SpatialSearch.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: SpatialSearch.C,v $
00013 * $Author: johns $ $Locker: $ $State: Exp $
00014 * $Revision: 1.26 $ $Date: 2020年07月08日 04:23:47 $
00015 *
00016 ***************************************************************************
00017 * DESCRIPTION:
00018 *
00019 * Distance based bond search code 
00020 *
00021 ***************************************************************************/
00022 
00023 #include <stdio.h>
00024 #include <math.h>
00025 #include <stdlib.h>
00026 #include "BondSearch.h"
00027 #include "Timestep.h"
00028 #include "BaseMolecule.h"
00029 #include "Molecule.h"
00030 #include "Inform.h"
00031 #include "WKFThreads.h"
00032 #include "WKFUtils.h"
00033 #include <ctype.h> // needed for isdigit()
00034 #include <string.h>
00035 
00036 static void add_link(GridSearchPair *link, int i, int j) {
00037 link->next = (GridSearchPair *) malloc(sizeof(GridSearchPair));
00038 link->next->ind1 = i;
00039 link->next->ind2 = j;
00040 link->next->next = NULL;
00041 }
00042 
00043 void find_minmax_all(const float *pos, int n, float *min, float *max) {
00044 float x1, x2, y1, y2, z1, z2;
00045 int i=0;
00046 
00047 // return immediately if there are no atoms, or no atoms are on.
00048 if (n < 1) return;
00049 
00050 // initialize min/max to first 'on' atom, and advance the counter.
00051 pos += 3L*i;
00052 x1 = x2 = pos[0];
00053 y1 = y2 = pos[1];
00054 z1 = z2 = pos[2];
00055 pos += 3;
00056 i++;
00057 
00058 for (; i < n; i++) {
00059 if (pos[0] < x1) x1 = pos[0];
00060 if (pos[0] > x2) x2 = pos[0];
00061 if (pos[1] < y1) y1 = pos[1];
00062 if (pos[1] > y2) y2 = pos[1];
00063 if (pos[2] < z1) z1 = pos[2];
00064 if (pos[2] > z2) z2 = pos[2];
00065 pos += 3;
00066 }
00067 min[0] = x1; min[1] = y1; min[2] = z1;
00068 max[0] = x2; max[1] = y2; max[2] = z2;
00069 }
00070 
00071 
00072 // find minmax of selected atoms. Return true if minmax was found;
00073 // false if all flags are zero.
00074 int find_minmax_selected(int n, const int *flgs, const float *pos,
00075 float &_xmin, float &_ymin, float &_zmin,
00076 float &_xmax, float &_ymax, float &_zmax) {
00077 int i;
00078 float xmin, xmax, ymin, ymax, zmin, zmax;
00079 for (i=0; i<n; i++) if (flgs[i]) break;
00080 if (i==n) return FALSE;
00081 pos += 3L*i;
00082 xmin=xmax=pos[0];
00083 ymin=ymax=pos[1];
00084 zmin=zmax=pos[2];
00085 pos += 3;
00086 i += 1;
00087 for (; i<n; i++, pos += 3) {
00088 if (!flgs[i]) continue;
00089 const float xi=pos[0];
00090 const float yi=pos[1];
00091 const float zi=pos[2];
00092 if (xmin>xi) xmin=xi;
00093 if (ymin>yi) ymin=yi;
00094 if (zmin>zi) zmin=zi;
00095 if (xmax<xi) xmax=xi;
00096 if (ymax<yi) ymax=yi;
00097 if (zmax<zi) zmax=zi;
00098 }
00099 
00100 // check for cases where NaN coordinates propagated to min/max bounds
00101 if (!(xmax >= xmin && ymax >= ymin && zmax >= zmin)) {
00102 msgErr << "find_minmax_selected: NaN coordinates in bounds!" << sendmsg;
00103 return FALSE;
00104 }
00105 
00106 _xmin=xmin;
00107 _ymin=ymin;
00108 _zmin=zmin;
00109 _xmax=xmax;
00110 _ymax=ymax;
00111 _zmax=zmax;
00112 return TRUE;
00113 }
00114 
00115 
00116 // old routine for finding min max of selected atoms
00117 void find_minmax(const float *pos, int n, const int *on, 
00118 float *min, float *max, int *oncount) {
00119 float x1, x2, y1, y2, z1, z2;
00120 int i, numon;
00121 
00122 // return immediately if there are no atoms, or no atoms are on.
00123 if (n < 1) return;
00124 
00125 // init on count
00126 numon = 0;
00127 
00128 // find first on atom
00129 for (i=0; i<n; i++) {
00130 if (on[i]) {
00131 numon++;
00132 break;
00133 }
00134 }
00135 if (i==n) {
00136 if (oncount != NULL) 
00137 *oncount = numon;
00138 return;
00139 }
00140 
00141 // initialize min/max to first 'on' atom, and advance the counter.
00142 pos += 3L*i;
00143 x1 = x2 = pos[0];
00144 y1 = y2 = pos[1];
00145 z1 = z2 = pos[2];
00146 pos += 3;
00147 i++;
00148 
00149 for (; i < n; i++) {
00150 if (on[i]) {
00151 if (pos[0] < x1) x1 = pos[0];
00152 else if (pos[0] > x2) x2 = pos[0];
00153 if (pos[1] < y1) y1 = pos[1];
00154 else if (pos[1] > y2) y2 = pos[1];
00155 if (pos[2] < z1) z1 = pos[2];
00156 else if (pos[2] > z2) z2 = pos[2];
00157 numon++;
00158 }
00159 pos += 3;
00160 }
00161 min[0] = x1; min[1] = y1; min[2] = z1;
00162 max[0] = x2; max[1] = y2; max[2] = z2;
00163 
00164 if (oncount != NULL) 
00165 *oncount = numon;
00166 }
00167 
00168 int make_neighborlist(int **nbrlist, int xb, int yb, int zb) {
00169 int xi, yi, zi, aindex, xytotb;
00170 
00171 if (nbrlist == NULL)
00172 return -1;
00173 
00174 xytotb = xb * yb;
00175 aindex = 0;
00176 for (zi=0; zi<zb; zi++) {
00177 for (yi=0; yi<yb; yi++) {
00178 for (xi=0; xi<xb; xi++) {
00179 int nbrs[15]; // 14 neighbors, and a -1 to mark end of the list
00180 int n=0; 
00181 nbrs[n++] = aindex; 
00182 if (xi < xb-1) nbrs[n++] = aindex + 1;
00183 if (yi < yb-1) nbrs[n++] = aindex + xb;
00184 if (zi < zb-1) nbrs[n++] = aindex + xytotb;
00185 if (xi < (xb-1) && yi < (yb-1)) nbrs[n++] = aindex + xb + 1;
00186 if (xi < (xb-1) && zi < (zb-1)) nbrs[n++] = aindex + xytotb + 1;
00187 if (yi < (yb-1) && zi < (zb-1)) nbrs[n++] = aindex + xytotb + xb;
00188 if (xi < (xb-1) && yi > 0) nbrs[n++] = aindex - xb + 1;
00189 if (xi > 0 && zi < (zb-1)) nbrs[n++] = aindex + xytotb - 1;
00190 if (yi > 0 && zi < (zb-1)) nbrs[n++] = aindex + xytotb - xb;
00191 if (xi < (xb-1) && yi < (yb-1) && zi < (zb-1))
00192 nbrs[n++] = aindex + xytotb + xb + 1;
00193 if (xi > 0 && yi < (yb-1) && zi < (zb-1))
00194 nbrs[n++] = aindex + xytotb + xb - 1;
00195 if (xi < (xb-1) && yi > 0 && zi < (zb-1))
00196 nbrs[n++] = aindex + xytotb - xb + 1;
00197 if (xi > 0 && yi > 0 && zi < (zb-1))
00198 nbrs[n++] = aindex + xytotb - xb - 1;
00199 nbrs[n++] = -1; // mark end of list
00200 
00201 int *lst = (int *) malloc(n*sizeof(int));
00202 if (lst == NULL)
00203 return -1; // return on failed allocations
00204 memcpy(lst, nbrs, n*sizeof(int));
00205 nbrlist[aindex] = lst;
00206 aindex++;
00207 }
00208 }
00209 }
00210 
00211 return 0;
00212 }
00213 
00214 
00215 // symetrical version of neighbourlist (each cell has 27 neighbours, incl. itself, instead of just 14)
00216 static int make_neighborlist_sym(int **nbrlist, int xb, int yb, int zb) {
00217 int xi, yi, zi, aindex, xytotb;
00218 
00219 if (nbrlist == NULL)
00220 return -1;
00221 
00222 xytotb = xb * yb;
00223 aindex = 0;
00224 for (zi=0; zi<zb; zi++) {
00225 for (yi=0; yi<yb; yi++) {
00226 for (xi=0; xi<xb; xi++) {
00227 int nbrs[28]; // 27 neighbors, and a -1 to mark end of the list
00228 int n=0; 
00229 nbrs[n++] = aindex; 
00230 if (xi < xb-1) nbrs[n++] = aindex + 1;
00231 if (yi < yb-1) nbrs[n++] = aindex + xb;
00232 if (zi < zb-1) nbrs[n++] = aindex + xytotb;
00233 if (xi < (xb-1) && yi < (yb-1)) nbrs[n++] = aindex + xb + 1;
00234 if (xi < (xb-1) && zi < (zb-1)) nbrs[n++] = aindex + xytotb + 1;
00235 if (yi < (yb-1) && zi < (zb-1)) nbrs[n++] = aindex + xytotb + xb;
00236 if (xi < (xb-1) && yi) nbrs[n++] = aindex - xb + 1;
00237 if (xi && zi < (zb-1)) nbrs[n++] = aindex + xytotb - 1;
00238 if (yi && zi < (zb-1)) nbrs[n++] = aindex + xytotb - xb;
00239 if (xi < (xb-1) && yi < (yb-1) && zi < (zb-1)) nbrs[n++] = aindex + xytotb + xb + 1;
00240 if (xi && yi < (yb-1) && zi < (zb-1)) nbrs[n++] = aindex + xytotb + xb - 1;
00241 if (xi < (xb-1) && yi && zi < (zb-1)) nbrs[n++] = aindex + xytotb - xb + 1;
00242 if (xi && yi && zi < (zb-1)) nbrs[n++] = aindex + xytotb - xb - 1;
00243 
00244 if (xi) nbrs[n++] = aindex - 1;
00245 if (yi) nbrs[n++] = aindex - xb;
00246 if (zi) nbrs[n++] = aindex - xytotb;
00247 if (xi && yi) nbrs[n++] = aindex - xb - 1;
00248 if (xi && zi) nbrs[n++] = aindex - xytotb - 1;
00249 if (yi && zi) nbrs[n++] = aindex - xytotb - xb;
00250 if (xi && yi < (yb-1)) nbrs[n++] = aindex + xb - 1; 
00251 if (xi < (xb-1) && zi) nbrs[n++] = aindex - xytotb + 1;
00252 if (yi < (yb-1) && zi) nbrs[n++] = aindex - xytotb + xb;
00253 if (xi && yi && zi) nbrs[n++] = aindex - xytotb - xb - 1; 
00254 if (xi < (xb-1) && yi && zi) nbrs[n++] = aindex - xytotb - xb + 1;
00255 if (xi && yi < (yb-1) && zi) nbrs[n++] = aindex - xytotb + xb - 1;
00256 if (xi < (xb-1) && yi < (yb-1) && zi) nbrs[n++] = aindex - xytotb + xb + 1;
00257 nbrs[n++] = -1; // mark end of list
00258 
00259 int *lst = (int *) malloc(n*sizeof(int));
00260 if (lst == NULL)
00261 return -1; // return on failed allocations
00262 memcpy(lst, nbrs, n*sizeof(int));
00263 nbrlist[aindex] = lst;
00264 aindex++;
00265 }
00266 }
00267 }
00268 
00269 return 0;
00270 }
00271 
00272 
00273 GridSearchPair *vmd_gridsearch1(const float *pos,int natoms, const int *on, 
00274 float pairdist, int allow_double_counting, int maxpairs) {
00275 float min[3]={0,0,0}, max[3]={0,0,0};
00276 float sqdist;
00277 int i, j, xb, yb, zb, xytotb, totb, aindex;
00278 int **boxatom, *numinbox, *maxinbox, **nbrlist;
00279 int numon = 0;
00280 float sidelen[3], volume;
00281 int paircount = 0;
00282 int maxpairsreached = 0;
00283 sqdist = pairdist * pairdist;
00284 
00285 // find bounding box for selected atoms, and number of atoms in selection.
00286 find_minmax(pos, natoms, on, min, max, &numon);
00287 
00288 // do sanity checks and complain if we've got bogus atom coordinates,
00289 // we shouldn't ever have density higher than 0.1 atom/A^3, but we'll
00290 // be generous and allow much higher densities. 
00291 vec_sub(sidelen, max, min);
00292 // include estimate for atom radius (1 Angstrom) in volume determination
00293 volume = fabsf((sidelen[0] + 2.0f) * (sidelen[1] + 2.0f) * (sidelen[2] + 2.0f));
00294 if (isnan(volume)) {
00295 msgErr << "vmd_gridsearch1: bounding volume yielded NaN" << sendmsg;
00296 return NULL;
00297 }
00298 
00299 if (maxpairs != -1) {
00300 if ((numon / volume) > 1.0) {
00301 msgWarn << "vmd_gridsearch1: insane atom density" << sendmsg;
00302 }
00303 }
00304 
00305 // I don't want the grid to get too large, otherwise I could run out
00306 // of memory. Octrees would be cool, but I'll just limit the grid size
00307 // and let the performance degrade a little for pathological systems.
00308 // Note that sqdist is what gets used for the actual distance checks;
00309 // from here on out pairdist is only used to set the grid size, so we 
00310 // can set it to anything larger than the original pairdist.
00311 const int MAXBOXES = 4000000;
00312 totb = MAXBOXES + 1;
00313 
00314 float newpairdist = pairdist;
00315 float xrange = max[0]-min[0];
00316 float yrange = max[1]-min[1];
00317 float zrange = max[2]-min[2];
00318 do {
00319 pairdist = newpairdist;
00320 const float invpairdist = 1.0f / pairdist;
00321 xb = ((int)(xrange*invpairdist))+1;
00322 yb = ((int)(yrange*invpairdist))+1;
00323 zb = ((int)(zrange*invpairdist))+1;
00324 xytotb = yb * xb;
00325 totb = xytotb * zb;
00326 newpairdist = pairdist * 1.26f; // cbrt(2) is about 1.26
00327 } while (totb > MAXBOXES || totb < 1); // check for integer wraparound too
00328 
00329 // 2. Sort each atom into appropriate bins
00330 boxatom = (int **) calloc(1, totb*sizeof(int *));
00331 numinbox = (int *) calloc(1, totb*sizeof(int));
00332 maxinbox = (int *) calloc(1, totb*sizeof(int));
00333 if (boxatom == NULL || numinbox == NULL || maxinbox == NULL) {
00334 if (boxatom != NULL)
00335 free(boxatom);
00336 if (numinbox != NULL)
00337 free(numinbox);
00338 if (maxinbox != NULL)
00339 free(maxinbox);
00340 msgErr << "Gridsearch memory allocation failed, bailing out" << sendmsg;
00341 return NULL; // ran out of memory, bail out!
00342 }
00343 
00344 const float invpairdist = 1.0f / pairdist;
00345 for (i=0; i<natoms; i++) {
00346 if (on[i]) {
00347 int axb, ayb, azb, aindex, num;
00348 
00349 // compute box index for new atom
00350 const float *loc = pos + 3L*i;
00351 axb = (int)((loc[0] - min[0])*invpairdist);
00352 ayb = (int)((loc[1] - min[1])*invpairdist);
00353 azb = (int)((loc[2] - min[2])*invpairdist);
00354 
00355 // clamp box indices to valid range in case of FP error
00356 if (axb >= xb) axb = xb-1;
00357 if (ayb >= yb) ayb = yb-1;
00358 if (azb >= zb) azb = zb-1;
00359 
00360 aindex = azb * xytotb + ayb * xb + axb;
00361 
00362 // grow box if necessary
00363 if ((num = numinbox[aindex]) == maxinbox[aindex]) {
00364 boxatom[aindex] = (int *) realloc(boxatom[aindex], (num+4)*sizeof(int));
00365 maxinbox[aindex] += 4;
00366 }
00367 
00368 // store atom index in box
00369 boxatom[aindex][num] = i;
00370 numinbox[aindex]++;
00371 }
00372 }
00373 free(maxinbox);
00374 
00375 nbrlist = (int **) calloc(1, totb*sizeof(int *));
00376 if (make_neighborlist(nbrlist, xb, yb, zb)) {
00377 if (boxatom != NULL) {
00378 for (i=0; i<totb; i++) {
00379 if (boxatom[i] != NULL) free(boxatom[i]);
00380 }
00381 free(boxatom);
00382 }
00383 if (nbrlist != NULL) {
00384 for (i=0; i<totb; i++) {
00385 if (nbrlist[i] != NULL) free(nbrlist[i]);
00386 }
00387 free(nbrlist);
00388 }
00389 free(numinbox);
00390 msgErr << "Gridsearch memory allocation failed, bailing out" << sendmsg;
00391 return NULL; // ran out of memory, bail out!
00392 }
00393 
00394 // setup head of pairlist
00395 GridSearchPair *head, *cur;
00396 head = (GridSearchPair *) malloc(sizeof(GridSearchPair));
00397 head->next = NULL;
00398 cur = head;
00399 
00400 wkfmsgtimer *msgt = wkf_msg_timer_create(5);
00401 for (aindex = 0; (aindex < totb) && (!maxpairsreached); aindex++) {
00402 int *tmpbox, *tmpnbr, *nbr;
00403 tmpbox = boxatom[aindex];
00404 tmpnbr = nbrlist[aindex];
00405 
00406 if (wkf_msg_timer_timeout(msgt)) {
00407 char tmpbuf[128] = { 0 };
00408 sprintf(tmpbuf, "%6.2f", (100.0f * aindex) / (float) totb);
00409 msgInfo << "vmd_gridsearch1: " << tmpbuf << "% complete" << sendmsg;
00410 }
00411 
00412 for (nbr = tmpnbr; (*nbr != -1) && (!maxpairsreached); nbr++) {
00413 int *nbrbox = boxatom[*nbr];
00414 for (i=0; (i<numinbox[aindex]) && (!maxpairsreached); i++) {
00415 int ind1 = tmpbox[i];
00416 if (!on[ind1]) 
00417 continue;
00418 const float *p1 = pos + 3L*ind1;
00419 int startj = 0;
00420 if (aindex == *nbr) startj = i+1;
00421 for (j=startj; (j<numinbox[*nbr]) && (!maxpairsreached); j++) {
00422 int ind2 = nbrbox[j];
00423 if (on[ind2]) {
00424 const float *p2 = pos + 3L*ind2;
00425 float ds2 = distance2(p1, p2);
00426 
00427 // ignore pairs between atoms with nearly identical coords
00428 if (ds2 < 0.001)
00429 continue;
00430 
00431 if (ds2 > sqdist) 
00432 continue;
00433 
00434 if (maxpairs > 0) {
00435 if (paircount >= maxpairs) {
00436 maxpairsreached = 1;
00437 continue; 
00438 } 
00439 }
00440 
00441 add_link(cur, ind1, ind2);
00442 paircount++;
00443 
00444 // XXX double-counting still ignores atoms with same coords...
00445 if (allow_double_counting) { 
00446 add_link(cur, ind2, ind1); 
00447 paircount++;
00448 }
00449 cur = cur->next;
00450 cur->next = NULL;
00451 }
00452 }
00453 }
00454 }
00455 } 
00456 
00457 for (i=0; i<totb; i++) {
00458 free(boxatom[i]);
00459 free(nbrlist[i]);
00460 }
00461 free(boxatom);
00462 free(nbrlist);
00463 free(numinbox);
00464 
00465 cur = head->next;
00466 free(head);
00467 
00468 if (maxpairsreached) 
00469 msgErr << "vmdgridsearch1: exceeded pairlist sanity check, aborted" << sendmsg;
00470 
00471 wkf_msg_timer_destroy(msgt);
00472 
00473 return cur;
00474 }
00475 
00476 GridSearchPair *vmd_gridsearch2(const float *pos,int natoms, 
00477 const int *A,const int *B, float pairdist, int maxpairs) {
00478 float min[3]={0,0,0}, max[3]={0,0,0};
00479 float sqdist;
00480 int i, j, xb, yb, zb, xytotb, totb, aindex;
00481 int **boxatom, *numinbox, *maxinbox, **nbrlist;
00482 float sidelen[3], volume;
00483 int numon = 0;
00484 int paircount = 0;
00485 int maxpairsreached = 0;
00486 sqdist = pairdist * pairdist;
00487 
00488 // on = union(A,B). We grid all atoms, then allow pairs only when one
00489 // atom is A and the other atom is B. An alternative scheme would be to
00490 // to grid the A and B atoms separately, and/or bin them separately, so
00491 // there would be fewer rejected pairs. 
00492 
00493 int *on = (int *) malloc(natoms*sizeof(int)); 
00494 for (i=0; i<natoms; i++) {
00495 if (A[i] || B[i]) {
00496 numon++;
00497 on[i] = 1;
00498 } else {
00499 on[i] = 0;
00500 }
00501 }
00502 
00503 // find bounding box for selected atoms
00504 find_minmax(pos, natoms, on, min, max, NULL);
00505 
00506 // do sanity checks and complain if we've got bogus atom coordinates,
00507 // we shouldn't ever have density higher than 0.1 atom/A^3, but we'll
00508 // be generous and allow much higher densities.
00509 vec_sub(sidelen, max, min);
00510 // include estimate for atom radius (1 Angstrom) in volume determination
00511 volume = fabsf((sidelen[0] + 2.0f) * (sidelen[1] + 2.0f) * (sidelen[2] + 2.0f));
00512 if (isnan(volume)) {
00513 msgErr << "vmd_gridsearch2: bounding volume yielded NaN" << sendmsg;
00514 return NULL;
00515 }
00516 
00517 if (maxpairs != -1) {
00518 if ((numon / volume) > 1.0) {
00519 msgWarn << "vmd_gridsearch2: insane atom density" << sendmsg;
00520 }
00521 }
00522 
00523 // I don't want the grid to get too large, otherwise I could run out
00524 // of memory. Octrees would be cool, but I'll just limit the grid size
00525 // and let the performance degrade a little for pathological systems.
00526 // Note that sqdist is what gets used for the actual distance checks;
00527 // from here on out pairdist is only used to set the grid size, so we 
00528 // can set it to anything larger than the original pairdist.
00529 const int MAXBOXES = 4000000;
00530 totb = MAXBOXES + 1;
00531 
00532 float newpairdist = pairdist;
00533 float xrange = max[0]-min[0];
00534 float yrange = max[1]-min[1];
00535 float zrange = max[2]-min[2];
00536 do {
00537 pairdist = newpairdist;
00538 const float invpairdist = 1.0f / pairdist;
00539 xb = ((int)(xrange*invpairdist))+1;
00540 yb = ((int)(yrange*invpairdist))+1;
00541 zb = ((int)(zrange*invpairdist))+1;
00542 xytotb = yb * xb;
00543 totb = xytotb * zb;
00544 newpairdist = pairdist * 1.26f; // cbrt(2) is about 1.26
00545 } while (totb > MAXBOXES || totb < 1); // check for integer wraparound too
00546 
00547 // 2. Sort each atom into appropriate bins
00548 boxatom = (int **) calloc(1, totb*sizeof(int *));
00549 numinbox = (int *) calloc(1, totb*sizeof(int));
00550 maxinbox = (int *) calloc(1, totb*sizeof(int));
00551 if (boxatom == NULL || numinbox == NULL || maxinbox == NULL) {
00552 if (boxatom != NULL)
00553 free(boxatom);
00554 if (numinbox != NULL)
00555 free(numinbox);
00556 if (maxinbox != NULL)
00557 free(maxinbox);
00558 msgErr << "Gridsearch memory allocation failed, bailing out" << sendmsg;
00559 return NULL; // ran out of memory, bail out!
00560 }
00561 
00562 const float invpairdist = 1.0f / pairdist;
00563 for (i=0; i<natoms; i++) {
00564 if (on[i]) {
00565 int axb, ayb, azb, aindex, num;
00566 
00567 // compute box index for new atom
00568 const float *loc = pos + 3L*i;
00569 axb = (int)((loc[0] - min[0])*invpairdist);
00570 ayb = (int)((loc[1] - min[1])*invpairdist);
00571 azb = (int)((loc[2] - min[2])*invpairdist);
00572 
00573 // clamp box indices to valid range in case of FP error
00574 if (axb >= xb) axb = xb-1;
00575 if (ayb >= yb) ayb = yb-1;
00576 if (azb >= zb) azb = zb-1;
00577 
00578 aindex = azb * xytotb + ayb * xb + axb;
00579 
00580 // grow box if necessary
00581 if ((num = numinbox[aindex]) == maxinbox[aindex]) {
00582 boxatom[aindex] = (int *) realloc(boxatom[aindex], (num+4)*sizeof(int));
00583 maxinbox[aindex] += 4;
00584 }
00585 
00586 // store atom index in box
00587 boxatom[aindex][num] = i;
00588 numinbox[aindex]++;
00589 }
00590 }
00591 free(maxinbox);
00592 free(on);
00593 
00594 nbrlist = (int **) calloc(1, totb*sizeof(int *));
00595 if (make_neighborlist(nbrlist, xb, yb, zb)) {
00596 if (boxatom != NULL) {
00597 for (i=0; i<totb; i++) {
00598 if (boxatom[i] != NULL) free(boxatom[i]);
00599 }
00600 free(boxatom);
00601 }
00602 if (nbrlist != NULL) {
00603 for (i=0; i<totb; i++) {
00604 if (nbrlist[i] != NULL) free(nbrlist[i]);
00605 }
00606 free(nbrlist);
00607 }
00608 free(numinbox);
00609 msgErr << "Gridsearch memory allocation failed, bailing out" << sendmsg;
00610 return NULL; // ran out of memory, bail out!
00611 }
00612 
00613 // setup head of pairlist
00614 GridSearchPair *head, *cur;
00615 head = (GridSearchPair *) malloc(sizeof(GridSearchPair));
00616 head->next = NULL;
00617 cur = head;
00618 
00619 for (aindex = 0; aindex < totb; aindex++) {
00620 int *tmpbox, *tmpnbr, *nbr;
00621 tmpbox = boxatom[aindex];
00622 tmpnbr = nbrlist[aindex];
00623 for (nbr = tmpnbr; (*nbr != -1) && (!maxpairsreached); nbr++) {
00624 int *nbrbox = boxatom[*nbr];
00625 for (i=0; (i<numinbox[aindex]) && (!maxpairsreached); i++) {
00626 const float *p1;
00627 int ind1, startj;
00628 ind1 = tmpbox[i];
00629 p1 = pos + 3L*ind1;
00630 startj = 0;
00631 if (aindex == *nbr) startj = i+1;
00632 for (j=startj; (j<numinbox[*nbr]) && (!maxpairsreached); j++) {
00633 const float *p2;
00634 int ind2;
00635 ind2 = nbrbox[j];
00636 if ((A[ind1] && B[ind2]) || (A[ind2] && B[ind1])) {
00637 p2 = pos + 3L*ind2;
00638 if (distance2(p1,p2) > sqdist) continue;
00639 
00640 if (maxpairs > 0) {
00641 if (paircount >= maxpairs) {
00642 maxpairsreached = 1;
00643 continue; 
00644 } 
00645 }
00646 
00647 if (A[ind1]) {
00648 add_link(cur, ind1, ind2);
00649 paircount++;
00650 } else {
00651 add_link(cur, ind2, ind1);
00652 paircount++;
00653 }
00654 cur = cur->next;
00655 cur->next = NULL;
00656 }
00657 }
00658 }
00659 }
00660 } 
00661 
00662 for (i=0; i<totb; i++) {
00663 free(boxatom[i]);
00664 free(nbrlist[i]);
00665 }
00666 free(boxatom);
00667 free(nbrlist);
00668 free(numinbox);
00669 
00670 cur = head->next;
00671 free(head);
00672 
00673 if (maxpairsreached) 
00674 msgErr << "vmdgridsearch2: exceeded pairlist sanity check, aborted" << sendmsg;
00675 
00676 return cur;
00677 }
00678 
00679 
00680 
00681 // Like vmd_gridsearch2, but pairs up atoms from different locations and molecules.
00682 // By default, if (posA == posB), all bonds are unique. Otherwise, double-counting is allowed.
00683 // This can be overridden by setting allow_double_counting (true, false, or default=-1).
00684 GridSearchPair *vmd_gridsearch3(const float *posA, int natomsA, const int *A, 
00685 const float *posB, int natomsB, const int *B, 
00686 float pairdist, int allow_double_counting, int maxpairs) {
00687 
00688 if (!natomsA || !natomsB) return NULL;
00689 
00690 if (allow_double_counting == -1) { //default
00691 if (posA == posB && natomsA == natomsB) 
00692 allow_double_counting = FALSE;
00693 else 
00694 allow_double_counting = TRUE;
00695 }
00696 
00697 // if same molecule and *A[] == *B[], it is a lot faster to use gridsearch1
00698 if (posA == posB && natomsA == natomsB) {
00699 bool is_equal = TRUE;
00700 for (int i=0; i<natomsA && is_equal; i++) {
00701 if (A[i] != B[i])
00702 is_equal = FALSE;
00703 }
00704 if (is_equal)
00705 return vmd_gridsearch1(posA, natomsA, A, pairdist, allow_double_counting, maxpairs);
00706 }
00707 
00708 float min[3], max[3], sqdist;
00709 float minB[3], maxB[3]; //tmp storage
00710 int i, j, xb, yb, zb, xytotb, totb, aindex;
00711 int **boxatomA, *numinboxA, *maxinboxA;
00712 int **boxatomB, *numinboxB, *maxinboxB;
00713 int **nbrlist;
00714 float sidelen[3], volume;
00715 int numonA = 0; int numonB = 0;
00716 int paircount = 0;
00717 int maxpairsreached = 0;
00718 sqdist = pairdist * pairdist;
00719 
00720 // 1. Find grid size for binning
00721 // find bounding box for selected atoms
00722 find_minmax(posA, natomsA, A, min, max, &numonA);
00723 find_minmax(posB, natomsB, B, minB, maxB, &numonB);
00724 
00725 // If no atoms were selected we don't have to go on
00726 if (!numonA || !numonB) return NULL;
00727 
00728 for (i=0; i<3; i++) {
00729 if (minB[i] < min[i]) min[i] = minB[i];
00730 if (maxB[i] > max[i]) max[i] = maxB[i];
00731 }
00732 
00733 // do sanity checks and complain if we've got bogus atom coordinates,
00734 // we shouldn't ever have density higher than 0.1 atom/A^3, but we'll
00735 // be generous and allow much higher densities. 
00736 vec_sub(sidelen, max, min);
00737 // include estimate for atom radius (1 Angstrom) in volume determination
00738 volume = fabsf((sidelen[0] + 2.0f) * (sidelen[1] + 2.0f) * (sidelen[2] + 2.0f));
00739 if (isnan(volume)) {
00740 msgErr << "vmd_gridsearch3: bounding volume yielded NaN" << sendmsg;
00741 return NULL;
00742 }
00743 
00744 if (maxpairs != -1) {
00745 if (((numonA + numonB) / volume) > 1.0) {
00746 msgWarn << "vmd_gridsearch3: insane atom density" << sendmsg;
00747 }
00748 } 
00749 
00750 // I don't want the grid to get too large, otherwise I could run out
00751 // of memory. Octrees would be cool, but I'll just limit the grid size
00752 // and let the performance degrade a little for pathological systems.
00753 // Note that sqdist is what gets used for the actual distance checks;
00754 // from here on out pairdist is only used to set the grid size, so we 
00755 // can set it to anything larger than the original pairdist.
00756 const int MAXBOXES = 4000000;
00757 totb = MAXBOXES + 1;
00758 
00759 float newpairdist = pairdist;
00760 float xrange = max[0]-min[0];
00761 float yrange = max[1]-min[1];
00762 float zrange = max[2]-min[2];
00763 do {
00764 pairdist = newpairdist;
00765 const float invpairdist = 1.0f / pairdist;
00766 xb = ((int)(xrange*invpairdist))+1;
00767 yb = ((int)(yrange*invpairdist))+1;
00768 zb = ((int)(zrange*invpairdist))+1;
00769 xytotb = yb * xb;
00770 totb = xytotb * zb;
00771 newpairdist = pairdist * 1.26f; // cbrt(2) is about 1.26
00772 } while (totb > MAXBOXES || totb < 1); // check for integer wraparound too
00773 
00774 // 2. Sort each atom into appropriate bins
00775 boxatomA = (int **) calloc(1, totb*sizeof(int *));
00776 numinboxA = (int *) calloc(1, totb*sizeof(int));
00777 maxinboxA = (int *) calloc(1, totb*sizeof(int));
00778 if (boxatomA == NULL || numinboxA == NULL || maxinboxA == NULL) {
00779 if (boxatomA != NULL)
00780 free(boxatomA);
00781 if (numinboxA != NULL)
00782 free(numinboxA);
00783 if (maxinboxA != NULL)
00784 free(maxinboxA);
00785 msgErr << "Gridsearch memory allocation failed, bailing out" << sendmsg;
00786 return NULL; // ran out of memory, bail out!
00787 }
00788 
00789 const float invpairdist = 1.0f / pairdist;
00790 for (i=0; i<natomsA; i++) {
00791 if (A[i]) {
00792 int axb, ayb, azb, aindex, num;
00793 
00794 // compute box index for new atom
00795 const float *loc = posA + 3L*i;
00796 axb = (int)((loc[0] - min[0])*invpairdist);
00797 ayb = (int)((loc[1] - min[1])*invpairdist);
00798 azb = (int)((loc[2] - min[2])*invpairdist);
00799 
00800 // clamp box indices to valid range in case of FP error
00801 if (axb >= xb) axb = xb-1;
00802 if (ayb >= yb) ayb = yb-1;
00803 if (azb >= zb) azb = zb-1;
00804 
00805 aindex = azb * xytotb + ayb * xb + axb;
00806 
00807 // grow box if necessary
00808 if ((num = numinboxA[aindex]) == maxinboxA[aindex]) {
00809 boxatomA[aindex] = (int *) realloc(boxatomA[aindex], (num+4)*sizeof(int));
00810 maxinboxA[aindex] += 4;
00811 }
00812 
00813 // store atom index in box
00814 boxatomA[aindex][num] = i;
00815 numinboxA[aindex]++;
00816 }
00817 }
00818 free(maxinboxA);
00819 
00820 boxatomB = (int **) calloc(1, totb*sizeof(int *));
00821 numinboxB = (int *) calloc(1, totb*sizeof(int));
00822 maxinboxB = (int *) calloc(1, totb*sizeof(int));
00823 if (boxatomB == NULL || numinboxB == NULL || maxinboxB == NULL) {
00824 if (boxatomB != NULL)
00825 free(boxatomB);
00826 if (numinboxB != NULL)
00827 free(numinboxB);
00828 if (maxinboxB != NULL)
00829 free(maxinboxB);
00830 msgErr << "Gridsearch memory allocation failed, bailing out" << sendmsg;
00831 return NULL; // ran out of memory, bail out!
00832 }
00833 
00834 for (i=0; i<natomsB; i++) {
00835 if (B[i]) {
00836 int axb, ayb, azb, aindex, num;
00837 
00838 // compute box index for new atom
00839 const float *loc = posB + 3L*i;
00840 axb = (int)((loc[0] - min[0])*invpairdist);
00841 ayb = (int)((loc[1] - min[1])*invpairdist);
00842 azb = (int)((loc[2] - min[2])*invpairdist);
00843 
00844 // clamp box indices to valid range in case of FP error
00845 if (axb >= xb) axb = xb-1;
00846 if (ayb >= yb) ayb = yb-1;
00847 if (azb >= zb) azb = zb-1;
00848 
00849 aindex = azb * xytotb + ayb * xb + axb;
00850 
00851 // grow box if necessary
00852 if ((num = numinboxB[aindex]) == maxinboxB[aindex]) {
00853 boxatomB[aindex] = (int *) realloc(boxatomB[aindex], (num+4)*sizeof(int));
00854 maxinboxB[aindex] += 4;
00855 }
00856 
00857 // store atom index in box
00858 boxatomB[aindex][num] = i;
00859 numinboxB[aindex]++;
00860 }
00861 }
00862 free(maxinboxB);
00863 
00864 
00865 // 3. Build pairlists of atoms less than sqrtdist apart
00866 nbrlist = (int **) calloc(1, totb*sizeof(int *));
00867 if (make_neighborlist_sym(nbrlist, xb, yb, zb)) {
00868 if (boxatomA != NULL) {
00869 for (i=0; i<totb; i++) {
00870 if (boxatomA[i] != NULL) free(boxatomA[i]);
00871 }
00872 free(boxatomA);
00873 }
00874 if (boxatomB != NULL) {
00875 for (i=0; i<totb; i++) {
00876 if (boxatomB[i] != NULL) free(boxatomB[i]);
00877 }
00878 free(boxatomB);
00879 }
00880 if (nbrlist != NULL) {
00881 for (i=0; i<totb; i++) {
00882 if (nbrlist[i] != NULL) free(nbrlist[i]);
00883 }
00884 free(nbrlist);
00885 }
00886 free(numinboxA);
00887 free(numinboxB);
00888 msgErr << "Gridsearch memory allocation failed, bailing out" << sendmsg;
00889 return NULL; // ran out of memory, bail out!
00890 }
00891 
00892 // setup head of pairlist
00893 GridSearchPair *head, *cur;
00894 head = (GridSearchPair *) malloc(sizeof(GridSearchPair));
00895 head->next = NULL;
00896 cur = head;
00897 
00898 for (aindex = 0; aindex < totb; aindex++) {
00899 int *tmpbox, *tmpnbr, *nbr;
00900 tmpbox = boxatomA[aindex];
00901 tmpnbr = nbrlist[aindex];
00902 for (nbr = tmpnbr; (*nbr != -1) && (!maxpairsreached); nbr++) {
00903 int *nbrboxB = boxatomB[*nbr];
00904 for (i=0; (i<numinboxA[aindex]) && (!maxpairsreached); i++) {
00905 const float *p1;
00906 int ind1 = tmpbox[i];
00907 p1 = posA + 3L*ind1;
00908 for (j=0; (j<numinboxB[*nbr]) && (!maxpairsreached); j++) { 
00909 const float *p2;
00910 int ind2 = nbrboxB[j];
00911 p2 = posB + 3L*ind2;
00912 if (!allow_double_counting && B[ind1] && A[ind2] && ind2<=ind1) continue; //don't double-count bonds XXX
00913 if (distance2(p1,p2) > sqdist) continue;
00914 
00915 if (maxpairs > 0) {
00916 if (paircount >= maxpairs) {
00917 maxpairsreached = 1;
00918 continue; 
00919 } 
00920 }
00921 
00922 add_link(cur, ind1, ind2);
00923 paircount++;
00924 cur = cur->next;
00925 cur->next = NULL;
00926 }
00927 }
00928 }
00929 } 
00930 
00931 for (i=0; i<totb; i++) {
00932 free(boxatomA[i]);
00933 free(boxatomB[i]);
00934 free(nbrlist[i]);
00935 }
00936 free(boxatomA);
00937 free(boxatomB);
00938 free(nbrlist);
00939 free(numinboxA);
00940 free(numinboxB);
00941 
00942 cur = head->next;
00943 free(head);
00944 
00945 if (maxpairsreached) 
00946 msgErr << "vmdgridsearch3: exceeded pairlist sanity check, aborted" << sendmsg;
00947 
00948 return cur;
00949 }
00950 
00951 
00952 //
00953 // Multithreaded implementation of find_within() 
00954 //
00955 extern "C" void * find_within_routine( void *v );
00956 
00957 struct AtomEntry {
00958 float x, y, z;
00959 int index;
00960 AtomEntry() {}
00961 AtomEntry(const float &_x, const float &_y, const float &_z, const int &_i)
00962 : x(_x), y(_y), z(_z), index(_i) {}
00963 };
00964 
00965 typedef ResizeArray<AtomEntry> atomlist;
00966 struct FindWithinData {
00967 int nthreads;
00968 int tid;
00969 int totb;
00970 int xytotb;
00971 int xb;
00972 int yb;
00973 int zb;
00974 float r2;
00975 const float * xyz;
00976 const atomlist * flgatoms;
00977 const atomlist * otheratoms;
00978 int * flgs;
00979 FindWithinData() : flgatoms(0), otheratoms(0), flgs(0) {}
00980 ~FindWithinData() { if (flgs) free(flgs); }
00981 };
00982 
00983 #define MAXGRIDDIM 31 
00984 void find_within(const float *xyz, int *flgs, const int *others, 
00985 int num, float r) {
00986 int i;
00987 float xmin, xmax, ymin, ymax, zmin, zmax;
00988 float oxmin, oymin, ozmin, oxmax, oymax, ozmax;
00989 float xwidth, ywidth, zwidth; 
00990 const float *pos;
00991 
00992 // find min/max bounds of atom coordinates in flgs
00993 if (!find_minmax_selected(num, flgs, xyz, xmin, ymin, zmin, xmax, ymax, zmax) ||
00994 !find_minmax_selected(num, others, xyz, oxmin, oymin, ozmin, oxmax, oymax, ozmax)) {
00995 memset(flgs, 0, num*sizeof(int));
00996 return;
00997 }
00998 
00999 // Find the set of atoms with the smallest extent; here we use the sum
01000 // of the box dimensions though other choices might be better.
01001 float size = xmax+ymax+zmax - (xmin+ymin+zmin);
01002 float osize = oxmax+oymax+ozmax - (oxmin+oymin+ozmin);
01003 if (osize < size) {
01004 xmin=oxmin;
01005 ymin=oymin;
01006 zmin=ozmin;
01007 xmax=oxmax;
01008 ymax=oymax;
01009 zmax=ozmax;
01010 }
01011 
01012 // Generate a grid of mesh size r based on the computed size of the molecule.
01013 // We limit the size of the grid cell dimensions so that we don't get too
01014 // many grid cells.
01015 xwidth = (xmax-xmin)/(MAXGRIDDIM-1);
01016 if (xwidth < r) xwidth = r;
01017 ywidth = (ymax-ymin)/(MAXGRIDDIM-1);
01018 if (ywidth < r) ywidth = r;
01019 zwidth = (zmax-zmin)/(MAXGRIDDIM-1);
01020 if (zwidth < r) zwidth = r;
01021 
01022 // Adjust the bounds so that we include atoms that are in the outermost
01023 // grid cells.
01024 xmin -= xwidth;
01025 xmax += xwidth;
01026 ymin -= ywidth;
01027 ymax += ywidth;
01028 zmin -= zwidth;
01029 zmax += zwidth;
01030 
01031 // The number of grid cells needed in each dimension is 
01032 // (int)((xmax-xmin)/xwidth) + 1
01033 const int xb = (int)((xmax-xmin)/xwidth) + 1;
01034 const int yb = (int)((ymax-ymin)/ywidth) + 1;
01035 const int zb = (int)((zmax-zmin)/zwidth) + 1;
01036 
01037 int xytotb = yb * xb;
01038 int totb = xytotb * zb; 
01039 
01040 atomlist* flgatoms = new atomlist[totb];
01041 atomlist* otheratoms = new atomlist[totb];
01042 
01043 float ixwidth = 1.0f/xwidth;
01044 float iywidth = 1.0f/ywidth;
01045 float izwidth = 1.0f/zwidth;
01046 for (i=0; i<num; i++) {
01047 if (!flgs[i] && !others[i]) continue;
01048 pos=xyz+3L*i;
01049 float xi = pos[0];
01050 float yi = pos[1];
01051 float zi = pos[2];
01052 if (xi<xmin || xi>xmax || yi<ymin || yi>ymax || zi<zmin || zi>zmax) {
01053 continue;
01054 }
01055 AtomEntry entry(xi,yi,zi,i);
01056 int axb = (int)((xi - xmin)*ixwidth);
01057 int ayb = (int)((yi - ymin)*iywidth);
01058 int azb = (int)((zi - zmin)*izwidth);
01059 
01060 // Due to floating point error in the calcuation of bin widths, we 
01061 // have to range clamp the computed box indices.
01062 if (axb==xb) axb=xb-1;
01063 if (ayb==yb) ayb=yb-1;
01064 if (azb==zb) azb=zb-1;
01065 
01066 int aindex = azb*xytotb + ayb*xb + axb;
01067 if (others[i]) otheratoms[aindex].append(entry);
01068 if ( flgs[i]) flgatoms[aindex].append(entry);
01069 }
01070 
01071 memset(flgs, 0, num*sizeof(int));
01072 const float r2 = (float) (r*r); 
01073 
01074 // set up workspace for multithreaded calculation
01075 int nthreads;
01076 #ifdef VMDTHREADS
01077 nthreads = wkf_thread_numprocessors();
01078 wkf_thread_t * threads = (wkf_thread_t *)calloc(nthreads, sizeof(wkf_thread_t));
01079 #else
01080 nthreads = 1;
01081 #endif
01082 FindWithinData *data = new FindWithinData[nthreads];
01083 for (i=0; i<nthreads; i++) {
01084 data[i].nthreads = nthreads;
01085 data[i].tid = i;
01086 data[i].totb = totb;
01087 data[i].xytotb = xytotb;
01088 data[i].xb = xb;
01089 data[i].yb = yb;
01090 data[i].zb = zb;
01091 data[i].r2 = r2;
01092 data[i].xyz = xyz;
01093 data[i].flgatoms = flgatoms;
01094 data[i].otheratoms = otheratoms;
01095 data[i].flgs = (int *)calloc(num, sizeof(int));
01096 }
01097 #ifdef VMDTHREADS
01098 for (i=0; i<nthreads; i++) {
01099 wkf_thread_create(threads+i, find_within_routine, data+i);
01100 }
01101 for (i=0; i<nthreads; i++) {
01102 wkf_thread_join(threads[i], NULL);
01103 }
01104 free(threads);
01105 #else
01106 find_within_routine(data);
01107 #endif
01108 
01109 // combine results
01110 for (i=0; i<nthreads; i++) {
01111 const int *tflg = data[i].flgs;
01112 for (int j=0; j<num; j++) {
01113 flgs[j] |= tflg[j];
01114 }
01115 }
01116 delete [] data;
01117 delete [] flgatoms;
01118 delete [] otheratoms;
01119 }
01120 
01121 extern "C" void * find_within_routine( void *v ) {
01122 FindWithinData *data = (FindWithinData *)v;
01123 const int nthreads = data->nthreads;
01124 const int tid = data->tid;
01125 const int totb = data->totb;
01126 const int xytotb = data->xytotb;
01127 const int xb = data->xb;
01128 const int yb = data->yb;
01129 const int zb = data->zb;
01130 const float r2 = data->r2;
01131 const atomlist * flgatoms = data->flgatoms;
01132 const atomlist * otheratoms = data->otheratoms;
01133 int * flgs = data->flgs;
01134 
01135 // Loop over boxes, checking for flg atoms and other atoms within one
01136 // box of each other. When one is found, mark the flag.
01137 for (int aindex = tid; aindex<totb; aindex += nthreads) {
01138 // Figure out the neighbors for this box
01139 int zi = aindex/xytotb;
01140 int ytmp = aindex - zi*xytotb;
01141 int yi = ytmp/xb;
01142 int xi = ytmp - yi*xb;
01143 int nbrs[14];
01144 int n=0;
01145 nbrs[n++] = aindex; // Always include self
01146 if (xi < xb-1) nbrs[n++] = aindex + 1;
01147 if (yi < yb-1) nbrs[n++] = aindex + xb;
01148 if (zi < zb-1) nbrs[n++] = aindex + xytotb;
01149 if (xi < (xb-1) && yi < (yb-1)) nbrs[n++] = aindex + xb + 1;
01150 if (xi < (xb-1) && zi < (zb-1)) nbrs[n++] = aindex + xytotb + 1;
01151 if (yi < (yb-1) && zi < (zb-1)) nbrs[n++] = aindex + xytotb + xb;
01152 if (xi < (xb-1) && yi > 0) nbrs[n++] = aindex - xb + 1;
01153 if (xi > 0 && zi < (zb-1)) nbrs[n++] = aindex + xytotb - 1;
01154 if (yi > 0 && zi < (zb-1)) nbrs[n++] = aindex + xytotb - xb;
01155 if (xi < (xb-1) && yi < (yb-1) && zi < (zb-1))
01156 nbrs[n++] = aindex + xytotb + xb + 1;
01157 if (xi > 0 && yi < (yb-1) && zi < (zb-1))
01158 nbrs[n++] = aindex + xytotb + xb - 1;
01159 if (xi < (xb-1) && yi > 0 && zi < (zb-1))
01160 nbrs[n++] = aindex + xytotb - xb + 1;
01161 if (xi > 0 && yi > 0 && zi < (zb-1))
01162 nbrs[n++] = aindex + xytotb - xb - 1;
01163 
01164 const atomlist& boxflg = flgatoms[aindex];
01165 // Compare the atoms in boxflg to those in nbrother
01166 int i;
01167 for (i=0; i<boxflg.num(); i++) {
01168 const AtomEntry &flgentry = boxflg[i];
01169 int flgind = flgentry.index;
01170 
01171 // Compare flag atoms in this box to other atoms in neighbor boxes,
01172 for (int inbr=0; inbr<n; inbr++) { // Loop over neighbors
01173 if (flgs[flgind]) break;
01174 
01175 // Fetch a neighboring otheratoms to compare to boxflg
01176 int nbr = nbrs[inbr];
01177 const atomlist& nbrother = otheratoms[nbr];
01178 for (int j=0; j<nbrother.num(); j++) {
01179 const AtomEntry &otherentry = nbrother[j];
01180 float dx2 = flgentry.x - otherentry.x; dx2 *= dx2;
01181 float dy2 = flgentry.y - otherentry.y; dy2 *= dy2;
01182 float dz2 = flgentry.z - otherentry.z; dz2 *= dz2;
01183 if (dx2 + dy2 + dz2 < r2) {
01184 flgs[flgind] = 1;
01185 break;
01186 }
01187 }
01188 }
01189 }
01190 
01191 // compare other atoms in this box to flag atoms in the neighbors.
01192 const atomlist& boxother = otheratoms[aindex];
01193 
01194 for (int inbr=0; inbr<n; inbr++) { // Loop over neighbors
01195 int nbr = nbrs[inbr];
01196 
01197 // Fetch a neighboring flgatoms to compare to boxother 
01198 const atomlist& nbrflg = flgatoms[nbr];
01199 
01200 // Compare the atoms in boxother to those in nbrflg
01201 for (i=0; i<nbrflg.num(); i++) {
01202 const AtomEntry &flgentry = nbrflg[i];
01203 int flgind = flgentry.index;
01204 // The next test helps a lot when the boxes are large, but hurts
01205 // a little when the boxes are small.
01206 if (flgs[flgind]) continue;
01207 for (int j=0; j<boxother.num(); j++) {
01208 const AtomEntry &otherentry = boxother[j];
01209 float dx2 = flgentry.x - otherentry.x; dx2 *= dx2;
01210 float dy2 = flgentry.y - otherentry.y; dy2 *= dy2;
01211 float dz2 = flgentry.z - otherentry.z; dz2 *= dz2;
01212 if (dx2 + dy2 + dz2 < r2) {
01213 flgs[flgind] = 1;
01214 break;
01215 }
01216 }
01217 }
01218 }
01219 } 
01220 return NULL;
01221 } 
01222 
01223