NAMD
msm_shortrng_sprec.c
Go to the documentation of this file.
1 
29 #include "msm_defn.h"
30 
31 static int setup_bin_data(NL_Msm *pm);
32 static int spatial_hashing(NL_Msm *pm);
33 static int bin_evaluation_1away(NL_Msm *pm);
34 static int bin_evaluation_k_away(NL_Msm *pm);
35 
36 
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 }
55 
56 
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 }
67 
68 
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 }
115 
116 
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 }
271 
272 
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 }
float cellcenter_f[3]
Definition: msm_defn.h:585
int msmflags
Definition: msm_defn.h:590
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
static int bin_evaluation_1away(NL_Msm *pm)
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
static int setup_bin_data(NL_Msm *pm)
double uelec
Definition: msm_defn.h:566
static Units next(Units u)
Definition: ParseOptions.C:48
int NL_msm_compute_short_range_sprec(NL_Msm *pm)
float * felec_f
Definition: msm_defn.h:571
#define Z
Definition: msm_defn.h:33
#define ASSERT(E)
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
float recipvec3_f[3]
Definition: msm_defn.h:588
#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
float recipvec1_f[3]
Definition: msm_defn.h:586
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
static int bin_evaluation_k_away(NL_Msm *pm)
int maxatoms
Definition: msm_defn.h:600
int * next
Definition: msm_defn.h:601
static int spatial_hashing(NL_Msm *pm)