NAMD
Functions
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, NL_Msm_t::msmflags, NL_Msm_t::nbx, NL_Msm_t::nby, NL_Msm_t::nbz, next(), NL_Msm_t::next, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, Q, split(), NL_Msm_t::split, SPOLY, NL_Msm_t::uelec, X, Y, and Z.

Referenced by NL_msm_compute_short_range().

117  {
118  enum { NUM_NBRBINS = 14 };
119  int nbrbin[3*NUM_NBRBINS] = {
120  0,0,0, 1,0,0, -1,1,0, 0,1,0, 1,1,0, -1,-1,1, 0,-1,1,
121  1,-1,1, -1,0,1, 0,0,1, 1,0,1, -1,1,1, 0,1,1, 1,1,1
122  };
123  int nbx = pm->nbx;
124  int nby = pm->nby;
125  int nbz = pm->nbz;
126  int aindex, bindex;
127  int ia, ja, ka, n, i, j;
128  int ispx = ((pm->msmflags & NL_MSM_PERIODIC_VEC1) != 0);
129  int ispy = ((pm->msmflags & NL_MSM_PERIODIC_VEC2) != 0);
130  int ispz = ((pm->msmflags & NL_MSM_PERIODIC_VEC3) != 0);
131  int split = pm->split;
132  const double *u = pm->cellvec1;
133  const double *v = pm->cellvec2;
134  const double *w = pm->cellvec3;
135  const int *bin = pm->bin;
136  const int *next = pm->next;
137  const double *atom = pm->atom;
138  double *force = pm->felec;
139  double a2 = pm->a * pm->a; /* cutoff^2 */
140  double a_1 = 1 / pm->a; /* 1 / cutoff */
141  double a_2 = a_1 * a_1; /* 1 / cutoff^2 */
142  double u_elec = 0; /* accumulate potential energy from electrostatics */
143 
144  for (aindex = 0, ka = 0; ka < nbz; ka++) {
145  for (ja = 0; ja < nby; ja++) {
146  for (ia = 0; ia < nbx; ia++, aindex++) { /* loop over bins A */
147 
148  for (n = 0; n < NUM_NBRBINS; n++) { /* loop B-bin neighborhood of A */
149 
150  double p[3] = { 0,0,0 }; /* periodic wrapping vector for B bin */
151 
152  int ib = ia + nbrbin[3*n + X]; /* index for B bin */
153  int jb = ja + nbrbin[3*n + Y]; /* index for B bin */
154  int kb = ka + nbrbin[3*n + Z]; /* index for B bin */
155 
156  /* do wrap around for bin index outside of periodic dimension range,
157  * short-circuit loop for bin index outside non-periodic range */
158  if (ispx) {
159  if (ib < 0) {
160  ib += nbx; p[X] -= u[X]; p[Y] -= u[Y]; p[Z] -= u[Z];
161  }
162  else if (ib >= nbx) {
163  ib -= nbx; p[X] += u[X]; p[Y] += u[Y]; p[Z] += u[Z];
164  }
165  }
166  else if (ib < 0 || ib >= nbx) continue;
167 
168  if (ispy) {
169  if (jb < 0) {
170  jb += nby; p[X] -= v[X]; p[Y] -= v[Y]; p[Z] -= v[Z];
171  }
172  else if (jb >= nby) {
173  jb -= nby; p[X] += v[X]; p[Y] += v[Y]; p[Z] += v[Z];
174  }
175  }
176  else if (jb < 0 || jb >= nby) continue;
177 
178  if (ispz) {
179  if (kb < 0) {
180  kb += nbz; p[X] -= w[X]; p[Y] -= w[Y]; p[Z] -= w[Z];
181  }
182  else if (kb >= nbz) {
183  kb -= nbz; p[X] += w[X]; p[Y] += w[Y]; p[Z] += w[Z];
184  }
185  }
186  else if (kb < 0 || kb >= nbz) continue;
187 
188  /* flat 1D index for B bin, after doing wrap around */
189  bindex = (kb*nby + jb)*nbx + ib;
190 
191  for (j = bin[bindex]; j != -1; j = next[j]) { /* loop over B bin */
192  double force_j[3] = { 0,0,0 };
193  double atom_j[3];
194  double qj = atom[4*j + Q];
195  int ihead;
196 
197  atom_j[X] = atom[4*j + X] + p[X];
198  atom_j[Y] = atom[4*j + Y] + p[Y];
199  atom_j[Z] = atom[4*j + Z] + p[Z];
200 
201  ihead = bin[aindex];
202  /* for self bin (A==B) interactions, visit each pair only once */
203  if (n == 0) ihead = next[j]; /* i.e. aindex == bindex */
204 
205  for (i = ihead; i != -1; i = next[i]) { /* loop over A bin */
206  double rij[3];
207  double r2;
208 
209  rij[X] = atom_j[X] - atom[4*i + X];
210  rij[Y] = atom_j[Y] - atom[4*i + Y];
211  rij[Z] = atom_j[Z] - atom[4*i + Z];
212 
213  r2 = rij[X] * rij[X] + rij[Y] * rij[Y] + rij[Z] * rij[Z];
214 
215  if (r2 < a2) {
216  double fij[3];
217  double qq = qj * atom[4*i + Q]; /* combined charge */
218  double r; /* length of vector r_ij */
219  double r_1; /* 1/r */
220  double r_2; /* 1/r^2 */
221  double r_a; /* r/a */
222  double g; /* normalized smoothing g(R), R=r/a */
223  double dg; /* (d/dR)g(R) */
224  double ue; /* U_elec(r) */
225  double due_r; /* (1/r)*(d/dr)U_elec(r) */
226 
227  r = sqrt(r2);
228  r_1 = 1/r;
229  r_2 = r_1 * r_1;
230 
231  /* calculate MSM splitting */
232  r_a = r * a_1;
233  SPOLY(&g, &dg, r_a, split);
234 
235  ue = qq * (r_1 - a_1 * g);
236  due_r = qq * r_1 * (-r_2 - a_2 * dg);
237 
238  fij[X] = -rij[X] * due_r;
239  fij[Y] = -rij[Y] * due_r;
240  fij[Z] = -rij[Z] * due_r;
241  force[3*i + X] -= fij[X];
242  force[3*i + Y] -= fij[Y];
243  force[3*i + Z] -= fij[Z];
244  force_j[X] += fij[X];
245  force_j[Y] += fij[Y];
246  force_j[Z] += fij[Z];
247 
248  u_elec += ue;
249 
250  /* XXX virial? */
251 
252  } /* end if r2 < cutoff2 */
253 
254  } /* end loop over A bin */
255 
256  force[3*j + X] += force_j[X];
257  force[3*j + Y] += force_j[Y];
258  force[3*j + Z] += force_j[Z];
259  } /* end loop over B bin */
260 
261  } /* end loop B-bin neighborhood of A */
262 
263  } /* end loop over bins A */
264  }
265  }
266 
267  pm->uelec += u_elec;
268 
269  return NL_MSM_SUCCESS;
270 }
int msmflags
Definition: msm_defn.h:590
double * felec
Definition: msm_defn.h:568
double cellvec2[3]
Definition: msm_defn.h:575
int * bin
Definition: msm_defn.h:597
#define X
Definition: msm_defn.h:29
double cellvec1[3]
Definition: msm_defn.h:574
double uelec
Definition: msm_defn.h:566
double cellvec3[3]
Definition: msm_defn.h:576
static Units next(Units u)
Definition: ParseOptions.C:48
#define Z
Definition: msm_defn.h:33
double a
Definition: msm_defn.h:635
const double * atom
Definition: msm_defn.h:569
std::vector< std::string > split(const std::string &text, std::string delimiter)
Definition: MoleculeQM.C:73
#define Q
Definition: msm_defn.h:35
int nbx
Definition: msm_defn.h:594
int split
Definition: msm_defn.h:643
int nby
Definition: msm_defn.h:595
#define SPOLY(pg, pdg, ra, split)
Definition: msm_defn.h:140
#define Y
Definition: msm_defn.h:31
int nbz
Definition: msm_defn.h:596
int * next
Definition: msm_defn.h:601
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, next(), NL_Msm_t::next, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, Q, split(), NL_Msm_t::split, SPOLY, NL_Msm_t::uelec, X, Y, and Z.

