00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
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 "VMDThreads.h"
00032 #include <ctype.h>
00033 #include <string.h>
00034
00035
00036 GridSearchPairlist *vmd_gridsearch_bonds(const float *pos, const float *radii,
00037 int natoms, float pairdist, int maxpairs) {
00038 GridSearchPairlist *head, *cur;
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 const float *loc;
00045 int paircount = 0;
00046
00047
00048 find_minmax_all(pos, natoms, min, max);
00049
00050
00051
00052
00053 if (maxpairs != -1) {
00054 vec_sub(sidelen, max, min);
00055
00056 volume = fabsf((sidelen[0] + 2.0f) * (sidelen[1] + 2.0f) * (sidelen[2] + 2.0f));
00057 if ((numon / volume) > 1.0) {
00058 msgWarn << "vmd_gridsearch_bonds: insane atom density" << sendmsg;
00059 }
00060 }
00061
00062
00063
00064
00065
00066
00067
00068 const int MAXBOXES = 4000000;
00069 totb = MAXBOXES + 1;
00070
00071 float newpairdist = pairdist;
00072 do {
00073 pairdist = newpairdist;
00074 xb = (int)((max[0]-min[0])/pairdist)+1;
00075 yb = (int)((max[1]-min[1])/pairdist)+1;
00076 zb = (int)((max[2]-min[2])/pairdist)+1;
00077 xytotb = yb * xb;
00078 totb = xytotb * zb;
00079 newpairdist = pairdist * 1.26f;
00080 } while (totb > MAXBOXES || totb < 1);
00081
00082
00083 boxatom = (int **) calloc(1, totb*sizeof(int *));
00084 numinbox = (int *) calloc(1, totb*sizeof(int));
00085 maxinbox = (int *) calloc(1, totb*sizeof(int));
00086 if (boxatom == NULL || numinbox == NULL || maxinbox == NULL) {
00087 if (boxatom != NULL)
00088 free(boxatom);
00089 if (numinbox != NULL)
00090 free(numinbox);
00091 if (maxinbox != NULL)
00092 free(maxinbox);
00093 msgErr << "Gridsearch memory allocation failed, bailing out" << sendmsg;
00094 return NULL;
00095 }
00096
00097 loc = pos;
00098 const float invpairdist = 1.0f / pairdist;
00099 for (i=0; i<natoms; i++) {
00100 int axb, ayb, azb, aindex, num;
00101
00102 axb = (int)((loc[0] - min[0])*invpairdist);
00103 ayb = (int)((loc[1] - min[1])*invpairdist);
00104 azb = (int)((loc[2] - min[2])*invpairdist);
00105 aindex = azb * xytotb + ayb * xb + axb;
00106 num = numinbox[aindex];
00107 if (maxinbox[aindex] == num) {
00108 void *tmp = realloc(boxatom[aindex], (num+4)*sizeof(int));
00109 boxatom[aindex] = (int *)tmp;
00110 maxinbox[aindex] += 4;
00111 }
00112 boxatom[aindex][num] = i;
00113 numinbox[aindex]++;
00114
00115 loc += 3;
00116 }
00117 free(maxinbox);
00118
00119 nbrlist = (int **) calloc(1, totb*sizeof(int *));
00120 if (make_neighborlist(nbrlist, xb, yb, zb)) {
00121 if (boxatom != NULL) {
00122 for (i=0; i<totb; i++) {
00123 if (boxatom[i] != NULL) free(boxatom[i]);
00124 }
00125 free(boxatom);
00126 }
00127 if (nbrlist != NULL) {
00128 for (i=0; i<totb; i++) {
00129 if (nbrlist[i] != NULL) free(nbrlist[i]);
00130 }
00131 free(nbrlist);
00132 }
00133 free(numinbox);
00134 msgErr << "Gridsearch memory allocation failed, bailing out" << sendmsg;
00135 return NULL;
00136 }
00137
00138
00139 if (maxpairs < 0) {
00140 maxpairs = 2147483647;
00141 }
00142
00143
00144 head = (GridSearchPairlist *) malloc(sizeof(GridSearchPairlist));
00145 head->next = NULL;
00146 paircount = vmd_bondsearch_thr(pos, radii, head, totb,
00147 boxatom, numinbox, nbrlist,
00148 maxpairs, pairdist);
00149
00150 for (i=0; i<totb; i++) {
00151 free(boxatom[i]);
00152 free(nbrlist[i]);
00153 }
00154 free(boxatom);
00155 free(nbrlist);
00156 free(numinbox);
00157
00158 cur = head->next;
00159 free(head);
00160
00161 if (paircount > maxpairs)
00162 msgErr << "vmdgridsearch_bonds: exceeded pairlist sanity check, aborted" << sendmsg;
00163
00164 return cur;
00165 }
00166
00167
00168
00169
00170 typedef struct {
00171 int threadid;
00172 int threadcount;
00173 vmd_mutex_t *pairlistmutex;
00174 GridSearchPairlist * head;
00175 float *pos;
00176 float *radii;
00177 int totb;
00178 int **boxatom;
00179 int *numinbox;
00180 int **nbrlist;
00181 int maxpairs;
00182 float pairdist;
00183 } bondsearchthrparms;
00184
00185
00186 static void * bondsearchthread(void *);
00187
00188
00189 int vmd_bondsearch_thr(const float *pos, const float *radii,
00190 GridSearchPairlist * head,
00191 int totb, int **boxatom,
00192 int *numinbox, int **nbrlist, int maxpairs,
00193 float pairdist) {
00194 int i;
00195 bondsearchthrparms *parms;
00196 vmd_thread_t * threads;
00197 vmd_mutex_t pairlistmutex;
00198 vmd_mutex_init(&pairlistmutex);
00199
00200 int numprocs = vmd_thread_numprocessors();
00201
00202
00203 threads = (vmd_thread_t *) calloc(numprocs * sizeof(vmd_thread_t), 1);
00204
00205
00206 parms = (bondsearchthrparms *) malloc(numprocs * sizeof(bondsearchthrparms));
00207 for (i=0; i<numprocs; i++) {
00208 parms[i].threadid = i;
00209 parms[i].threadcount = numprocs;
00210 parms[i].pairlistmutex = &pairlistmutex;
00211 parms[i].head = NULL;
00212 parms[i].pos = (float *) pos;
00213 parms[i].radii = (float *) radii;
00214 parms[i].totb = totb;
00215 parms[i].boxatom = boxatom;
00216 parms[i].numinbox = numinbox;
00217 parms[i].nbrlist = nbrlist;
00218 parms[i].maxpairs = maxpairs;
00219 parms[i].pairdist = pairdist;
00220 }
00221
00222 #if defined(VMDTHREADS)
00223
00224 for (i=0; i<numprocs; i++) {
00225 vmd_thread_create(&threads[i], bondsearchthread, &parms[i]);
00226 }
00227
00228
00229 for (i=0; i<numprocs; i++) {
00230 vmd_thread_join(threads[i], NULL);
00231 }
00232 #else
00233 bondsearchthread(&parms[0]);
00234 #endif
00235
00236
00237 for (i=0; i<numprocs; i++) {
00238 if (parms[i].head != NULL) {
00239 GridSearchPairlist *tmp = head->next;
00240 head->next = parms[i].head;
00241 parms[i].head->next = tmp;
00242 }
00243 }
00244
00245 vmd_mutex_destroy(&pairlistmutex);
00246
00247
00248 free(parms);
00249 free(threads);
00250
00251 return 0;
00252 }
00253
00254 static void * bondsearchthread(void *voidparms) {
00255 int i, j, aindex;
00256 int paircount = 0;
00257
00258 bondsearchthrparms *parms = (bondsearchthrparms *) voidparms;
00259
00260 const int threadid = parms->threadid;
00261 const int threadcount = parms->threadcount;
00262 vmd_mutex_t *pairlistmutex = parms->pairlistmutex;
00263 const float *pos = parms->pos;
00264 const float *radii = parms->radii;
00265 const int totb = parms->totb;
00266 const int **boxatom = (const int **) parms->boxatom;
00267 const int *numinbox = parms->numinbox;
00268 const int **nbrlist = (const int **) parms->nbrlist;
00269 const int maxpairs = parms->maxpairs;
00270 const float pairdist = parms->pairdist;
00271
00272 ResizeArray<int> *pairs = new ResizeArray<int>;
00273 float sqdist = pairdist * pairdist;
00274
00275 msgtimer *msgt = msg_timer_create(5);
00276 for (aindex = threadid; (aindex < totb) && (paircount < maxpairs); aindex+=threadcount) {
00277 const int *tmpbox, *nbr;
00278 tmpbox = boxatom[aindex];
00279 int anbox = numinbox[aindex];
00280
00281 if (((aindex & 255) == 0) && msg_timer_timeout(msgt)) {
00282 char tmpbuf[128];
00283 sprintf(tmpbuf, "%6.2f", (100.0f * aindex) / (float) totb);
00284
00285
00286
00287 printf("vmd_gridsearch_bonds (thread %d): %s %% complete\n",
00288 threadid, tmpbuf);
00289 }
00290
00291 for (nbr = nbrlist[aindex]; (*nbr != -1) && (paircount < maxpairs); nbr++) {
00292 int nnbr=*nbr;
00293 const int *nbrbox = boxatom[nnbr];
00294 int nbox=numinbox[nnbr];
00295 int self = (aindex == nnbr) ? 1 : 0;
00296
00297 for (i=0; (i<anbox) && (paircount < maxpairs); i++) {
00298 int ind1 = tmpbox[i];
00299 const float *p1 = pos + 3*ind1;
00300
00301
00302 int startj = (self) ? i+1 : 0;
00303
00304 for (j=startj; (j<nbox) && (paircount < maxpairs); j++) {
00305 int ind2 = nbrbox[j];
00306 const float *p2 = pos + 3*ind2;
00307 float dx = p1[0] - p2[0];
00308 float dy = p1[1] - p2[1];
00309 float dz = p1[2] - p2[2];
00310 float ds2 = dx*dx + dy*dy + dz*dz;
00311
00312
00313
00314 if ((ds2 > sqdist) || (ds2 < 0.001))
00315 continue;
00316
00317 if (radii) {
00318 float cut = 0.6f * (radii[ind1] + radii[ind2]);
00319 if (ds2 > cut*cut)
00320 continue;
00321 }
00322
00323 pairs->append(ind1);
00324 pairs->append(ind2);
00325 paircount++;
00326 }
00327 }
00328 }
00329 }
00330
00331
00332 GridSearchPairlist *head;
00333 head = (GridSearchPairlist *) malloc(sizeof(GridSearchPairlist));
00334 head->next = NULL;
00335 head->pairlist = pairs;
00336
00337 vmd_mutex_lock(pairlistmutex);
00338 parms->head = head;
00339 vmd_mutex_unlock(pairlistmutex);
00340
00341 msg_timer_destroy(msgt);
00342
00343 return NULL;
00344 }
00345
00346
00347
00348
00349
00350
00351
00352 int vmd_bond_search(BaseMolecule *mol, const Timestep *ts,
00353 float cutoff, int dupcheck) {
00354 const float *pos;
00355 int natoms;
00356 int i;
00357 const float *radius = mol->radius();
00358
00359 if (ts == NULL) {
00360 msgErr << "Internal Error: NULL Timestep in vmd_bond_search" << sendmsg;
00361 return 0;
00362 }
00363
00364 natoms = mol->nAtoms;
00365 if (natoms == 0 || cutoff == 0.0)
00366 return 1;
00367
00368 msgInfo << "Determining bond structure from distance search ..." << sendmsg;
00369
00370 if (dupcheck)
00371 msgInfo << "Eliminating bonds duplicated from existing structure..." << sendmsg;
00372
00373
00374 float dist = cutoff;
00375 if (cutoff < 0.0) {
00376 for (i=0; i<natoms; i++) {
00377 float rad = radius[i];
00378 if (rad > dist) dist = rad;
00379 }
00380 dist = 1.2f * dist;
00381 }
00382
00383 pos = ts->pos;
00384
00385
00386
00387
00388 GridSearchPairlist *pairlist = vmd_gridsearch_bonds(pos,
00389 (cutoff < 0) ? radius : NULL,
00390 natoms, dist, natoms * 27);
00391
00392
00393 GridSearchPairlist *p, *tmp;
00394 for (p = pairlist; p != NULL; p = tmp) {
00395 int numpairs = p->pairlist->num() / 2;
00396
00397 for (i=0; i<numpairs; i++) {
00398 int ind1 = (*p->pairlist)[i*2 ];
00399 int ind2 = (*p->pairlist)[i*2+1];
00400
00401 MolAtom *atom1 = mol->atom(ind1);
00402 MolAtom *atom2 = mol->atom(ind2);
00403
00404
00405
00406 if ((atom1->altlocindex != atom2->altlocindex) &&
00407 ((mol->altlocNames.name(atom1->altlocindex)[0] != '\0') &&
00408 (mol->altlocNames.name(atom2->altlocindex)[0] != '\0'))) {
00409 continue;
00410 }
00411
00412
00413 #if 1
00414
00415
00416
00417
00418
00419
00420 if (!IS_HYDROGEN(mol->atomNames.name(atom1->nameindex)) ||
00421 !IS_HYDROGEN(mol->atomNames.name(atom2->nameindex)) ) {
00422 #else
00423
00424 if (!(atom1->atomType == ATOMHYDROGEN) ||
00425 !(atom2->atomType == ATOMHYDROGEN)) {
00426 #endif
00427
00428 if (dupcheck)
00429 mol->add_bond_dupcheck(ind1, ind2, 1);
00430 else
00431 mol->add_bond(ind1, ind2, 1);
00432 }
00433 }
00434
00435
00436 tmp = p->next;
00437 delete p->pairlist;
00438 free(p);
00439 }
00440
00441 return 1;
00442 }
00443
00444
00445
00446