NAMD: msm_shortrng.c Source File

00001 
00029 #include "msm_defn.h"
00030 
00031 static int setup_bin_data(NL_Msm *pm);
00032 static int spatial_hashing(NL_Msm *pm);
00033 static int bin_evaluation_1away(NL_Msm *pm);
00034 static int bin_evaluation_k_away(NL_Msm *pm);
00035 
00036 
00037 int NL_msm_compute_short_range(NL_Msm *pm) {
00038   int rc = 0;  /* return code */
00039 
00040   if (pm->maxatoms < pm->numatoms) {
00041     rc = setup_bin_data(pm);
00042     if (rc) return rc;
00043   }
00044   rc = spatial_hashing(pm);
00045   if (rc) return rc;
00046   if (pm->msmflags & NL_MSM_COMPUTE_1AWAY) {
00047     rc = bin_evaluation_1away(pm);
00048   }
00049   else {
00050     rc = bin_evaluation_k_away(pm);
00051   }
00052   if (rc) return rc;
00053   return NL_MSM_SUCCESS;
00054 }
00055 
00056 
00057 int setup_bin_data(NL_Msm *pm) {
00058   void *vptr = NULL;
00059 
00060   ASSERT(pm->maxatoms < pm->numatoms);
00061   vptr = malloc(pm->numatoms * sizeof(int));
00062   if (vptr == NULL) return NL_MSM_ERROR_MALLOC;
00063   pm->next = (int *) vptr;
00064   pm->maxatoms = pm->numatoms;
00065   return NL_MSM_SUCCESS;
00066 }
00067 
00068 
00069 int spatial_hashing(NL_Msm *pm) {
00070   double d[3], s[3];
00071   const double *atom = pm->atom;  /* stored x/y/z/q */
00072   const double *c = pm->cellcenter;
00073   const double *ru = pm->recipvec1;
00074   const double *rv = pm->recipvec2;
00075   const double *rw = pm->recipvec3;
00076   int *bin = pm->bin;
00077   int *next = pm->next;
00078   int numbins = pm->numbins;
00079   int nbx = pm->nbx;
00080   int nby = pm->nby;
00081   int nbz = pm->nbz;
00082   int numatoms = pm->numatoms;
00083   int n, ib, jb, kb, index;
00084 
00085   /* reset all bin and next IDs */
00086   for (n = 0;  n < numbins;  n++)  bin[n] = -1;
00087   for (n = 0;  n < numatoms;  n++)  next[n] = -1;
00088 
00089   for (n = 0;  n < numatoms;  n++, atom += 4) {
00090     /* transform coordinates to reciprocal space */
00091     d[X] = atom[X] - c[X];
00092     d[Y] = atom[Y] - c[Y];
00093     d[Z] = atom[Z] - c[Z];
00094     s[X] = ru[X] * d[X] + ru[Y] * d[Y] + ru[Z] * d[Z];
00095     s[Y] = rv[X] * d[X] + rv[Y] * d[Y] + rv[Z] * d[Z];
00096     s[Z] = rw[X] * d[X] + rw[Y] * d[Y] + rw[Z] * d[Z];
00097     /* determine bin indexing in 3D */
00098     ib = (int) floor((s[X] + 0.5) * nbx);
00099     jb = (int) floor((s[Y] + 0.5) * nby);
00100     kb = (int) floor((s[Z] + 0.5) * nbz);
00101     /* assume coordinate is within defined cell;
00102      * adjust bin index in case of roundoff error */
00103     if      (ib < 0)    ib = 0;
00104     else if (ib >= nbx) ib = nbx - 1;
00105     if      (jb < 0)    jb = 0;
00106     else if (jb >= nby) jb = nby - 1;
00107     if      (kb < 0)    kb = 0;
00108     else if (kb >= nbz) kb = nbz - 1;
00109     index = (kb*nby + jb)*nbx + ib;  /* flatten 3D indexing into 1D */
00110     next[n] = bin[index];  /* attach atom index to front of linked list */
00111     bin[index] = n;
00112   }
00113   return NL_MSM_SUCCESS;
00114 }
00115 
00116 
00117 int bin_evaluation_1away(NL_Msm *pm) {
00118   enum { NUM_NBRBINS = 14 };
00119   int nbrbin[3*NUM_NBRBINS] = {
00120     0,0,0,  1,0,0,  -1,1,0,  0,1,0,  1,1,0,  -1,-1,1,  0,-1,1,
00121     1,-1,1,  -1,0,1,  0,0,1,  1,0,1,  -1,1,1,  0,1,1,  1,1,1
00122   };
00123   int nbx = pm->nbx;
00124   int nby = pm->nby;
00125   int nbz = pm->nbz;
00126   int aindex, bindex;
00127   int ia, ja, ka, n, i, j;
00128   int ispx = ((pm->msmflags & NL_MSM_PERIODIC_VEC1) != 0);
00129   int ispy = ((pm->msmflags & NL_MSM_PERIODIC_VEC2) != 0);
00130   int ispz = ((pm->msmflags & NL_MSM_PERIODIC_VEC3) != 0);
00131   int split = pm->split;
00132   const double *u = pm->cellvec1;
00133   const double *v = pm->cellvec2;
00134   const double *w = pm->cellvec3;
00135   const int *bin = pm->bin;
00136   const int *next = pm->next;
00137   const double *atom = pm->atom;
00138   double *force = pm->felec;
00139   double a2 = pm->a * pm->a;  /* cutoff^2 */
00140   double a_1 = 1 / pm->a;     /* 1 / cutoff */
00141   double a_2 = a_1 * a_1;     /* 1 / cutoff^2 */
00142   double u_elec = 0;  /* accumulate potential energy from electrostatics */
00143 
00144   for (aindex = 0, ka = 0;  ka < nbz;  ka++) {
00145     for (ja = 0;  ja < nby;  ja++) {
00146       for (ia = 0;  ia < nbx;  ia++, aindex++) {  /* loop over bins A */
00147 
00148         for (n = 0;  n < NUM_NBRBINS;  n++) { /* loop B-bin neighborhood of A */
00149 
00150           double p[3] = { 0,0,0 };  /* periodic wrapping vector for B bin */
00151 
00152           int ib = ia + nbrbin[3*n + X];  /* index for B bin */
00153           int jb = ja + nbrbin[3*n + Y];  /* index for B bin */
00154           int kb = ka + nbrbin[3*n + Z];  /* index for B bin */
00155 
00156           /* do wrap around for bin index outside of periodic dimension range,
00157            * short-circuit loop for bin index outside non-periodic range */
00158           if (ispx) {
00159             if (ib < 0) {
00160               ib += nbx;  p[X] -= u[X];  p[Y] -= u[Y];  p[Z] -= u[Z];
00161             }
00162             else if (ib >= nbx) {
00163               ib -= nbx;  p[X] += u[X];  p[Y] += u[Y];  p[Z] += u[Z];
00164             }
00165           }
00166           else if (ib < 0 || ib >= nbx) continue;
00167 
00168           if (ispy) {
00169             if (jb < 0) {
00170               jb += nby;  p[X] -= v[X];  p[Y] -= v[Y];  p[Z] -= v[Z];
00171             }
00172             else if (jb >= nby) {
00173               jb -= nby;  p[X] += v[X];  p[Y] += v[Y];  p[Z] += v[Z];
00174             }
00175           }
00176           else if (jb < 0 || jb >= nby) continue;
00177 
00178           if (ispz) {
00179             if (kb < 0) {
00180               kb += nbz;  p[X] -= w[X];  p[Y] -= w[Y];  p[Z] -= w[Z];
00181             }
00182             else if (kb >= nbz) {
00183               kb -= nbz;  p[X] += w[X];  p[Y] += w[Y];  p[Z] += w[Z];
00184             }
00185           }
00186           else if (kb < 0 || kb >= nbz) continue;
00187 
00188           /* flat 1D index for B bin, after doing wrap around */
00189           bindex = (kb*nby + jb)*nbx + ib;
00190 
00191           for (j = bin[bindex];  j != -1;  j = next[j]) { /* loop over B bin */
00192             double force_j[3] = { 0,0,0 };
00193             double atom_j[3];
00194             double qj = atom[4*j + Q];
00195             int ihead;
00196 
00197             atom_j[X] = atom[4*j + X] + p[X];
00198             atom_j[Y] = atom[4*j + Y] + p[Y];
00199             atom_j[Z] = atom[4*j + Z] + p[Z];
00200 
00201             ihead = bin[aindex];
00202             /* for self bin (A==B) interactions, visit each pair only once */
00203             if (n == 0) ihead = next[j];  /* i.e. aindex == bindex */
00204 
00205             for (i = ihead;  i != -1;  i = next[i]) { /* loop over A bin */
00206               double rij[3];
00207               double r2;
00208 
00209               rij[X] = atom_j[X] - atom[4*i + X];
00210               rij[Y] = atom_j[Y] - atom[4*i + Y];
00211               rij[Z] = atom_j[Z] - atom[4*i + Z];
00212 
00213               r2 = rij[X] * rij[X] + rij[Y] * rij[Y] + rij[Z] * rij[Z];
00214 
00215               if (r2 < a2) {
00216                 double fij[3];
00217                 double qq = qj * atom[4*i + Q];  /* combined charge */
00218                 double r;      /* length of vector r_ij */
00219                 double r_1;    /* 1/r */
00220                 double r_2;    /* 1/r^2 */
00221                 double r_a;    /* r/a */
00222                 double g;      /* normalized smoothing g(R), R=r/a */
00223                 double dg;     /* (d/dR)g(R) */
00224                 double ue;     /* U_elec(r) */
00225                 double due_r;  /* (1/r)*(d/dr)U_elec(r) */
00226 
00227                 r = sqrt(r2);
00228                 r_1 = 1/r;
00229                 r_2 = r_1 * r_1;
00230 
00231                 /* calculate MSM splitting */
00232                 r_a = r * a_1;
00233                 SPOLY(&g, &dg, r_a, split);
00234 
00235                 ue = qq * (r_1 - a_1 * g);
00236                 due_r = qq * r_1 * (-r_2 - a_2 * dg);
00237 
00238                 fij[X] = -rij[X] * due_r;
00239                 fij[Y] = -rij[Y] * due_r;
00240                 fij[Z] = -rij[Z] * due_r;
00241                 force[3*i + X] -= fij[X];
00242                 force[3*i + Y] -= fij[Y];
00243                 force[3*i + Z] -= fij[Z];
00244                 force_j[X] += fij[X];
00245                 force_j[Y] += fij[Y];
00246                 force_j[Z] += fij[Z];
00247 
00248                 u_elec += ue;
00249 
00250                 /* XXX virial? */
00251 
00252               } /* end if r2 < cutoff2 */
00253 
00254             } /* end loop over A bin */
00255 
00256             force[3*j + X] += force_j[X];
00257             force[3*j + Y] += force_j[Y];
00258             force[3*j + Z] += force_j[Z];
00259           } /* end loop over B bin */
00260 
00261         } /* end loop B-bin neighborhood of A */
00262 
00263       } /* end loop over bins A */
00264     }
00265   }
00266 
00267   pm->uelec += u_elec;
00268 
00269   return NL_MSM_SUCCESS;
00270 }
00271 
00272 
00273 int bin_evaluation_k_away(NL_Msm *pm) {
00274   int *nbrhd = pm->nbrhd;
00275   int nbrhdlen = pm->nbrhdlen;
00276   int nbx = pm->nbx;
00277   int nby = pm->nby;
00278   int nbz = pm->nbz;
00279   int aindex, bindex;
00280   int ia, ja, ka, n, i, j;
00281   int ispx = ((pm->msmflags & NL_MSM_PERIODIC_VEC1) != 0);
00282   int ispy = ((pm->msmflags & NL_MSM_PERIODIC_VEC2) != 0);
00283   int ispz = ((pm->msmflags & NL_MSM_PERIODIC_VEC3) != 0);
00284   int split = pm->split;
00285   const double *u = pm->cellvec1;
00286   const double *v = pm->cellvec2;
00287   const double *w = pm->cellvec3;
00288   const int *bin = pm->bin;
00289   const int *next = pm->next;
00290   const double *atom = pm->atom;
00291   double *force = pm->felec;
00292   double a2 = pm->a * pm->a;  /* cutoff^2 */
00293   double a_1 = 1 / pm->a;     /* 1 / cutoff */
00294   double a_2 = a_1 * a_1;     /* 1 / cutoff^2 */
00295   double u_elec = 0;  /* accumulate potential energy from electrostatics */
00296 
00297   for (aindex = 0, ka = 0;  ka < nbz;  ka++) {
00298     for (ja = 0;  ja < nby;  ja++) {
00299       for (ia = 0;  ia < nbx;  ia++, aindex++) {  /* loop over bins A */
00300 
00301         for (n = 0;  n < nbrhdlen;  n++) { /* loop B-bin neighborhood of A */
00302 
00303           double p[3] = { 0,0,0 };  /* periodic wrapping vector for B bin */
00304 
00305           int ib = ia + nbrhd[3*n + X];  /* index for B bin */
00306           int jb = ja + nbrhd[3*n + Y];  /* index for B bin */
00307           int kb = ka + nbrhd[3*n + Z];  /* index for B bin */
00308 
00309           /* do wrap around for bin index outside of periodic dimension range,
00310            * short-circuit loop for bin index outside non-periodic range */
00311           if (ispx) {
00312             if (ib < 0) {
00313               ib += nbx;  p[X] -= u[X];  p[Y] -= u[Y];  p[Z] -= u[Z];
00314             }
00315             else if (ib >= nbx) {
00316               ib -= nbx;  p[X] += u[X];  p[Y] += u[Y];  p[Z] += u[Z];
00317             }
00318           }
00319           else if (ib < 0 || ib >= nbx) continue;
00320 
00321           if (ispy) {
00322             if (jb < 0) {
00323               jb += nby;  p[X] -= v[X];  p[Y] -= v[Y];  p[Z] -= v[Z];
00324             }
00325             else if (jb >= nby) {
00326               jb -= nby;  p[X] += v[X];  p[Y] += v[Y];  p[Z] += v[Z];
00327             }
00328           }
00329           else if (jb < 0 || jb >= nby) continue;
00330 
00331           if (ispz) {
00332             if (kb < 0) {
00333               kb += nbz;  p[X] -= w[X];  p[Y] -= w[Y];  p[Z] -= w[Z];
00334             }
00335             else if (kb >= nbz) {
00336               kb -= nbz;  p[X] += w[X];  p[Y] += w[Y];  p[Z] += w[Z];
00337             }
00338           }
00339           else if (kb < 0 || kb >= nbz) continue;
00340 
00341           /* flat 1D index for B bin, after doing wrap around */
00342           bindex = (kb*nby + jb)*nbx + ib;
00343 
00344           for (j = bin[bindex];  j != -1;  j = next[j]) { /* loop over B bin */
00345             double force_j[3] = { 0,0,0 };
00346             double atom_j[3];
00347             double qj = atom[4*j + Q];
00348             int ihead;
00349 
00350             atom_j[X] = atom[4*j + X] + p[X];
00351             atom_j[Y] = atom[4*j + Y] + p[Y];
00352             atom_j[Z] = atom[4*j + Z] + p[Z];
00353 
00354             ihead = bin[aindex];
00355             /* for self bin (A==B) interactions, visit each pair only once */
00356             if (n == 0) ihead = next[j];  /* i.e. aindex == bindex */
00357 
00358             for (i = ihead;  i != -1;  i = next[i]) { /* loop over A bin */
00359               double rij[3];
00360               double r2;
00361 
00362               rij[X] = atom_j[X] - atom[4*i + X];
00363               rij[Y] = atom_j[Y] - atom[4*i + Y];
00364               rij[Z] = atom_j[Z] - atom[4*i + Z];
00365 
00366               r2 = rij[X] * rij[X] + rij[Y] * rij[Y] + rij[Z] * rij[Z];
00367 
00368               if (r2 < a2) {
00369                 double fij[3];
00370                 double qq = qj * atom[4*i + Q];  /* combined charge */
00371                 double r;      /* length of vector r_ij */
00372                 double r_1;    /* 1/r */
00373                 double r_2;    /* 1/r^2 */
00374                 double r_a;    /* r/a */
00375                 double g;      /* normalized smoothing g(R), R=r/a */
00376                 double dg;     /* (d/dR)g(R) */
00377                 double ue;     /* U_elec(r) */
00378                 double due_r;  /* (1/r)*(d/dr)U_elec(r) */
00379 
00380                 r = sqrt(r2);
00381                 r_1 = 1/r;
00382                 r_2 = r_1 * r_1;
00383 
00384                 /* calculate MSM splitting */
00385                 r_a = r * a_1;
00386                 SPOLY(&g, &dg, r_a, split);
00387 
00388                 ue = qq * (r_1 - a_1 * g);
00389                 due_r = qq * r_1 * (-r_2 - a_2 * dg);
00390 
00391                 fij[X] = -rij[X] * due_r;
00392                 fij[Y] = -rij[Y] * due_r;
00393                 fij[Z] = -rij[Z] * due_r;
00394                 force[3*i + X] -= fij[X];
00395                 force[3*i + Y] -= fij[Y];
00396                 force[3*i + Z] -= fij[Z];
00397                 force_j[X] += fij[X];
00398                 force_j[Y] += fij[Y];
00399                 force_j[Z] += fij[Z];
00400 
00401                 u_elec += ue;
00402 
00403                 /* XXX virial? */
00404 
00405               } /* end if r2 < cutoff2 */
00406 
00407             } /* end loop over A bin */
00408 
00409             force[3*j + X] += force_j[X];
00410             force[3*j + Y] += force_j[Y];
00411             force[3*j + Z] += force_j[Z];
00412           } /* end loop over B bin */
00413 
00414         } /* end loop B-bin neighborhood of A */
00415 
00416       } /* end loop over bins A */
00417     }
00418   }
00419 
00420   pm->uelec += u_elec;
00421 
00422   return NL_MSM_SUCCESS;
00423 }