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

BondSearch.C

Go to the documentation of this file.
00001 /***************************************************************************
00002  *cr                                                                       
00003  *cr            (C) Copyright 1995-2011 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: BondSearch.C,v $
00013  *      $Author: johns $        $Locker:  $             $State: Exp $
00014  *      $Revision: 1.71 $       $Date: 2013/01/30 21:01:20 $
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 
00037 GridSearchPairlist *vmd_gridsearch_bonds(const float *pos, const float *radii,
00038                                    int natoms, float pairdist, int maxpairs) {
00039   float min[3], max[3];
00040   int i, xb, yb, zb, xytotb, totb;
00041   int **boxatom, *numinbox, *maxinbox, **nbrlist;
00042   int numon = 0;
00043   float sidelen[3], volume;
00044   int paircount = 0;
00045 
00046   // find bounding box for selected atoms, and number of atoms in selection.
00047 #if 1
00048   minmax_3fv_aligned(pos, natoms, min, max);
00049 #else
00050   find_minmax_all(pos, natoms, min, max);
00051 #endif
00052 
00053   // do sanity checks and complain if we've got bogus atom coordinates,
00054   // we shouldn't ever have density higher than 0.1 atom/A^3, but we'll
00055   // be generous and allow much higher densities.  
00056   if (maxpairs != -1) {
00057     vec_sub(sidelen, max, min);
00058     // include estimate for atom radius (1 Angstrom) in volume determination
00059     volume = fabsf((sidelen[0] + 2.0f) * (sidelen[1] + 2.0f) * (sidelen[2] + 2.0f));
00060     if ((numon / volume) > 1.0) {
00061       msgWarn << "vmd_gridsearch_bonds: insane atom density" << sendmsg;
00062     }
00063   }
00064 
00065   // I don't want the grid to get too large, otherwise I could run out
00066   // of memory.  Octrees would be cool, but I'll just limit the grid size
00067   // and let the performance degrade a little for pathological systems.
00068   // Note that pairdist^2 is what gets used for the actual distance checks;
00069   // from here on out pairdist is only used to set the grid size, so we 
00070   // can set it to anything larger than the original pairdist.
00071   const int MAXBOXES = 4000000;
00072   totb = MAXBOXES + 1;
00073 
00074   float newpairdist = pairdist;
00075   float xrange = max[0]-min[0];
00076   float yrange = max[1]-min[1];
00077   float zrange = max[2]-min[2];
00078   do {
00079     pairdist = newpairdist;
00080     const float invpairdist = 1.0f / pairdist; 
00081     xb = ((int)(xrange*invpairdist))+1;
00082     yb = ((int)(yrange*invpairdist))+1;
00083     zb = ((int)(zrange*invpairdist))+1;
00084     xytotb = yb * xb;
00085     totb = xytotb * zb;
00086     newpairdist = pairdist * 1.26f; // cbrt(2) is about 1.26
00087   } while (totb > MAXBOXES || totb < 1); // check for integer wraparound too
00088  
00089   // 2. Sort each atom into appropriate bins
00090   boxatom = (int **) calloc(1, totb*sizeof(int *));
00091   numinbox = (int *) calloc(1, totb*sizeof(int));
00092   maxinbox = (int *) calloc(1, totb*sizeof(int));
00093   if (boxatom == NULL || numinbox == NULL || maxinbox == NULL) {
00094     if (boxatom != NULL)
00095       free(boxatom);
00096     if (numinbox != NULL)
00097       free(numinbox);
00098     if (maxinbox != NULL)
00099       free(maxinbox);
00100     msgErr << "Bondsearch memory allocation failed, bailing out" << sendmsg;
00101     return NULL; // ran out of memory, bail out!
00102   }
00103 
00104   const float invpairdist = 1.0f / pairdist; 
00105   for (i=0; i<natoms; i++) {
00106     int axb, ayb, azb, aindex, num;
00107 
00108     // compute box index for new atom
00109     const float *loc = pos + 3*i;
00110     axb = (int)((loc[0] - min[0])*invpairdist);
00111     ayb = (int)((loc[1] - min[1])*invpairdist);
00112     azb = (int)((loc[2] - min[2])*invpairdist);
00113 
00114     // clamp box indices to valid range in case of FP error
00115     if (axb >= xb) axb = xb-1;
00116     if (ayb >= yb) ayb = yb-1;
00117     if (azb >= zb) azb = zb-1;
00118 
00119     aindex = azb * xytotb + ayb * xb + axb;
00120 
00121     // grow box if necessary 
00122     if ((num = numinbox[aindex]) == maxinbox[aindex]) {
00123       boxatom[aindex] = (int *) realloc(boxatom[aindex], (num+4)*sizeof(int));
00124       maxinbox[aindex] += 4;
00125     }
00126 
00127     // store atom index in box
00128     boxatom[aindex][num] = i;
00129     numinbox[aindex]++;
00130   }
00131   free(maxinbox);
00132  
00133   nbrlist = (int **) calloc(1, totb*sizeof(int *));
00134   if (make_neighborlist(nbrlist, xb, yb, zb)) {
00135     if (boxatom != NULL) {
00136       for (i=0; i<totb; i++) {
00137         if (boxatom[i] != NULL) free(boxatom[i]);
00138       }
00139       free(boxatom);
00140     }
00141     if (nbrlist != NULL) {
00142       for (i=0; i<totb; i++) {
00143         if (nbrlist[i] != NULL) free(nbrlist[i]);
00144       }
00145       free(nbrlist);
00146     }
00147     free(numinbox);
00148     msgErr << "Bondsearch memory allocation failed, bailing out" << sendmsg;
00149     return NULL; // ran out of memory, bail out!
00150   }
00151 
00152   // if maxpairs is "unlimited", set it to the biggest positive int
00153   if (maxpairs < 0) {
00154     maxpairs = 2147483647;
00155   }
00156 
00157   // setup head of pairlist
00158   GridSearchPairlist *head, *cur;
00159   head = (GridSearchPairlist *) malloc(sizeof(GridSearchPairlist));
00160   head->next = NULL;
00161   paircount = vmd_bondsearch_thr(pos, radii, head, totb, 
00162                                  boxatom, numinbox, nbrlist, 
00163                                  maxpairs, pairdist);
00164 
00165   for (i=0; i<totb; i++) {
00166     free(boxatom[i]);
00167     free(nbrlist[i]);
00168   }
00169   free(boxatom);
00170   free(nbrlist);
00171   free(numinbox);
00172 
00173   cur = head->next;
00174   free(head);
00175 
00176   if (paircount > maxpairs) 
00177     msgErr << "vmdgridsearch_bonds: exceeded pairlist sanity check, aborted" << sendmsg;
00178 
00179   return cur;
00180 }
00181 
00182 
00183 
00184 // bond search thread parameter structure
00185 typedef struct {
00186   int threadid;
00187   int threadcount;
00188   wkf_mutex_t *pairlistmutex;
00189   GridSearchPairlist * head;
00190   float *pos;
00191   float *radii;
00192   int totb;
00193   int **boxatom;
00194   int *numinbox;
00195   int **nbrlist;  
00196   int maxpairs;
00197   float pairdist;
00198 } bondsearchthrparms;
00199 
00200 // thread prototype 
00201 extern "C" void * bondsearchthread(void *);
00202 
00203 // setup and launch bond search threads
00204 int vmd_bondsearch_thr(const float *pos, const float *radii,
00205                        GridSearchPairlist * head, 
00206                        int totb, int **boxatom, 
00207                        int *numinbox, int **nbrlist, int maxpairs, 
00208                        float pairdist) {
00209   int i;
00210   bondsearchthrparms *parms;
00211   wkf_thread_t * threads;
00212   wkf_mutex_t pairlistmutex; 
00213   wkf_mutex_init(&pairlistmutex); // init mutex before use
00214 
00215   int numprocs = wkf_thread_numprocessors();
00216 
00217   /* allocate array of threads */
00218   threads = (wkf_thread_t *) calloc(numprocs * sizeof(wkf_thread_t), 1);
00219 
00220   /* allocate and initialize array of thread parameters */
00221   parms = (bondsearchthrparms *) malloc(numprocs * sizeof(bondsearchthrparms));
00222   for (i=0; i<numprocs; i++) {
00223     parms[i].threadid = i;
00224     parms[i].threadcount = numprocs;
00225     parms[i].pairlistmutex = &pairlistmutex;
00226     parms[i].head = NULL;
00227     parms[i].pos = (float *) pos;
00228     parms[i].radii = (float *) radii;
00229     parms[i].totb = totb;
00230     parms[i].boxatom = boxatom;
00231     parms[i].numinbox = numinbox;
00232     parms[i].nbrlist = nbrlist;  
00233     parms[i].maxpairs = maxpairs;
00234     parms[i].pairdist = pairdist;
00235   }
00236 
00237 #if defined(VMDTHREADS)
00238   /* spawn child threads to do the work */
00239   for (i=0; i<numprocs; i++) {
00240     wkf_thread_create(&threads[i], bondsearchthread, &parms[i]);
00241   }
00242 
00243   /* join the threads after work is done */
00244   for (i=0; i<numprocs; i++) {
00245     wkf_thread_join(threads[i], NULL);
00246   }
00247 #else
00248   bondsearchthread(&parms[0]); // single-threaded code
00249 #endif
00250 
00251   // assemble final pairlist from sublists
00252   for (i=0; i<numprocs; i++) {
00253     if (parms[i].head != NULL) {
00254       GridSearchPairlist *tmp = head->next;
00255       head->next = parms[i].head;
00256       parms[i].head->next = tmp;
00257     }
00258   }
00259 
00260   wkf_mutex_destroy(&pairlistmutex); // destroy mutex when finished
00261 
00262   /* free thread parms */
00263   free(parms);
00264   free(threads);
00265 
00266   return 0;
00267 }
00268 
00269 extern "C" void * bondsearchthread(void *voidparms) {
00270   int i, j, aindex;
00271   int paircount = 0;
00272 
00273   bondsearchthrparms *parms = (bondsearchthrparms *) voidparms;
00274 
00275   const int threadid = parms->threadid;
00276   const int threadcount = parms->threadcount;
00277   wkf_mutex_t *pairlistmutex = parms->pairlistmutex;
00278   const float *pos = parms->pos;
00279   const float *radii = parms->radii;
00280   const int totb = parms->totb;
00281   const int **boxatom = (const int **) parms->boxatom;
00282   const int *numinbox = parms->numinbox;
00283   const int **nbrlist = (const int **) parms->nbrlist; 
00284   const int maxpairs = parms->maxpairs;
00285   const float pairdist = parms->pairdist;
00286 
00287   ResizeArray<int> *pairs = new ResizeArray<int>;
00288   float sqdist = pairdist * pairdist;
00289 
00290   wkfmsgtimer *msgt = wkf_msg_timer_create(5);
00291   for (aindex = threadid; (aindex < totb) && (paircount < maxpairs); aindex+=threadcount) {
00292     const int *tmpbox, *nbr;
00293     tmpbox = boxatom[aindex];
00294     int anbox = numinbox[aindex];
00295 
00296     if (((aindex & 255) == 0) && wkf_msg_timer_timeout(msgt)) {
00297       char tmpbuf[128];
00298       sprintf(tmpbuf, "%6.2f", (100.0f * aindex) / (float) totb);
00299 //  XXX: we have to use printf here as msgInfo is not thread-safe.
00300 //  msgInfo << "vmd_gridsearch_bonds (thread " << threadid << "): " 
00301 //          << tmpbuf << "% complete" << sendmsg;
00302       printf("vmd_gridsearch_bonds (thread %d): %s %% complete\n", 
00303         threadid, tmpbuf); 
00304     }
00305 
00306     for (nbr = nbrlist[aindex]; (*nbr != -1) && (paircount < maxpairs); nbr++) {
00307       int nnbr=*nbr;
00308       const int *nbrbox = boxatom[nnbr];
00309       int nbox=numinbox[nnbr];
00310       int self = (aindex == nnbr) ? 1 : 0;
00311 
00312       for (i=0; (i<anbox) && (paircount < maxpairs); i++) {
00313         int ind1 = tmpbox[i];
00314         const float *p1 = pos + 3*ind1;
00315 
00316         // skip over self and already-tested atoms
00317         int startj = (self) ? i+1 : 0;
00318 
00319         for (j=startj; (j<nbox) && (paircount < maxpairs); j++) {
00320           int ind2 = nbrbox[j];
00321           const float *p2 = pos + 3*ind2;
00322           float dx = p1[0] - p2[0];
00323           float dy = p1[1] - p2[1];
00324           float dz = p1[2] - p2[2];
00325           float ds2 = dx*dx + dy*dy + dz*dz;
00326 
00327           // perform distance test, but ignore pairs between atoms
00328           // with nearly identical coords
00329           if ((ds2 > sqdist) || (ds2 < 0.001))
00330             continue;
00331 
00332           if (radii) { // Do atom-specific distance check
00333             float cut = 0.6f * (radii[ind1] + radii[ind2]);
00334             if (ds2 > cut*cut)
00335               continue;
00336           }
00337 
00338           pairs->append(ind1);
00339           pairs->append(ind2);
00340           paircount++;
00341         }
00342       }
00343     }
00344   }
00345 
00346   // setup results pairlist node
00347   GridSearchPairlist *head;
00348   head = (GridSearchPairlist *) malloc(sizeof(GridSearchPairlist));
00349   head->next = NULL;
00350   head->pairlist = pairs;
00351 
00352   wkf_mutex_lock(pairlistmutex);   // lock pairlist before update
00353   parms->head = head;
00354   wkf_mutex_unlock(pairlistmutex); // unlock pairlist after update
00355 
00356   wkf_msg_timer_destroy(msgt);
00357 
00358   return NULL;
00359 }
00360 
00361 
00362 
00363 
00364 
00365 // determine bonds from position of atoms previously read.
00366 // If cutoff < 0, use vdw radius to determine if bonded.
00367 int vmd_bond_search(BaseMolecule *mol, const Timestep *ts, 
00368                     float cutoff, int dupcheck) {
00369   const float *pos;
00370   int natoms;
00371   int i;
00372   const float *radius = mol->radius();
00373  
00374   if (ts == NULL) {
00375     msgErr << "Internal Error: NULL Timestep in vmd_bond_search" << sendmsg;
00376     return 0;
00377   }
00378 
00379   natoms = mol->nAtoms; 
00380   if (natoms == 0 || cutoff == 0.0)
00381     return 1;
00382 
00383   msgInfo << "Determining bond structure from distance search ..." << sendmsg;
00384 
00385   if (dupcheck)
00386     msgInfo << "Eliminating bonds duplicated from existing structure..." << sendmsg;
00387 
00388   // Set box distance to either the cutoff, or 1.2 times the largest VDW radius
00389   float dist = cutoff; 
00390   if (cutoff < 0.0) {
00391     // set minimum cutoff distance for the case when the loaded molecule
00392     // has radii set to zero.  This must be >0.0 or the grid search will hang.
00393     dist = 0.833f;
00394     for (i=0; i<natoms; i++) {  
00395       float rad = radius[i];
00396       if (rad > dist) dist = rad;
00397     }
00398     dist = 1.2f * dist;
00399   }
00400   
00401   pos = ts->pos; 
00402  
00403   // Call the bond search to get all atom pairs within the specified distance
00404   // XXX set maxpairs to 27 bonds per atom, which ought to be ridiculously high
00405   //     for any real structure someone would load, but well below N^2
00406   GridSearchPairlist *pairlist = vmd_gridsearch_bonds(pos, 
00407                                                  (cutoff < 0) ? radius : NULL,
00408                                                  natoms, dist, natoms * 27);
00409 
00410   // Go through the pairlist adding validated bonds freeing nodes as we go. 
00411   GridSearchPairlist *p, *tmp; 
00412   for (p = pairlist; p != NULL; p = tmp) {
00413     int numpairs = p->pairlist->num() / 2;
00414     
00415     for (i=0; i<numpairs; i++) {
00416       int ind1 = (*p->pairlist)[i*2  ]; 
00417       int ind2 = (*p->pairlist)[i*2+1];
00418 
00419       MolAtom *atom1 = mol->atom(ind1);
00420       MolAtom *atom2 = mol->atom(ind2);
00421 
00422       // don't bond atoms that aren't part of the same conformation
00423       // or that aren't in the all-conformations part of the structure
00424       if ((atom1->altlocindex != atom2->altlocindex) &&
00425           ((mol->altlocNames.name(atom1->altlocindex)[0] != '\0') && 
00426            (mol->altlocNames.name(atom2->altlocindex)[0] != '\0'))) {
00427         continue;
00428       }
00429 
00430       // Prevent hydrogens from bonding with each other.
00431 #if 1
00432       // XXX must use the atom name strings presently because the
00433       // hydrogen flags aren't necessarily set by the time the bond search
00434       // code executes.  It may soon be time to do something a different
00435       // with per-atom flag storage so that the atom types can be setup
00436       // during the earliest phase of structure analysis, eliminating this
00437       // and other potential gotchas.
00438       if (!IS_HYDROGEN(mol->atomNames.name(atom1->nameindex)) ||
00439           !IS_HYDROGEN(mol->atomNames.name(atom2->nameindex)) ) {
00440 #else
00441       // Use atomType info derived during initial molecule analysis for speed.
00442       if (!(atom1->atomType == ATOMHYDROGEN) ||
00443           !(atom2->atomType == ATOMHYDROGEN)) {
00444 #endif
00445         // Add a bond, bondorder defaults to 1, bond type to -1
00446         if (dupcheck)
00447           mol->add_bond_dupcheck(ind1, ind2, 1, -1);
00448         else
00449           mol->add_bond(ind1, ind2, 1, -1);
00450       }
00451     }
00452 
00453     // free this pairlist node and its ResizeArray of pairs
00454     tmp = p->next;
00455     delete p->pairlist;
00456     free(p);
00457   }
00458 
00459   return 1;
00460 }
00461 
00462 
00463 
00464 

Generated on Sun May 19 01:50:09 2013 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002