msm_shortrng.c File Reference

#include "msm_defn.h"

Go to the source code of this file.

Functions

static int setup_bin_data (NL_Msm *pm)
static int spatial_hashing (NL_Msm *pm)
static int bin_evaluation_1away (NL_Msm *pm)
static int bin_evaluation_k_away (NL_Msm *pm)
int NL_msm_compute_short_range (NL_Msm *pm)


Function Documentation

int bin_evaluation_1away ( NL_Msm pm  )  [static]

Definition at line 117 of file msm_shortrng.c.

References NL_Msm_t::a, NL_Msm_t::atom, NL_Msm_t::bin, NL_Msm_t::cellvec1, NL_Msm_t::cellvec2, NL_Msm_t::cellvec3, NL_Msm_t::felec, j, NL_Msm_t::msmflags, NL_Msm_t::nbx, NL_Msm_t::nby, NL_Msm_t::nbz, NL_Msm_t::next, next(), NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, Q, NL_Msm_t::split, eabffunc::split(), SPOLY, NL_Msm_t::uelec, X, Y, and Z.

Referenced by NL_msm_compute_short_range(), and NL_msm_compute_short_range_sprec().

00117                                      {
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 }

int bin_evaluation_k_away ( NL_Msm pm  )  [static]

Definition at line 273 of file msm_shortrng.c.

References NL_Msm_t::a, NL_Msm_t::atom, NL_Msm_t::bin, NL_Msm_t::cellvec1, NL_Msm_t::cellvec2, NL_Msm_t::cellvec3, NL_Msm_t::felec, NL_Msm_t::msmflags, NL_Msm_t::nbrhd, NL_Msm_t::nbrhdlen, NL_Msm_t::nbx, NL_Msm_t::nby, NL_Msm_t::nbz, NL_Msm_t::next, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, Q, NL_Msm_t::split, SPOLY, NL_Msm_t::uelec, X, Y, and Z.

Referenced by NL_msm_compute_short_range(), and NL_msm_compute_short_range_sprec().

00273                                       {
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 }

int NL_msm_compute_short_range ( NL_Msm pm  ) 

Definition at line 37 of file msm_shortrng.c.

References bin_evaluation_1away(), bin_evaluation_k_away(), NL_Msm_t::maxatoms, NL_Msm_t::msmflags, NL_MSM_COMPUTE_1AWAY, NL_MSM_SUCCESS, NL_Msm_t::numatoms, setup_bin_data(), and spatial_hashing().

Referenced by NL_msm_compute_force().

00037                                            {
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 }

int setup_bin_data ( NL_Msm pm  )  [static]

msm_shortrng.c

Compute the short-range part of MSM forces.

Perform spatial hashing of atoms into bins. The implementation uses the simplest approach where the bin lengths are at least as large as the cutoff distance, so only nearest neighbor bins need be checked. The bins are represented by a simple linked list of atom indices using fixed array storage: an array of int for the bins giving the index of a "head" atom in the bin and an array of int for the "next" atom in the bin, where -1 denotes the end of the list.

The formulation of bins is determined for an arbitrary set of basis vectors and the spatial hashing is performed in reciprocal space, so this part of the MSM code is ready to support non-orthogonal basis vectors.

XXX The virial is not calculated.

XXX The idea of excluded interaction pairs is totally disregarded. This should be fine for calculating normal systems of atoms, where the full 1/r interaction is later subtracted out for excluded pairs without causing excessive numerical error. However, doing so will be a problem for systems containing Drude particles, possibly also for lone pairs, and is unsuitable for calculating Lennard-Jones.

Definition at line 57 of file msm_shortrng.c.

References ASSERT, NL_Msm_t::maxatoms, NL_Msm_t::next, NL_MSM_ERROR_MALLOC, NL_MSM_SUCCESS, and NL_Msm_t::numatoms.

Referenced by NL_msm_compute_short_range(), and NL_msm_compute_short_range_sprec().

00057                                {
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 }

int spatial_hashing ( NL_Msm pm  )  [static]

Definition at line 69 of file msm_shortrng.c.

References NL_Msm_t::atom, NL_Msm_t::bin, NL_Msm_t::cellcenter, NL_Msm_t::nbx, NL_Msm_t::nby, NL_Msm_t::nbz, NL_Msm_t::next, next(), NL_MSM_SUCCESS, NL_Msm_t::numatoms, numatoms, NL_Msm_t::numbins, NL_Msm_t::recipvec1, NL_Msm_t::recipvec2, NL_Msm_t::recipvec3, X, Y, and Z.

Referenced by NL_msm_compute_short_range(), and NL_msm_compute_short_range_sprec().

00069                                 {
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 }


Generated on Tue Nov 21 01:17:16 2017 for NAMD by  doxygen 1.4.7