NAMD
Functions
msm_shortrng_sprec.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_sprec (NL_Msm *pm)
 

Function Documentation

int bin_evaluation_1away ( NL_Msm pm)
static

Definition at line 117 of file msm_shortrng_sprec.c.

References NL_Msm_t::a_f, NL_Msm_t::atom_f, NL_Msm_t::bin, NL_Msm_t::cellvec1_f, NL_Msm_t::cellvec2_f, NL_Msm_t::cellvec3_f, NL_Msm_t::felec_f, 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_SPREC, NL_Msm_t::uelec, X, Y, and Z.

Referenced by NL_msm_compute_short_range_sprec().

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 float *u = pm->cellvec1_f;
133  const float *v = pm->cellvec2_f;
134  const float *w = pm->cellvec3_f;
135  const int *bin = pm->bin;
136  const int *next = pm->next;
137  const float *atom = pm->atom_f;
138  float *force = pm->felec_f;
139  float a2 = pm->a_f * pm->a_f; /* cutoff^2 */
140  float a_1 = 1 / pm->a_f; /* 1 / cutoff */
141  float a_2 = a_1 * a_1; /* 1 / cutoff^2 */
142  double u_elec = 0; /* need double precision for accumulating potential */
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  float 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  float force_j[3] = { 0,0,0 };
193  float atom_j[3];
194  float 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  float rij[3];
207  float 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  float fij[3];
217  float qq = qj * atom[4*i + Q]; /* combined charge */
218  float r; /* length of vector r_ij */
219  float r_1; /* 1/r */
220  float r_2; /* 1/r^2 */
221  float r_a; /* r/a */
222  float g; /* normalized smoothing g(R), R=r/a */
223  float dg; /* (d/dR)g(R) */
224  float ue; /* U_elec(r) */
225  float due_r; /* (1/r)*(d/dr)U_elec(r) */
226 
227  r = sqrtf(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_SPREC(&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
int * bin
Definition: msm_defn.h:597
const float * atom_f
Definition: msm_defn.h:572
float cellvec1_f[3]
Definition: msm_defn.h:582
#define X
Definition: msm_defn.h:29
#define SPOLY_SPREC(pg, pdg, ra, split)
Definition: msm_defn.h:351
double uelec
Definition: msm_defn.h:566
static Units next(Units u)
Definition: ParseOptions.C:48
float * felec_f
Definition: msm_defn.h:571
#define Z
Definition: msm_defn.h:33
float cellvec3_f[3]
Definition: msm_defn.h:584
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
float a_f
Definition: msm_defn.h:638
int nby
Definition: msm_defn.h:595
float cellvec2_f[3]
Definition: msm_defn.h:583
#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_sprec.c.

References NL_Msm_t::a_f, NL_Msm_t::atom_f, NL_Msm_t::bin, NL_Msm_t::cellvec1_f, NL_Msm_t::cellvec2_f, NL_Msm_t::cellvec3_f, NL_Msm_t::felec_f, 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_SPREC, NL_Msm_t::uelec, X, Y, and Z.

Referenced by NL_msm_compute_short_range_sprec().

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 float *u = pm->cellvec1_f;
286  const float *v = pm->cellvec2_f;
287  const float *w = pm->cellvec3_f;
288  const int *bin = pm->bin;
289  const int *next = pm->next;
290  const float *atom = pm->atom_f;
291  float *force = pm->felec_f;
292  float a2 = pm->a_f * pm->a_f; /* cutoff^2 */
293  float a_1 = 1 / pm->a_f; /* 1 / cutoff */
294  float a_2 = a_1 * a_1; /* 1 / cutoff^2 */
295  double u_elec = 0; /* need double precision for accumulating potential */
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  float 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  float force_j[3] = { 0,0,0 };
346  float atom_j[3];
347  float 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  float rij[3];
360  float 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  float fij[3];
370  float qq = qj * atom[4*i + Q]; /* combined charge */
371  float r; /* length of vector r_ij */
372  float r_1; /* 1/r */
373  float r_2; /* 1/r^2 */
374  float r_a; /* r/a */
375  float g; /* normalized smoothing g(R), R=r/a */
376  float dg; /* (d/dR)g(R) */
377  float ue; /* U_elec(r) */
378  float due_r; /* (1/r)*(d/dr)U_elec(r) */
379 
380  r = sqrtf(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_SPREC(&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
int * bin
Definition: msm_defn.h:597
const float * atom_f
Definition: msm_defn.h:572
float cellvec1_f[3]
Definition: msm_defn.h:582
#define X
Definition: msm_defn.h:29
#define SPOLY_SPREC(pg, pdg, ra, split)
Definition: msm_defn.h:351
double uelec
Definition: msm_defn.h:566
static Units next(Units u)
Definition: ParseOptions.C:48
float * felec_f
Definition: msm_defn.h:571
#define Z
Definition: msm_defn.h:33
float cellvec3_f[3]
Definition: msm_defn.h:584
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
float a_f
Definition: msm_defn.h:638
int nby
Definition: msm_defn.h:595
float cellvec2_f[3]
Definition: msm_defn.h:583
#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_sprec ( NL_Msm pm)

Definition at line 37 of file msm_shortrng_sprec.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_sprec().

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 bin_evaluation_1away(NL_Msm *pm)
static int setup_bin_data(NL_Msm *pm)
static int bin_evaluation_k_away(NL_Msm *pm)
int maxatoms
Definition: msm_defn.h:600
static int spatial_hashing(NL_Msm *pm)
int setup_bin_data ( NL_Msm pm)
static

msm_shortrng_sprec.c

Compute the short-range part of MSM forces, using single precision.

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_sprec.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_sprec().

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_sprec.c.

References NL_Msm_t::atom_f, NL_Msm_t::bin, NL_Msm_t::cellcenter_f, 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_f, NL_Msm_t::recipvec2_f, NL_Msm_t::recipvec3_f, X, Y, and Z.

Referenced by NL_msm_compute_short_range_sprec().

69  {
70  double d[3], s[3];
71  const float *atom = pm->atom_f; /* stored x/y/z/q */
72  const float *c = pm->cellcenter_f;
73  const float *ru = pm->recipvec1_f;
74  const float *rv = pm->recipvec2_f;
75  const float *rw = pm->recipvec3_f;
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) floorf((s[X] + 0.5f) * nbx);
99  jb = (int) floorf((s[Y] + 0.5f) * nby);
100  kb = (int) floorf((s[Z] + 0.5f) * 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 }
float cellcenter_f[3]
Definition: msm_defn.h:585
int numbins
Definition: msm_defn.h:592
int * bin
Definition: msm_defn.h:597
const float * atom_f
Definition: msm_defn.h:572
float recipvec2_f[3]
Definition: msm_defn.h:587
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
#define Z
Definition: msm_defn.h:33
float recipvec3_f[3]
Definition: msm_defn.h:588
int nbx
Definition: msm_defn.h:594
float recipvec1_f[3]
Definition: msm_defn.h:586
int nby
Definition: msm_defn.h:595
#define Y
Definition: msm_defn.h:31
int nbz
Definition: msm_defn.h:596
int * next
Definition: msm_defn.h:601