Referenced by NL_msm_compute_short_range().

273  {
274  int *nbrhd = pm->nbrhd;
275  int nbrhdlen = pm->nbrhdlen;
276  int nbx = pm->nbx;
277  int nby = pm->nby;
278  int nbz = pm->nbz;
279  int aindex, bindex;
280  int ia, ja, ka, n, i, j;
281  int ispx = ((pm->msmflags & NL_MSM_PERIODIC_VEC1) != 0);
282  int ispy = ((pm->msmflags & NL_MSM_PERIODIC_VEC2) != 0);
283  int ispz = ((pm->msmflags & NL_MSM_PERIODIC_VEC3) != 0);
284  int split = pm->split;
285  const double *u = pm->cellvec1;
286  const double *v = pm->cellvec2;
287  const double *w = pm->cellvec3;
288  const int *bin = pm->bin;
289  const int *next = pm->next;
290  const double *atom = pm->atom;
291  double *force = pm->felec;
292  double a2 = pm->a * pm->a; /* cutoff^2 */
293  double a_1 = 1 / pm->a; /* 1 / cutoff */
294  double a_2 = a_1 * a_1; /* 1 / cutoff^2 */
295  double u_elec = 0; /* accumulate potential energy from electrostatics */
296 
297  for (aindex = 0, ka = 0; ka < nbz; ka++) {
298  for (ja = 0; ja < nby; ja++) {
299  for (ia = 0; ia < nbx; ia++, aindex++) { /* loop over bins A */
300 
301  for (n = 0; n < nbrhdlen; n++) { /* loop B-bin neighborhood of A */
302 
303  double p[3] = { 0,0,0 }; /* periodic wrapping vector for B bin */
304 
305  int ib = ia + nbrhd[3*n + X]; /* index for B bin */
306  int jb = ja + nbrhd[3*n + Y]; /* index for B bin */
307  int kb = ka + nbrhd[3*n + Z]; /* index for B bin */
308 
309  /* do wrap around for bin index outside of periodic dimension range,
310  * short-circuit loop for bin index outside non-periodic range */
311  if (ispx) {
312  if (ib < 0) {
313  ib += nbx; p[X] -= u[X]; p[Y] -= u[Y]; p[Z] -= u[Z];
314  }
315  else if (ib >= nbx) {
316  ib -= nbx; p[X] += u[X]; p[Y] += u[Y]; p[Z] += u[Z];
317  }
318  }
319  else if (ib < 0 || ib >= nbx) continue;
320 
321  if (ispy) {
322  if (jb < 0) {
323  jb += nby; p[X] -= v[X]; p[Y] -= v[Y]; p[Z] -= v[Z];
324  }
325  else if (jb >= nby) {
326  jb -= nby; p[X] += v[X]; p[Y] += v[Y]; p[Z] += v[Z];
327  }
328  }
329  else if (jb < 0 || jb >= nby) continue;
330 
331  if (ispz) {
332  if (kb < 0) {
333  kb += nbz; p[X] -= w[X]; p[Y] -= w[Y]; p[Z] -= w[Z];
334  }
335  else if (kb >= nbz) {
336  kb -= nbz; p[X] += w[X]; p[Y] += w[Y]; p[Z] += w[Z];
337  }
338  }
339  else if (kb < 0 || kb >= nbz) continue;
340 
341  /* flat 1D index for B bin, after doing wrap around */
342  bindex = (kb*nby + jb)*nbx + ib;
343 
344  for (j = bin[bindex]; j != -1; j = next[j]) { /* loop over B bin */
345  double force_j[3] = { 0,0,0 };
346  double atom_j[3];
347  double qj = atom[4*j + Q];
348  int ihead;
349 
350  atom_j[X] = atom[4*j + X] + p[X];
351  atom_j[Y] = atom[4*j + Y] + p[Y];
352  atom_j[Z] = atom[4*j + Z] + p[Z];
353 
354  ihead = bin[aindex];
355  /* for self bin (A==B) interactions, visit each pair only once */
356  if (n == 0) ihead = next[j]; /* i.e. aindex == bindex */
357 
358  for (i = ihead; i != -1; i = next[i]) { /* loop over A bin */
359  double rij[3];
360  double r2;
361 
362  rij[X] = atom_j[X] - atom[4*i + X];
363  rij[Y] = atom_j[Y] - atom[4*i + Y];
364  rij[Z] = atom_j[Z] - atom[4*i + Z];
365 
366  r2 = rij[X] * rij[X] + rij[Y] * rij[Y] + rij[Z] * rij[Z];
367 
368  if (r2 < a2) {
369  double fij[3];
370  double qq = qj * atom[4*i + Q]; /* combined charge */
371  double r; /* length of vector r_ij */
372  double r_1; /* 1/r */
373  double r_2; /* 1/r^2 */
374  double r_a; /* r/a */
375  double g; /* normalized smoothing g(R), R=r/a */
376  double dg; /* (d/dR)g(R) */
377  double ue; /* U_elec(r) */
378  double due_r; /* (1/r)*(d/dr)U_elec(r) */
379 
380  r = sqrt(r2);
381  r_1 = 1/r;
382  r_2 = r_1 * r_1;
383 
384  /* calculate MSM splitting */
385  r_a = r * a_1;
386  SPOLY(&g, &dg, r_a, split);
387 
388  ue = qq * (r_1 - a_1 * g);
389  due_r = qq * r_1 * (-r_2 - a_2 * dg);
390 
391  fij[X] = -rij[X] * due_r;
392  fij[Y] = -rij[Y] * due_r;
393  fij[Z] = -rij[Z] * due_r;
394  force[3*i + X] -= fij[X];
395  force[3*i + Y] -= fij[Y];
396  force[3*i + Z] -= fij[Z];
397  force_j[X] += fij[X];
398  force_j[Y] += fij[Y];
399  force_j[Z] += fij[Z];
400 
401  u_elec += ue;
402 
403  /* XXX virial? */
404 
405  } /* end if r2 < cutoff2 */
406 
407  } /* end loop over A bin */
408 
409  force[3*j + X] += force_j[X];
410  force[3*j + Y] += force_j[Y];
411  force[3*j + Z] += force_j[Z];
412  } /* end loop over B bin */
413 
414  } /* end loop B-bin neighborhood of A */
415 
416  } /* end loop over bins A */
417  }
418  }
419 
420  pm->uelec += u_elec;
421 
422  return NL_MSM_SUCCESS;
423 }
int msmflags
Definition: msm_defn.h:590
double * felec
Definition: msm_defn.h:568
double cellvec2[3]
Definition: msm_defn.h:575
int * bin
Definition: msm_defn.h:597
#define X
Definition: msm_defn.h:29
double cellvec1[3]
Definition: msm_defn.h:574
double uelec
Definition: msm_defn.h:566
double cellvec3[3]
Definition: msm_defn.h:576
static Units next(Units u)
Definition: ParseOptions.C:48
#define Z
Definition: msm_defn.h:33
double a
Definition: msm_defn.h:635
const double * atom
Definition: msm_defn.h:569
std::vector< std::string > split(const std::string &text, std::string delimiter)
Definition: MoleculeQM.C:73
#define Q
Definition: msm_defn.h:35
int nbx
Definition: msm_defn.h:594
int split
Definition: msm_defn.h:643
int * nbrhd
Definition: msm_defn.h:612
int nby
Definition: msm_defn.h:595
#define SPOLY(pg, pdg, ra, split)
Definition: msm_defn.h:140
#define Y
Definition: msm_defn.h:31
int nbz
Definition: msm_defn.h:596
int nbrhdlen
Definition: msm_defn.h:613
int * next
Definition: msm_defn.h:601
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().

