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

Generated on Sat May 26 01:47:44 2012 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002