37  {
38  int rc = 0; /* return code */
39 
40  if (pm->maxatoms < pm->numatoms) {
41  rc = setup_bin_data(pm);
42  if (rc) return rc;
43  }
44  rc = spatial_hashing(pm);
45  if (rc) return rc;
46  if (pm->msmflags & NL_MSM_COMPUTE_1AWAY) {
47  rc = bin_evaluation_1away(pm);
48  }
49  else {
50  rc = bin_evaluation_k_away(pm);
51  }
52  if (rc) return rc;
53  return NL_MSM_SUCCESS;
54 }
int msmflags
Definition: msm_defn.h:590
int numatoms
Definition: msm_defn.h:599
static int spatial_hashing(NL_Msm *pm)
Definition: msm_shortrng.c:69
static int setup_bin_data(NL_Msm *pm)
Definition: msm_shortrng.c:57
static int bin_evaluation_k_away(NL_Msm *pm)
Definition: msm_shortrng.c:273
static int bin_evaluation_1away(NL_Msm *pm)
Definition: msm_shortrng.c:117
int maxatoms
Definition: msm_defn.h:600
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().

57  {
58  void *vptr = NULL;
59 
60  ASSERT(pm->maxatoms < pm->numatoms);
61  vptr = malloc(pm->numatoms * sizeof(int));
62  if (vptr == NULL) return NL_MSM_ERROR_MALLOC;
63  pm->next = (int *) vptr;
64  pm->maxatoms = pm->numatoms;
65  return NL_MSM_SUCCESS;
66 }
int numatoms
Definition: msm_defn.h:599
#define ASSERT(E)
int maxatoms
Definition: msm_defn.h:600
int * next
Definition: msm_defn.h:601
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, next(), NL_Msm_t::next, NL_MSM_SUCCESS, numatoms, NL_Msm_t::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().

69  {
70  double d[3], s[3];
71  const double *atom = pm->atom; /* stored x/y/z/q */
72  const double *c = pm->cellcenter;
73  const double *ru = pm->recipvec1;
74  const double *rv = pm->recipvec2;
75  const double *rw = pm->recipvec3;
76  int *bin = pm->bin;
77  int *next = pm->next;
78  int numbins = pm->numbins;
79  int nbx = pm->nbx;
80  int nby = pm->nby;
81  int nbz = pm->nbz;
82  int numatoms = pm->numatoms;
83  int n, ib, jb, kb, index;
84 
85  /* reset all bin and next IDs */
86  for (n = 0; n < numbins; n++) bin[n] = -1;
87  for (n = 0; n < numatoms; n++) next[n] = -1;
88 
89  for (n = 0; n < numatoms; n++, atom += 4) {
90  /* transform coordinates to reciprocal space */
91  d[X] = atom[X] - c[X];
92  d[Y] = atom[Y] - c[Y];
93  d[Z] = atom[Z] - c[Z];
94  s[X] = ru[X] * d[X] + ru[Y] * d[Y] + ru[Z] * d[Z];
95  s[Y] = rv[X] * d[X] + rv[Y] * d[Y] + rv[Z] * d[Z];
96  s[Z] = rw[X] * d[X] + rw[Y] * d[Y] + rw[Z] * d[Z];
97  /* determine bin indexing in 3D */
98  ib = (int) floor((s[X] + 0.5) * nbx);
99  jb = (int) floor((s[Y] + 0.5) * nby);
100  kb = (int) floor((s[Z] + 0.5) * nbz);
101  /* assume coordinate is within defined cell;
102  * adjust bin index in case of roundoff error */
103  if (ib < 0) ib = 0;
104  else if (ib >= nbx) ib = nbx - 1;
105  if (jb < 0) jb = 0;
106  else if (jb >= nby) jb = nby - 1;
107  if (kb < 0) kb = 0;
108  else if (kb >= nbz) kb = nbz - 1;
109  index = (kb*nby + jb)*nbx + ib; /* flatten 3D indexing into 1D */
110  next[n] = bin[index]; /* attach atom index to front of linked list */
111  bin[index] = n;
112  }
113  return NL_MSM_SUCCESS;
114 }
double recipvec2[3]
Definition: msm_defn.h:579
double recipvec3[3]
Definition: msm_defn.h:580
int numbins
Definition: msm_defn.h:592
int * bin
Definition: msm_defn.h:597
int numatoms
Definition: msm_defn.h:599
static int numatoms
Definition: ScriptTcl.C:64
#define X
Definition: msm_defn.h:29
static Units next(Units u)
Definition: ParseOptions.C:48
double cellcenter[3]
Definition: msm_defn.h:577
#define Z
Definition: msm_defn.h:33
const double * atom
Definition: msm_defn.h:569
int nbx
Definition: msm_defn.h:594
int nby
Definition: msm_defn.h:595
#define Y
Definition: msm_defn.h:31
double recipvec1[3]
Definition: msm_defn.h:578
int nbz
Definition: msm_defn.h:596
int * next
Definition: msm_defn.h:601