NAMD
Classes | Typedefs | Functions | Variables
msm_setup.c File Reference
#include "msm_defn.h"

Go to the source code of this file.

Classes

struct  InterpParams_t
 

Typedefs

typedef struct InterpParams_t InterpParams
 

Functions

void NL_msm_cleanup (NL_Msm *pm)
 
static int setup_bins_1away (NL_Msm *pm)
 
static int setup_bins_k_away (NL_Msm *pm)
 
static int setup_cell_vectors (NL_Msm *pm)
 
static int setup_grids (NL_Msm *pm)
 
static int setup_hgrid_1d (NL_Msm *pm, double len, double *hh, int *nn, int *aindex, int *bindex, int isperiodic)
 
static int print_status (NL_Msm *msm)
 
int NL_msm_setup (NL_Msm *pm, double cutoff, double cellvec1[3], double cellvec2[3], double cellvec3[3], double cellcenter[3], int msmflags)
 

Variables

static InterpParams INTERP_PARAMS []
 

Typedef Documentation

typedef struct InterpParams_t InterpParams

Function Documentation

void NL_msm_cleanup ( NL_Msm pm)

Definition at line 7 of file msm_setup.c.

References NL_Msm_t::eh, NL_Msm_t::eh_f, NL_Msm_t::gc, NL_Msm_t::gc_f, GRID_DONE, NL_Msm_t::lyzd, NL_Msm_t::lyzd_f, NL_Msm_t::lzd, NL_Msm_t::lzd_f, NL_Msm_t::maxlevels, NL_Msm_t::msmflags, NL_MSM_COMPUTE_CUDA_GRID_CUTOFF, NL_MSM_COMPUTE_SPREC, NL_Msm_t::qh, and NL_Msm_t::qh_f.

Referenced by NL_msm_destroy().

7  {
8  int i;
9 #ifdef NL_USE_CUDA
11  NL_msm_cuda_cleanup_gridcutoff(pm);
12  }
13 #endif /* NL_USE_CUDA */
14  if (pm->msmflags & NL_MSM_COMPUTE_SPREC) {
15  for (i = 0; i < pm->maxlevels; i++) {
16  GRID_DONE( &(pm->qh_f[i]) );
17  GRID_DONE( &(pm->eh_f[i]) );
18  GRID_DONE( &(pm->gc_f[i]) );
19  }
20  }
21  else {
22  for (i = 0; i < pm->maxlevels; i++) {
23  GRID_DONE( &(pm->qh[i]) );
24  GRID_DONE( &(pm->eh[i]) );
25  GRID_DONE( &(pm->gc[i]) );
26  }
27  }
28  free(pm->lzd);
29  free(pm->lyzd);
30  free(pm->lzd_f);
31  free(pm->lyzd_f);
32  free(pm->qh);
33  free(pm->eh);
34  free(pm->gc);
35  free(pm->qh_f);
36  free(pm->eh_f);
37  free(pm->gc_f);
38 }
int msmflags
Definition: msm_defn.h:590
NL_Msmgrid_double * eh
Definition: msm_defn.h:667
#define GRID_DONE(_p)
Definition: msm_defn.h:61
NL_Msmgrid_double * gc
Definition: msm_defn.h:668
double * lzd
Definition: msm_defn.h:677
float * lyzd_f
Definition: msm_defn.h:681
NL_Msmgrid_float * gc_f
Definition: msm_defn.h:672
double * lyzd
Definition: msm_defn.h:678
NL_Msmgrid_float * eh_f
Definition: msm_defn.h:671
NL_Msmgrid_double * qh
Definition: msm_defn.h:666
NL_Msmgrid_float * qh_f
Definition: msm_defn.h:670
int maxlevels
Definition: msm_defn.h:674
float * lzd_f
Definition: msm_defn.h:680
int NL_msm_setup ( NL_Msm msm,
double  cutoff,
double  cellvec1[3],
double  cellvec2[3],
double  cellvec3[3],
double  cellcenter[3],
int  msmflags 
)

Setup the MSM solver for the molecular system. These parameters come directly from the simulation. Here the cutoff provides some control over the accuracy. The region of space containing the atoms is defined as a parallelepiped, and the "isperiodic" flag determines whether or not each individual cell dimension has periodic boundary conditions. The defined cell is expected to contain the atoms regardless of the choice of boundary conditions. The cell volume must contain the cutoff-sphere.

XXX For now, the MSM solver permits only cell vectors that align with the x-, y-, and z-axis, respectively.

Parameters
pmthe MSM solver object
cutoffcutoff distance for short-range part
cellvec1cell basis vector 1
cellvec2cell basis vector 2
cellvec3cell basis vector 3
cellcentercenter of cell
msmflagsflags for periodicity and calculation

Definition at line 62 of file msm_setup.c.

References NL_Msm_t::a, NL_Msm_t::a_f, NL_Msm_t::cellcenter, NL_Msm_t::cellcenter_f, NL_Msm_t::cellvec1, NL_Msm_t::cellvec1_f, NL_Msm_t::cellvec2, NL_Msm_t::cellvec2_f, NL_Msm_t::cellvec3, NL_Msm_t::cellvec3_f, NL_Msm_t::gx, NL_Msm_t::gx_f, NL_Msm_t::gy, NL_Msm_t::gy_f, NL_Msm_t::gz, NL_Msm_t::gz_f, NL_Msm_t::hx, NL_Msm_t::hx_f, NL_Msm_t::hy, NL_Msm_t::hy_f, NL_Msm_t::hz, NL_Msm_t::hz_f, NL_Msm_t::msmflags, NL_MSM_ALL_FLAGS, NL_MSM_COMPUTE_1AWAY, NL_MSM_COMPUTE_ALL, NL_MSM_COMPUTE_CUDA_FALL_BACK, NL_MSM_COMPUTE_CUDA_GRID_CUTOFF, NL_MSM_COMPUTE_LONG_RANGE, NL_MSM_COMPUTE_SHORT_RANGE, NL_MSM_COMPUTE_SPREC, NL_MSM_ERROR_PARAM, NL_MSM_ERROR_SUPPORT, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_REPORT_TIMINGS, NL_MSM_SUCCESS, print_status(), NL_Msm_t::recipvec1, NL_Msm_t::recipvec1_f, NL_Msm_t::recipvec2, NL_Msm_t::recipvec2_f, NL_Msm_t::recipvec3, NL_Msm_t::recipvec3_f, NL_Msm_t::report_timings, setup_bins_1away(), setup_bins_k_away(), setup_cell_vectors(), and setup_grids().

Referenced by ComputeMsmSerialMgr::recvCoord().

70  {
71 
72  int rc = 0;
73 
74  /* check sanity of msmflags */
75  if ((~NL_MSM_ALL_FLAGS & msmflags) != 0 ||
76  (NL_MSM_COMPUTE_ALL & msmflags) == 0) {
77  return NL_MSM_ERROR_PARAM;
78  }
79 
80  /* report timings? */
81  pm->report_timings = ((msmflags & NL_MSM_REPORT_TIMINGS) != 0);
82 
83  /* for now, permit only orthogonal cell aligned with coordinate axes */
84  if (cellvec1[1] != 0 || cellvec1[2] != 0 ||
85  cellvec2[0] != 0 || cellvec2[2] != 0 ||
86  cellvec3[0] != 0 || cellvec3[1] != 0 ||
87  cellvec1[0] <= 0 || cellvec2[1] <= 0 || cellvec3[2] <= 0) {
88  return NL_MSM_ERROR_SUPPORT;
89  }
90 
91  /* check sanity of cutoff wrt expected cell;
92  * XXX cell widths must be at least the cutoff distance length */
93  if (cutoff <= 0 ||
94  (cutoff > cellvec1[0] && (msmflags & NL_MSM_PERIODIC_VEC1)) ||
95  (cutoff > cellvec2[1] && (msmflags & NL_MSM_PERIODIC_VEC2)) ||
96  (cutoff > cellvec3[2] && (msmflags & NL_MSM_PERIODIC_VEC3))) {
97  return NL_MSM_ERROR_PARAM;
98  }
99 
100  pm->msmflags = msmflags;
101  pm->cellvec1[0] = cellvec1[0];
102  pm->cellvec1[1] = cellvec1[1];
103  pm->cellvec1[2] = cellvec1[2];
104  pm->cellvec2[0] = cellvec2[0];
105  pm->cellvec2[1] = cellvec2[1];
106  pm->cellvec2[2] = cellvec2[2];
107  pm->cellvec3[0] = cellvec3[0];
108  pm->cellvec3[1] = cellvec3[1];
109  pm->cellvec3[2] = cellvec3[2];
110  pm->cellcenter[0] = cellcenter[0];
111  pm->cellcenter[1] = cellcenter[1];
112  pm->cellcenter[2] = cellcenter[2];
113 
114  pm->a = cutoff;
115 
116  rc = setup_cell_vectors(pm);
117  if (rc) return rc;
118 
119  if (msmflags & NL_MSM_COMPUTE_SHORT_RANGE) {
120  /* set up bins for short-range part */
121  if (msmflags & NL_MSM_COMPUTE_1AWAY) {
122  rc = setup_bins_1away(pm);
123  }
124  else {
125  rc = setup_bins_k_away(pm);
126  }
127  if (rc) return rc;
128  }
129 
130  if (msmflags & NL_MSM_COMPUTE_LONG_RANGE) {
131  /* set up grid hierarchy for long-range part */
132  rc = setup_grids(pm);
133  if (rc) return rc;
134  }
135 
136  if (msmflags & NL_MSM_COMPUTE_SPREC) {
137  /* fill out single precision data if needed */
138  pm->cellvec1_f[0] = (float) pm->cellvec1[0];
139  pm->cellvec1_f[1] = (float) pm->cellvec1[1];
140  pm->cellvec1_f[2] = (float) pm->cellvec1[2];
141  pm->cellvec2_f[0] = (float) pm->cellvec2[0];
142  pm->cellvec2_f[1] = (float) pm->cellvec2[1];
143  pm->cellvec2_f[2] = (float) pm->cellvec2[2];
144  pm->cellvec3_f[0] = (float) pm->cellvec3[0];
145  pm->cellvec3_f[1] = (float) pm->cellvec3[1];
146  pm->cellvec3_f[2] = (float) pm->cellvec3[2];
147  pm->cellcenter_f[0] = (float) pm->cellcenter[0];
148  pm->cellcenter_f[1] = (float) pm->cellcenter[1];
149  pm->cellcenter_f[2] = (float) pm->cellcenter[2];
150  pm->recipvec1_f[0] = (float) pm->recipvec1[0];
151  pm->recipvec1_f[1] = (float) pm->recipvec1[1];
152  pm->recipvec1_f[2] = (float) pm->recipvec1[2];
153  pm->recipvec2_f[0] = (float) pm->recipvec2[0];
154  pm->recipvec2_f[1] = (float) pm->recipvec2[1];
155  pm->recipvec2_f[2] = (float) pm->recipvec2[2];
156  pm->recipvec3_f[0] = (float) pm->recipvec3[0];
157  pm->recipvec3_f[1] = (float) pm->recipvec3[1];
158  pm->recipvec3_f[2] = (float) pm->recipvec3[2];
159  pm->hx_f = pm->hx;
160  pm->hy_f = pm->hy;
161  pm->hz_f = pm->hz;
162  pm->a_f = pm->a;
163  pm->gx_f = pm->gx;
164  pm->gy_f = pm->gy;
165  pm->gz_f = pm->gz;
166  }
167 
168  print_status(pm);
169 
170 #ifdef NL_USE_CUDA
171  if (msmflags & NL_MSM_COMPUTE_CUDA_GRID_CUTOFF) {
172  rc = NL_msm_cuda_setup_gridcutoff(pm);
173  if (rc == NL_MSM_SUCCESS) {
174  printf("Using CUDA for grid cutoff computation\n");
175  }
176  else {
177  printf("Unable to set up CUDA for grid cutoff computation\n");
178  if (msmflags & NL_MSM_COMPUTE_CUDA_FALL_BACK) {
179  NL_msm_cuda_cleanup_gridcutoff(pm);
180  printf("Falling back on CPU\n");
181  pm->msmflags &= ~NL_MSM_COMPUTE_CUDA_GRID_CUTOFF;
182  }
183  else return rc;
184  }
185  }
186 #else
187  if (msmflags & NL_MSM_COMPUTE_CUDA_GRID_CUTOFF) {
188  if (msmflags & NL_MSM_COMPUTE_CUDA_FALL_BACK) {
189  printf("Falling back on CPU\n");
190  pm->msmflags &= ~NL_MSM_COMPUTE_CUDA_GRID_CUTOFF;
191  }
192  else return NL_MSM_ERROR_SUPPORT;
193  }
194 #endif /* NL_USE_CUDA */
195 
196  return NL_MSM_SUCCESS;
197 }
float cellcenter_f[3]
Definition: msm_defn.h:585
double gz
Definition: msm_defn.h:647
double recipvec2[3]
Definition: msm_defn.h:579
double recipvec3[3]
Definition: msm_defn.h:580
int msmflags
Definition: msm_defn.h:590
double cellvec2[3]
Definition: msm_defn.h:575
float recipvec2_f[3]
Definition: msm_defn.h:587
float gz_f
Definition: msm_defn.h:648
float cellvec1_f[3]
Definition: msm_defn.h:582
static int setup_grids(NL_Msm *pm)
Definition: msm_setup.c:707
float hy_f
Definition: msm_defn.h:637
double cellvec1[3]
Definition: msm_defn.h:574
double hy
Definition: msm_defn.h:634
float hx_f
Definition: msm_defn.h:637
double cellvec3[3]
Definition: msm_defn.h:576
double cellcenter[3]
Definition: msm_defn.h:577
float hz_f
Definition: msm_defn.h:637
static int setup_bins_k_away(NL_Msm *pm)
Definition: msm_setup.c:420
static int setup_bins_1away(NL_Msm *pm)
Definition: msm_setup.c:318
static int print_status(NL_Msm *msm)
Definition: msm_setup.c:218
double a
Definition: msm_defn.h:635
double gy
Definition: msm_defn.h:647
float cellvec3_f[3]
Definition: msm_defn.h:584
double gx
Definition: msm_defn.h:647
float recipvec3_f[3]
Definition: msm_defn.h:588
double hz
Definition: msm_defn.h:634
static int setup_cell_vectors(NL_Msm *pm)
Definition: msm_setup.c:281
double hx
Definition: msm_defn.h:634
float a_f
Definition: msm_defn.h:638
float recipvec1_f[3]
Definition: msm_defn.h:586
float gx_f
Definition: msm_defn.h:648
float cellvec2_f[3]
Definition: msm_defn.h:583
double recipvec1[3]
Definition: msm_defn.h:578
float gy_f
Definition: msm_defn.h:648
int report_timings
Definition: msm_defn.h:683
int print_status ( NL_Msm msm)
static

Definition at line 218 of file msm_setup.c.

References NL_Msm_t::a, NL_Msm_t::approx, NL_Msm_t::cellvec1, NL_Msm_t::cellvec2, NL_Msm_t::cellvec3, NL_Msm_t::gridspacing, NL_Msm_t::gx, NL_Msm_t::gy, NL_Msm_t::gz, NL_Msm_t::hx, NL_Msm_t::hy, NL_Msm_t::hz, NL_Msm_t::msmflags, NL_msm_approx_name(), NL_MSM_COMPUTE_SPREC, NL_MSM_PERIODIC_ALL, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_msm_split_name(), NL_MSM_SUCCESS, NL_Msm_t::nlevels, NL_Msm_t::nx, NL_Msm_t::ny, NL_Msm_t::nz, NL_Msm_t::qh, NL_Msm_t::qh_f, and NL_Msm_t::split.

Referenced by NL_msm_setup().

218  {
219  int k;
220  int ispx = (pm->msmflags & NL_MSM_PERIODIC_VEC1);
221  int ispy = (pm->msmflags & NL_MSM_PERIODIC_VEC2);
222  int ispz = (pm->msmflags & NL_MSM_PERIODIC_VEC3);
223  int ispany = (pm->msmflags & NL_MSM_PERIODIC_ALL);
224  int ispall = (ispany == NL_MSM_PERIODIC_ALL);
225 
226  const double xlen = pm->cellvec1[0]; /* XXX */
227  const double ylen = pm->cellvec2[1];
228  const double zlen = pm->cellvec3[2];
229 
230  printf("#MSM SETUP\n");
231  printf("# cell lengths= %g %g %g\n", xlen, ylen, zlen);
232  printf("# grid origin= %g %g %g\n", pm->gx, pm->gy, pm->gz);
233  if (ispall) {
234  printf("# periodic boundaries\n");
235  }
236  else if (!ispany) {
237  printf("# non-periodic boundaries\n");
238  }
239  else {
240  printf("# periodic boundaries in:%s%s%s\n",
241  (ispx ? " x" : ""), (ispy ? " y" : ""), (ispz ? " z" : ""));
242  }
243  printf("# cutoff= %g\n", pm->a);
244  printf("# grid spacing= %g\n", pm->gridspacing);
245  printf("# hx, hy, hz= %g, %g, %g\n", pm->hx, pm->hy, pm->hz);
246  printf("# h-grid size= %d x %d x %d\n", pm->nx, pm->ny, pm->nz);
247  printf("# number of levels= %d\n", pm->nlevels);
248  printf("# approximation= %s\n", NL_msm_approx_name(pm->approx));
249  printf("# splitting= %s\n", NL_msm_split_name(pm->split));
250  printf("# grid hierarchy:\n");
251  if (pm->msmflags & NL_MSM_COMPUTE_SPREC) {
252  for (k = 0; k < pm->nlevels; k++) {
253  NL_Msmgrid_float *qh = &(pm->qh_f[k]);
254  int ia = qh->i0;
255  int ib = ia + qh->ni - 1;
256  int ja = qh->j0;
257  int jb = ja + qh->nj - 1;
258  int ka = qh->k0;
259  int kb = ka + qh->nk - 1;
260  printf("# level= %d: [%d..%d] x [%d..%d] x [%d..%d]\n",
261  k, ia, ib, ja, jb, ka, kb);
262  }
263  }
264  else {
265  for (k = 0; k < pm->nlevels; k++) {
266  NL_Msmgrid_double *qh = &(pm->qh[k]);
267  int ia = qh->i0;
268  int ib = ia + qh->ni - 1;
269  int ja = qh->j0;
270  int jb = ja + qh->nj - 1;
271  int ka = qh->k0;
272  int kb = ka + qh->nk - 1;
273  printf("# level= %d: [%d..%d] x [%d..%d] x [%d..%d]\n",
274  k, ia, ib, ja, jb, ka, kb);
275  }
276  }
277  return NL_MSM_SUCCESS;
278 }
const char * NL_msm_split_name(int split)
Definition: msm.c:220
const char * NL_msm_approx_name(int approx)
Definition: msm.c:170
int setup_bins_1away ( NL_Msm pm)
static

Definition at line 318 of file msm_setup.c.

References NL_Msm_t::a, NL_Msm_t::bin, NL_Msm_t::cellvec1, NL_Msm_t::cellvec2, NL_Msm_t::cellvec3, NL_Msm_t::maxbins, NL_Msm_t::msmflags, NL_Msm_t::nbx, NL_Msm_t::nby, NL_Msm_t::nbz, NL_MSM_ERROR_MALLOC, NL_MSM_ERROR_PARAM, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, NL_Msm_t::numbins, NL_Msm_t::recipvec1, NL_Msm_t::recipvec2, NL_Msm_t::recipvec3, X, Y, and Z.

Referenced by NL_msm_setup().

318  {
319  const double *u = pm->cellvec1;
320  const double *v = pm->cellvec2;
321  const double *w = pm->cellvec3;
322  double *ru = pm->recipvec1;
323  double *rv = pm->recipvec2;
324  double *rw = pm->recipvec3;
325  double p[3];
326  double pulen, pvlen, pwlen, s;
327  double cutoff = pm->a;
328  int nbx, nby, nbz, numbins;
329  int ispx = ((pm->msmflags & NL_MSM_PERIODIC_VEC1) != 0);
330  int ispy = ((pm->msmflags & NL_MSM_PERIODIC_VEC2) != 0);
331  int ispz = ((pm->msmflags & NL_MSM_PERIODIC_VEC3) != 0);
332 
333  /* calculate the number of atom bins in each basis vector dimension,
334  * such that each bin (a parallelepiped) inscribes the cutoff-cube;
335  * for periodic boundaries, we have to choose equal sized bins with
336  * length >= cutoff that tile the cell; for non-periodic boundaries,
337  * we can have bins of length = cutoff */
338 
339  /* find the largest orthogonal box inscribed within parallelepiped cell
340  * by taking orthogonal projections onto cross products of basis vectors */
341 
342  p[X] = v[Y]*w[Z] - v[Z]*w[Y]; /* p = v CROSS w */
343  p[Y] = v[Z]*w[X] - v[X]*w[Z];
344  p[Z] = v[X]*w[Y] - v[Y]*w[X];
345  /* s = (u DOT p) / (p DOT p) */
346  s = (u[X]*p[X] + u[Y]*p[Y] + u[Z]*p[Z]) / (p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
347  p[X] *= s; /* p is orthogonal projection of u onto v CROSS w */
348  p[Y] *= s;
349  p[Z] *= s;
350  pulen = sqrt(p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
351 
352  p[X] = w[Y]*u[Z] - w[Z]*u[Y]; /* p = w CROSS u */
353  p[Y] = w[Z]*u[X] - w[X]*u[Z];
354  p[Z] = w[X]*u[Y] - w[Y]*u[X];
355  /* s = (v DOT p) / (p DOT p) */
356  s = (v[X]*p[X] + v[Y]*p[Y] + v[Z]*p[Z]) / (p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
357  p[X] *= s; /* p is orthogonal projection of v onto w CROSS u */
358  p[Y] *= s;
359  p[Z] *= s;
360  pvlen = sqrt(p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
361 
362  p[X] = u[Y]*v[Z] - u[Z]*v[Y]; /* p = u CROSS v */
363  p[Y] = u[Z]*v[X] - u[X]*v[Z];
364  p[Z] = u[X]*v[Y] - u[Y]*v[X];
365  /* s = (w DOT p) / (p DOT p) */
366  s = (w[X]*p[X] + w[Y]*p[Y] + w[Z]*p[Z]) / (p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
367  p[X] *= s; /* p is orthogonal projection of w onto u CROSS v */
368  p[Y] *= s;
369  p[Z] *= s;
370  pwlen = sqrt(p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
371 
372  /* find the number of cells along each dimension by dividing the
373  * orthogonal projection vector lengths by the cutoff distance */
374  nbx = (ispx ? (int) floor(pulen / cutoff) : (int) ceil(pulen / cutoff));
375  nby = (ispy ? (int) floor(pvlen / cutoff) : (int) ceil(pvlen / cutoff));
376  nbz = (ispz ? (int) floor(pwlen / cutoff) : (int) ceil(pwlen / cutoff));
377  if (nbx == 0 || nby == 0 || nbz == 0) {
378  return NL_MSM_ERROR_PARAM; /* cutoff is too large for cell basis */
379  }
380  numbins = nbx * nby * nbz;
381  /* we don't know the number of atoms until compute */
382 
383  /* allocate one atom index for each bin */
384  if (pm->maxbins < numbins) {
385  void *vptr = malloc(numbins * sizeof(int));
386  if (vptr == NULL) return NL_MSM_ERROR_MALLOC;
387  pm->bin = (int *) vptr;
388  pm->maxbins = numbins;
389  }
390  pm->nbx = nbx;
391  pm->nby = nby;
392  pm->nbz = nbz;
393  pm->numbins = numbins;
394 
395  /* must allocate next index array when we know number of atoms */
396 
397  /* rescale recip space vectors for non-periodic expansion */
398  if (ispx == 0) {
399  s = pulen / (nbx * cutoff);
400  ru[X] *= s;
401  ru[Y] *= s;
402  ru[Z] *= s;
403  }
404  if (ispy == 0) {
405  s = pvlen / (nby * cutoff);
406  rv[X] *= s;
407  rv[Y] *= s;
408  rv[Z] *= s;
409  }
410  if (ispz == 0) {
411  s = pwlen / (nbz * cutoff);
412  rw[X] *= s;
413  rw[Y] *= s;
414  rw[Z] *= s;
415  }
416  return NL_MSM_SUCCESS;
417 }
double recipvec2[3]
Definition: msm_defn.h:579
double recipvec3[3]
Definition: msm_defn.h:580
int msmflags
Definition: msm_defn.h:590
double cellvec2[3]
Definition: msm_defn.h:575
int numbins
Definition: msm_defn.h:592
int * bin
Definition: msm_defn.h:597
#define X
Definition: msm_defn.h:29
double cellvec1[3]
Definition: msm_defn.h:574
double cellvec3[3]
Definition: msm_defn.h:576
#define Z
Definition: msm_defn.h:33
double a
Definition: msm_defn.h:635
int maxbins
Definition: msm_defn.h:593
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 setup_bins_k_away ( NL_Msm pm)
static

Definition at line 420 of file msm_setup.c.

References NL_Msm_t::a, NL_Msm_t::bin, NL_Msm_t::binfill, NL_Msm_t::bu, NL_Msm_t::bv, NL_Msm_t::bw, NL_Msm_t::cellvec1, NL_Msm_t::cellvec2, NL_Msm_t::cellvec3, NL_Msm_t::density, NL_Msm_t::maxbins, NL_Msm_t::msmflags, NL_Msm_t::nbinslots, NL_Msm_t::nbrhd, NL_Msm_t::nbrhdlen, NL_Msm_t::nbrhdmax, NL_Msm_t::nbx, NL_Msm_t::nby, NL_Msm_t::nbz, NL_MSM_ERROR_MALLOC, NL_MSM_ERROR_PARAM, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, NL_Msm_t::nradius, NL_Msm_t::numbins, NL_Msm_t::recipvec1, NL_Msm_t::recipvec2, NL_Msm_t::recipvec3, X, Y, and Z.

Referenced by NL_msm_setup().

420  {
421  double *u = pm->cellvec1;
422  double *v = pm->cellvec2;
423  double *w = pm->cellvec3;
424  double *ru = pm->recipvec1;
425  double *rv = pm->recipvec2;
426  double *rw = pm->recipvec3;
427  double *bu = pm->bu;
428  double *bv = pm->bv;
429  double *bw = pm->bw;
430  double p[3];
431  double pulen, pvlen, pwlen, s;
432  double cutoff = pm->a;
433  int nbx, nby, nbz, numbins;
434  int ispx = ((pm->msmflags & NL_MSM_PERIODIC_VEC1) != 0);
435  int ispy = ((pm->msmflags & NL_MSM_PERIODIC_VEC2) != 0);
436  int ispz = ((pm->msmflags & NL_MSM_PERIODIC_VEC3) != 0);
437  int nbrhdmax, nradius, ndiameter, index, i, j, k;
438  double volume;
439  double binvol, abincnt, binlen, c;
440  double min2;
441 
442  /* find the largest orthogonal box inscribed within parallelepiped cell
443  * by taking orthogonal projections onto cross products of basis vectors */
444  p[X] = v[Y]*w[Z] - v[Z]*w[Y]; /* p = v CROSS w */
445  p[Y] = v[Z]*w[X] - v[X]*w[Z];
446  p[Z] = v[X]*w[Y] - v[Y]*w[X];
447  /* s = (u DOT p) / (p DOT p) */
448  s = (u[X]*p[X] + u[Y]*p[Y] + u[Z]*p[Z]) / (p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
449  p[X] *= s; /* p is orthogonal projection of u onto v CROSS w */
450  p[Y] *= s;
451  p[Z] *= s;
452  pulen = sqrt(p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
453 
454  p[X] = w[Y]*u[Z] - w[Z]*u[Y]; /* p = w CROSS u */
455  p[Y] = w[Z]*u[X] - w[X]*u[Z];
456  p[Z] = w[X]*u[Y] - w[Y]*u[X];
457  /* s = (v DOT p) / (p DOT p) */
458  s = (v[X]*p[X] + v[Y]*p[Y] + v[Z]*p[Z]) / (p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
459  p[X] *= s; /* p is orthogonal projection of v onto w CROSS u */
460  p[Y] *= s;
461  p[Z] *= s;
462  pvlen = sqrt(p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
463 
464  p[X] = u[Y]*v[Z] - u[Z]*v[Y]; /* p = u CROSS v */
465  p[Y] = u[Z]*v[X] - u[X]*v[Z];
466  p[Z] = u[X]*v[Y] - u[Y]*v[X];
467  volume = w[X]*p[X] + w[Y]*p[Y] + w[Z]*p[Z];
468  /* s = (w DOT p) / (p DOT p) */
469  s = volume / (p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
470  p[X] *= s; /* p is orthogonal projection of w onto u CROSS v */
471  p[Y] *= s;
472  p[Z] *= s;
473  pwlen = sqrt(p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
474 
475  if ((ispx && cutoff > pulen) ||
476  (ispy && cutoff > pvlen) ||
477  (ispz && cutoff > pwlen)) {
478  printf("Cutoff %.3g too big for cell along basis%s%s%s\n",
479  cutoff, (cutoff > pulen ? " x" : ""), (cutoff > pvlen ? " y" : ""),
480  (cutoff > pwlen ? " z" : ""));
481  return NL_MSM_ERROR_PARAM;
482  }
483 
484  /* calculate the ideal bin volume based on a fixed bin size (nbinslots),
485  * the particle density (density), and a desired fill ratio (binfill) */
486  binvol = pm->binfill * pm->nbinslots / pm->density;
487  abincnt = volume / binvol;
488  binlen = pow(binvol, 1./3); /* ideal side length - use for nonperiodic */
489 
490  /* factor "abincnt" into 3 parts, each part proportional to the
491  * lengths of the orthogonal projection vectors calculated above */
492  c = pow(abincnt / (pulen*pvlen*pwlen), 1./3); /* proportionality const */
493  nbx = (int) ceil(c * pulen);
494  nby = (int) ceil(c * pvlen);
495  nbz = (int) ceil(c * pwlen);
496  numbins = nbx * nby * nbz;
497 
498  printf("nbx=%d nby=%d nbz=%d numbins=%d\n", nbx, nby, nbz, numbins);
499 
500  /* allocate one atom index for each bin */
501  if (pm->maxbins < numbins) {
502  void *vptr = malloc(numbins * sizeof(int));
503  if (vptr == NULL) return NL_MSM_ERROR_MALLOC;
504  pm->bin = (int *) vptr;
505  pm->maxbins = numbins;
506  }
507  pm->nbx = nbx;
508  pm->nby = nby;
509  pm->nbz = nbz;
510  pm->numbins = numbins;
511 
512  /* rescale basis and recip space vectors for non-periodic expansion */
513  if ( ! ispx) {
514  s = pulen / (nbx * binlen);
515  ru[X] *= s;
516  ru[Y] *= s;
517  ru[Z] *= s;
518  s = 1./s;
519  u[X] *= s;
520  u[Y] *= s;
521  u[Z] *= s;
522  }
523  if ( ! ispy) {
524  s = pvlen / (nby * binlen);
525  rv[X] *= s;
526  rv[Y] *= s;
527  rv[Z] *= s;
528  s = 1./s;
529  v[X] *= s;
530  v[Y] *= s;
531  v[Z] *= s;
532  }
533  if ( ! ispz) {
534  s = pwlen / (nbz * binlen);
535  rw[X] *= s;
536  rw[Y] *= s;
537  rw[Z] *= s;
538  s = 1./s;
539  w[X] *= s;
540  w[Y] *= s;
541  w[Z] *= s;
542  }
543 
544  /* scale basis vectors by number of bins to get bin basis vectors */
545  s = 1./nbx;
546  bu[X] = s * u[X];
547  bu[Y] = s * u[Y];
548  bu[Z] = s * u[Z];
549  s = 1./nby;
550  bv[X] = s * v[X];
551  bv[Y] = s * v[Y];
552  bv[Z] = s * v[Z];
553  s = 1./nbz;
554  bw[X] = s * w[X];
555  bw[Y] = s * w[Y];
556  bw[Z] = s * w[Z];
557 
558  /* determine the neighborhood of bins */
559 
560  /* first find minimum width of cell */
561  min2 = (bu[X]*bu[X] + bu[Y]*bu[Y] + bu[Z]*bu[Z]);
562  s = (bv[X]*bv[X] + bv[Y]*bv[Y] + bv[Z]*bv[Z]);
563  if (min2 > s) min2 = s;
564  s = (bw[X]*bw[X] + bw[Y]*bw[Y] + bw[Z]*bw[Z]);
565  if (min2 > s) min2 = s;
566 
567  /* also find the minimum of the four major diagonals */
568  p[X] = bu[X] + bv[X] + bw[X];
569  p[Y] = bu[Y] + bv[Y] + bw[Y];
570  p[Z] = bu[Z] + bv[Z] + bw[Z];
571  s = (p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
572  if (min2 > s) min2 = s;
573  p[X] = -bu[X] + bv[X] + bw[X];
574  p[Y] = -bu[Y] + bv[Y] + bw[Y];
575  p[Z] = -bu[Z] + bv[Z] + bw[Z];
576  s = (p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
577  if (min2 > s) min2 = s;
578  p[X] = bu[X] - bv[X] + bw[X];
579  p[Y] = bu[Y] - bv[Y] + bw[Y];
580  p[Z] = bu[Z] - bv[Z] + bw[Z];
581  s = (p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
582  if (min2 > s) min2 = s;
583  p[X] = bu[X] + bv[X] - bw[X];
584  p[Y] = bu[Y] + bv[Y] - bw[Y];
585  p[Z] = bu[Z] + bv[Z] - bw[Z];
586  s = (p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
587  if (min2 > s) min2 = s;
588 
589  /* the neighborhood is contained in the cube of nradius bins,
590  * store only upper half of bin neighborhood */
591  nradius = (int) ceil(sqrt(cutoff*cutoff / min2));
592  ndiameter = 2 * nradius + 1;
593  nbrhdmax = (ndiameter * ndiameter * ndiameter) / 2 + 1;
594 
595 #if 0
596  printf("Neighborhood radius = %d\n", nradius);
597  printf("min2 = %g\n", min2);
598 #endif
599 
600  /* allocate space for the entire cube */
601  if (pm->nbrhd) free(pm->nbrhd);
602  pm->nbrhd = (int *) calloc(3*nbrhdmax, sizeof(int));
603  if (pm->nbrhd == NULL) return NL_MSM_ERROR_MALLOC;
604  pm->nbrhdmax = nbrhdmax;
605  pm->nradius = nradius; /* save neighborhood radius for diagnostic purposes */
606 
607  /* trim the neighborhood */
608  index = 0;
609  for (k = -nradius; k <= nradius; k++) {
610  for (j = -nradius; j <= nradius; j++) {
611  for (i = -nradius; i <= nradius; i++) {
612  /* XXX should we remove (0,0,0) from bin neighborhood? */
613  int ip, jp, kp, iq, jq, kq;
614  int is_within = 0;
615  int binindex = (k * ndiameter + j) * ndiameter + i;
616  if (binindex < 0) continue;
617  for (kp = 0; kp <= 1; kp++) {
618  for (jp = 0; jp <= 1; jp++) {
619  for (ip = 0; ip <= 1; ip++) {
620  p[X] = (i+ip)*bu[X] + (j+jp)*bv[X] + (k+kp)*bw[X];
621  p[Y] = (i+ip)*bu[Y] + (j+jp)*bv[Y] + (k+kp)*bw[Y];
622  p[Z] = (i+ip)*bu[Z] + (j+jp)*bv[Z] + (k+kp)*bw[Z];
623  for (kq = 0; kq <= 1; kq++) {
624  for (jq = 0; jq <= 1; jq++) {
625  for (iq = 0; iq <= 1; iq++) {
626  double q[3];
627  double dx, dy, dz, r2;
628  q[X] = iq*bu[X] + jq*bv[X] + kq*bw[X];
629  q[Y] = iq*bu[Y] + jq*bv[Y] + kq*bw[Y];
630  q[Z] = iq*bu[Z] + jq*bv[Z] + kq*bw[Z];
631  dx = p[X] - q[X];
632  dy = p[Y] - q[Y];
633  dz = p[Z] - q[Z];
634  r2 = dx*dx + dy*dy + dz*dz;
635  is_within |= (r2 < cutoff*cutoff);
636  }
637  }
638  } /* end q loop */
639 
640  }
641  }
642  } /* end p loop */
643  if (is_within) {
644  pm->nbrhd[index++] = i;
645  pm->nbrhd[index++] = j;
646  pm->nbrhd[index++] = k;
647  }
648 
649  }
650  }
651  }
652  pm->nbrhdlen = index / 3;
653 
654  return NL_MSM_SUCCESS;
655 }
double recipvec2[3]
Definition: msm_defn.h:579
double recipvec3[3]
Definition: msm_defn.h:580
int msmflags
Definition: msm_defn.h:590
double cellvec2[3]
Definition: msm_defn.h:575
double bw[3]
Definition: msm_defn.h:610
int numbins
Definition: msm_defn.h:592
int * bin
Definition: msm_defn.h:597
double bu[3]
Definition: msm_defn.h:608
#define X
Definition: msm_defn.h:29
double cellvec1[3]
Definition: msm_defn.h:574
int nradius
Definition: msm_defn.h:615
double density
Definition: msm_defn.h:605
double cellvec3[3]
Definition: msm_defn.h:576
int nbinslots
Definition: msm_defn.h:606
#define Z
Definition: msm_defn.h:33
double a
Definition: msm_defn.h:635
double binfill
Definition: msm_defn.h:604
int maxbins
Definition: msm_defn.h:593
int nbx
Definition: msm_defn.h:594
int * nbrhd
Definition: msm_defn.h:612
int nbrhdmax
Definition: msm_defn.h:614
double bv[3]
Definition: msm_defn.h:609
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 nbrhdlen
Definition: msm_defn.h:613
int setup_cell_vectors ( NL_Msm pm)
static

Definition at line 281 of file msm_setup.c.

References NL_Msm_t::cellvec1, NL_Msm_t::cellvec2, NL_Msm_t::cellvec3, NL_MSM_SUCCESS, NL_Msm_t::recipvec1, NL_Msm_t::recipvec2, NL_Msm_t::recipvec3, X, Y, and Z.

Referenced by NL_msm_setup().

281  {
282  double *u = pm->cellvec1;
283  double *v = pm->cellvec2;
284  double *w = pm->cellvec3;
285  double *bu = pm->recipvec1;
286  double *bv = pm->recipvec2;
287  double *bw = pm->recipvec3;
288  double c[3], s;
289 
290  c[X] = v[Y]*w[Z] - v[Z]*w[Y]; /* v CROSS w */
291  c[Y] = v[Z]*w[X] - v[X]*w[Z];
292  c[Z] = v[X]*w[Y] - v[Y]*w[X];
293  s = 1 / (u[X]*c[X] + u[Y]*c[Y] + u[Z]*c[Z]); /* 1 / (c DOT u) */
294  bu[X] = s*c[X];
295  bu[Y] = s*c[Y];
296  bu[Z] = s*c[Z];
297 
298  c[X] = w[Y]*u[Z] - w[Z]*u[Y]; /* w CROSS u */
299  c[Y] = w[Z]*u[X] - w[X]*u[Z];
300  c[Z] = w[X]*u[Y] - w[Y]*u[X];
301  s = 1 / (v[X]*c[X] + v[Y]*c[Y] + v[Z]*c[Z]); /* 1 / (c DOT v) */
302  bv[X] = s*c[X];
303  bv[Y] = s*c[Y];
304  bv[Z] = s*c[Z];
305 
306  c[X] = u[Y]*v[Z] - u[Z]*v[Y]; /* u CROSS v */
307  c[Y] = u[Z]*v[X] - u[X]*v[Z];
308  c[Z] = u[X]*v[Y] - u[Y]*v[X];
309  s = 1 / (w[X]*c[X] + w[Y]*c[Y] + w[Z]*c[Z]); /* 1 / (c DOT w) */
310  bw[X] = s*c[X];
311  bw[Y] = s*c[Y];
312  bw[Z] = s*c[Z];
313 
314  return NL_MSM_SUCCESS;
315 }
double recipvec2[3]
Definition: msm_defn.h:579
double recipvec3[3]
Definition: msm_defn.h:580
double cellvec2[3]
Definition: msm_defn.h:575
#define X
Definition: msm_defn.h:29
double cellvec1[3]
Definition: msm_defn.h:574
double cellvec3[3]
Definition: msm_defn.h:576
#define Z
Definition: msm_defn.h:33
#define Y
Definition: msm_defn.h:31
double recipvec1[3]
Definition: msm_defn.h:578
int setup_grids ( NL_Msm pm)
static

Definition at line 707 of file msm_setup.c.

References NL_Msm_t::a, NL_Msm_t::approx, ASSERT, NL_Msm_t::cellcenter, NL_Msm_t::cellvec1, NL_Msm_t::cellvec2, NL_Msm_t::cellvec3, NL_Msm_t::eh, NL_Msm_t::eh_f, NL_Msm_t::gc, NL_Msm_t::gc_f, GRID_INDEX, GRID_INDEX_CHECK, GRID_INIT, GRID_RESIZE, NL_Msm_t::gx, NL_Msm_t::gy, NL_Msm_t::gz, NL_Msm_t::gzero, NL_Msm_t::gzero_f, NL_Msm_t::hx, NL_Msm_t::hy, NL_Msm_t::hz, NL_Msm_t::lyzd, NL_Msm_t::lyzd_f, NL_Msm_t::lzd, NL_Msm_t::lzd_f, NL_Msm_t::max_lyzd, NL_Msm_t::max_lzd, NL_Msm_t::maxlevels, NL_Msm_t::msmflags, NL_MSM_COMPUTE_SPREC, NL_MSM_ERROR_MALLOC, NL_MSM_PERIODIC_ALL, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, NL_Msm_t::nlevels, InterpParams_t::nu, NL_Msm_t::nx, NL_Msm_t::ny, NL_Msm_t::nz, InterpParams_t::omega, NL_Msm_t::qh, NL_Msm_t::qh_f, setup_hgrid_1d(), split(), NL_Msm_t::split, and SPOLY.

Referenced by NL_msm_setup().

707  {
708  const int nu = INTERP_PARAMS[pm->approx].nu;
709  const int omega = INTERP_PARAMS[pm->approx].omega;
710  const int split = pm->split;
711  const int ispx = (pm->msmflags & NL_MSM_PERIODIC_VEC1);
712  const int ispy = (pm->msmflags & NL_MSM_PERIODIC_VEC2);
713  const int ispz = (pm->msmflags & NL_MSM_PERIODIC_VEC3);
714  const int ispany = (pm->msmflags & NL_MSM_PERIODIC_ALL);
715 
716  const int issprec = (pm->msmflags & NL_MSM_COMPUTE_SPREC);
717 
718  const double xlen = pm->cellvec1[0]; /* XXX */
719  const double ylen = pm->cellvec2[1];
720  const double zlen = pm->cellvec3[2];
721 
722  const double a = pm->a;
723  double hx, hy, hz;
724  double scaling;
725  double d; /* temporary for SPOLY derivative */
726 
727  NL_Msmgrid_double *p = NULL;
728  NL_Msmgrid_float *p_f = NULL;
729  int ia, ib, ja, jb, ka, kb, ni, nj, nk;
730  int nx, ny, nz; /* counts the grid points that span just the domain */
731 
732  int i, j, k, n;
733  int index;
734  int level, toplevel, nlevels, maxlevels;
735  int lastnelems = 1;
736  int isclamped = 0;
737  int done, alldone;
738 
739  int rc = 0;
740 
741  rc = setup_hgrid_1d(pm, xlen, &hx, &nx, &ia, &ib, ispx);
742  if (rc) return rc;
743 
744  rc = setup_hgrid_1d(pm, ylen, &hy, &ny, &ja, &jb, ispy);
745  if (rc) return rc;
746 
747  rc = setup_hgrid_1d(pm, zlen, &hz, &nz, &ka, &kb, ispz);
748  if (rc) return rc;
749 
750  pm->hx = hx;
751  pm->hy = hy;
752  pm->hz = hz;
753 
754  /* XXX set coordinate for h-grid (0,0,0) point */
755  pm->gx = pm->cellcenter[0] - ((nx >> 1) * hx);
756  pm->gy = pm->cellcenter[1] - ((ny >> 1) * hy);
757  pm->gz = pm->cellcenter[2] - ((nz >> 1) * hz);
758 
759  pm->nx = nx;
760  pm->ny = ny;
761  pm->nz = nz;
762 
763  ni = ib - ia + 1;
764  nj = jb - ja + 1;
765  nk = kb - ka + 1;
766 
767  /* allocate temp buffer space for factored grid transfer */
768  n = (nk > omega ? nk : omega); /* row along z-dimension */
769  if (pm->max_lzd < n) {
770  if (issprec) {
771  float *t;
772  t = (float *) realloc(pm->lzd, n * sizeof(float));
773  if (NULL == t) return NL_MSM_ERROR_MALLOC;
774  pm->lzd_f = t;
775  }
776  else {
777  double *t;
778  t = (double *) realloc(pm->lzd, n * sizeof(double));
779  if (NULL == t) return NL_MSM_ERROR_MALLOC;
780  pm->lzd = t;
781  }
782  pm->max_lzd = n;
783  }
784  n *= (nj > omega ? nj : omega); /* plane along yz-dimensions */
785  if (pm->max_lyzd < n) {
786  if (issprec) {
787  float *t;
788  t = (float *) realloc(pm->lyzd, n * sizeof(float));
789  if (NULL == t) return NL_MSM_ERROR_MALLOC;
790  pm->lyzd_f = t;
791  }
792  else {
793  double *t;
794  t = (double *) realloc(pm->lyzd, n * sizeof(double));
795  if (NULL == t) return NL_MSM_ERROR_MALLOC;
796  pm->lyzd = t;
797  }
798  pm->max_lyzd = n;
799  }
800 
801  nlevels = pm->nlevels;
802  if (nlevels <= 0) {
803  /* automatically set number of levels */
804  n = ni;
805  if (n < nj) n = nj;
806  if (n < nk) n = nk;
807  for (maxlevels = 1; n > 0; n >>= 1) maxlevels++;
808  nlevels = maxlevels;
809  if (ispany == 0) { /* no periodicity */
810  int omega3 = omega * omega * omega;
811  int nhalf = (int) sqrtf(ni*nj*nk); /* scale down for performance? */
812  lastnelems = (nhalf > omega3 ? nhalf : omega3);
813  isclamped = 1;
814  }
815  }
816  else {
817  /* user-defined number of levels */
818  maxlevels = nlevels;
819  }
820 
821  /* allocate any additional levels that may be needed */
822  if (pm->maxlevels < maxlevels) {
823  void *vqh, *veh, *vgc;
824  if (issprec) {
825  vqh = realloc(pm->qh, maxlevels * sizeof(NL_Msmgrid_float));
826  if (NULL == vqh) return NL_MSM_ERROR_MALLOC;
827  veh = realloc(pm->eh, maxlevels * sizeof(NL_Msmgrid_float));
828  if (NULL == veh) return NL_MSM_ERROR_MALLOC;
829  vgc = realloc(pm->gc, maxlevels * sizeof(NL_Msmgrid_float));
830  if (NULL == vgc) return NL_MSM_ERROR_MALLOC;
831  pm->qh_f = (NL_Msmgrid_float *) vqh;
832  pm->eh_f = (NL_Msmgrid_float *) veh;
833  pm->gc_f = (NL_Msmgrid_float *) vgc;
834  /* initialize the newest grids appended to array */
835  for (level = pm->maxlevels; level < maxlevels; level++) {
836  GRID_INIT( &(pm->qh_f[level]) );
837  GRID_INIT( &(pm->eh_f[level]) );
838  GRID_INIT( &(pm->gc_f[level]) );
839  }
840  }
841  else {
842  vqh = realloc(pm->qh, maxlevels * sizeof(NL_Msmgrid_double));
843  if (NULL == vqh) return NL_MSM_ERROR_MALLOC;
844  veh = realloc(pm->eh, maxlevels * sizeof(NL_Msmgrid_double));
845  if (NULL == veh) return NL_MSM_ERROR_MALLOC;
846  vgc = realloc(pm->gc, maxlevels * sizeof(NL_Msmgrid_double));
847  if (NULL == vgc) return NL_MSM_ERROR_MALLOC;
848  pm->qh = (NL_Msmgrid_double *) vqh;
849  pm->eh = (NL_Msmgrid_double *) veh;
850  pm->gc = (NL_Msmgrid_double *) vgc;
851  /* initialize the newest grids appended to array */
852  for (level = pm->maxlevels; level < maxlevels; level++) {
853  GRID_INIT( &(pm->qh[level]) );
854  GRID_INIT( &(pm->eh[level]) );
855  GRID_INIT( &(pm->gc[level]) );
856  }
857  }
858  pm->maxlevels = maxlevels;
859  }
860 
861  level = 0;
862  done = 0;
863  alldone = 0;
864  do {
865  if (issprec) {
866  GRID_RESIZE( &(pm->qh_f[level]), float, ia, ni, ja, nj, ka, nk);
867  GRID_RESIZE( &(pm->eh_f[level]), float, ia, ni, ja, nj, ka, nk);
868  }
869  else {
870  GRID_RESIZE( &(pm->qh[level]), double, ia, ni, ja, nj, ka, nk);
871  GRID_RESIZE( &(pm->eh[level]), double, ia, ni, ja, nj, ka, nk);
872  }
873 
874  if (++level == nlevels) done |= 0x07; /* user limit on levels */
875 
876  alldone = (done == 0x07); /* make sure all dimensions are done */
877 
878  if (isclamped) {
879  int nelems = ni * nj * nk;
880  if (nelems <= lastnelems) done |= 0x07;
881  }
882 
883  if (ispx) {
884  ni >>= 1;
885  ib = ni-1;
886  if (ni & 1) done |= 0x07; /* == 3 or 1 */
887  else if (ni == 2) done |= 0x01; /* can do one more */
888  }
889  else {
890  ia = -((-ia+1)/2) - nu;
891  ib = (ib+1)/2 + nu;
892  ni = ib - ia + 1;
893  if (ni <= omega) done |= 0x01; /* can do more restrictions */
894  }
895 
896  if (ispy) {
897  nj >>= 1;
898  jb = nj-1;
899  if (nj & 1) done |= 0x07; /* == 3 or 1 */
900  else if (nj == 2) done |= 0x02; /* can do one more */
901  }
902  else {
903  ja = -((-ja+1)/2) - nu;
904  jb = (jb+1)/2 + nu;
905  nj = jb - ja + 1;
906  if (nj <= omega) done |= 0x02; /* can do more restrictions */
907  }
908 
909  if (ispz) {
910  nk >>= 1;
911  kb = nk-1;
912  if (nk & 1) done |= 0x07; /* == 3 or 1 */
913  else if (nk == 2) done |= 0x04; /* can do one more */
914  }
915  else {
916  ka = -((-ka+1)/2) - nu;
917  kb = (kb+1)/2 + nu;
918  nk = kb - ka + 1;
919  if (nk <= omega) done |= 0x04; /* can do more restrictions */
920  }
921 
922  } while ( ! alldone );
923  pm->nlevels = level;
924 
925  toplevel = (ispany ? pm->nlevels : pm->nlevels - 1);
926 
927  /* ellipsoid axes for grid cutoff weights */
928  ni = (int) ceil(2*a/hx) - 1;
929  nj = (int) ceil(2*a/hy) - 1;
930  nk = (int) ceil(2*a/hz) - 1;
931  scaling = 1; /* initially there is no scaling */
932  for (level = 0; level < toplevel; level++) {
933  if (issprec) {
934  p_f = &(pm->gc_f[level]);
935  GRID_RESIZE(p_f, float, -ni, 2*ni+1, -nj, 2*nj+1, -nk, 2*nk+1);
936  }
937  else {
938  p = &(pm->gc[level]);
939  GRID_RESIZE(p, double, -ni, 2*ni+1, -nj, 2*nj+1, -nk, 2*nk+1);
940  }
941 
942  if (0 == level) {
943  index = 0;
944  for (k = -nk; k <= nk; k++) {
945  for (j = -nj; j <= nj; j++) {
946  for (i = -ni; i <= ni; i++) {
947  double s, t, gs, gt, g;
948  s = sqrt((i*hx)*(i*hx) + (j*hy)*(j*hy) + (k*hz)*(k*hz)) / a;
949  t = 0.5 * s;
950  if (t >= 1) {
951  g = 0;
952  }
953  else if (s >= 1) {
954  gs = 1/s;
955  SPOLY(&gt, &d, t, split);
956  g = (gs - 0.5 * gt) / a;
957  }
958  else {
959  SPOLY(&gs, &d, s, split);
960  SPOLY(&gt, &d, t, split);
961  g = (gs - 0.5 * gt) / a;
962  }
963  if (issprec) {
964  GRID_INDEX_CHECK(p_f, i, j, k);
965  ASSERT(p_f->buffer + index == p_f->data + GRID_INDEX(p_f,i,j,k));
966  p_f->buffer[index] = (float) g;
967  }
968  else {
969  GRID_INDEX_CHECK(p, i, j, k);
970  ASSERT( p->buffer + index == p->data + GRID_INDEX(p, i, j, k) );
971  p->buffer[index] = g;
972  }
973  index++;
974  }
975  }
976  } /* end loops over k-j-i */
977  }
978  else {
979  /* set each level as scaling of h-level */
980  if (issprec) {
981  const NL_Msmgrid_float *first = &(pm->gc_f[0]);
982  index = 0;
983  for (k = -nk; k <= nk; k++) {
984  for (j = -nj; j <= nj; j++) {
985  for (i = -ni; i <= ni; i++) {
986  GRID_INDEX_CHECK(p_f, i, j, k);
987  ASSERT(p_f->buffer + index == p_f->data + GRID_INDEX(p_f,i,j,k));
988  p_f->buffer[index] = (float) (scaling * first->buffer[index]);
989  index++;
990  }
991  }
992  } /* for loops */
993  } /* if isprec */
994  else {
995  const NL_Msmgrid_double *first = &(pm->gc[0]);
996  index = 0;
997  for (k = -nk; k <= nk; k++) {
998  for (j = -nj; j <= nj; j++) {
999  for (i = -ni; i <= ni; i++) {
1000  GRID_INDEX_CHECK(p, i, j, k);
1001  ASSERT( p->buffer + index == p->data + GRID_INDEX(p, i, j, k) );
1002  p->buffer[index] = scaling * first->buffer[index];
1003  index++;
1004  }
1005  }
1006  } /* for loops */
1007  } /* else ! issprec */
1008  }
1009  scaling *= 0.5; /* adjust scaling here to also accommodate toplevel */
1010  } /* end loop over levels */
1011 
1012  if (toplevel < pm->nlevels) {
1013  /* nonperiodic in all dimensions,
1014  * calculate top level weights, ellipsoid axes are length of grid */
1015  if (issprec) {
1016  const NL_Msmgrid_float *qhtop_f = &(pm->qh_f[toplevel]);
1017  ni = qhtop_f->ni - 1;
1018  nj = qhtop_f->nj - 1;
1019  nk = qhtop_f->nk - 1;
1020  }
1021  else {
1022  const NL_Msmgrid_double *qhtop = &(pm->qh[toplevel]);
1023  ni = qhtop->ni - 1;
1024  nj = qhtop->nj - 1;
1025  nk = qhtop->nk - 1;
1026  }
1027  if (issprec) {
1028  p_f = &(pm->gc_f[toplevel]);
1029  GRID_RESIZE(p_f, float, -ni, 2*ni+1, -nj, 2*nj+1, -nk, 2*nk+1);
1030  }
1031  else {
1032  p = &(pm->gc[toplevel]);
1033  GRID_RESIZE(p, double, -ni, 2*ni+1, -nj, 2*nj+1, -nk, 2*nk+1);
1034  }
1035  index = 0;
1036  for (k = -nk; k <= nk; k++) {
1037  for (j = -nj; j <= nj; j++) {
1038  for (i = -ni; i <= ni; i++) {
1039  double s, gs;
1040  s = sqrt((i*hx)*(i*hx) + (j*hy)*(j*hy) + (k*hz)*(k*hz)) / a;
1041  if (s >= 1) {
1042  gs = 1/s;
1043  }
1044  else {
1045  SPOLY(&gs, &d, s, split);
1046  }
1047  if (issprec) {
1048  GRID_INDEX_CHECK(p_f, i, j, k);
1049  ASSERT(p_f->buffer + index == p_f->data + GRID_INDEX(p_f,i,j,k));
1050  p_f->buffer[index] = (float) (scaling * gs/a);
1051  }
1052  else {
1053  GRID_INDEX_CHECK(p, i, j, k);
1054  ASSERT( p->buffer + index == p->data + GRID_INDEX(p, i, j, k) );
1055  p->buffer[index] = scaling * gs/a;
1056  }
1057  index++;
1058  }
1059  }
1060  } /* end loops over k-j-i for coarsest level weights */
1061  }
1062 
1063  /* calculate self energy factor for splitting */
1064  if (1) {
1065  double s, gs;
1066  s = 0;
1067  SPOLY(&gs, &d, s, split);
1068  if (issprec) {
1069  pm->gzero_f = (float) (gs/a);
1070  }
1071  else {
1072  pm->gzero = gs/a;
1073  }
1074  }
1075 
1076  return NL_MSM_SUCCESS;
1077 }
double gz
Definition: msm_defn.h:647
double gzero
Definition: msm_defn.h:650
int msmflags
Definition: msm_defn.h:590
#define GRID_RESIZE(_p, TYPE, __i0, __ni, __j0, __nj, __k0, __nk)
Definition: msm_defn.h:77
static InterpParams INTERP_PARAMS[]
Definition: msm_setup.c:206
NL_Msmgrid_double * eh
Definition: msm_defn.h:667
NL_Msmgrid_double * gc
Definition: msm_defn.h:668
double cellvec2[3]
Definition: msm_defn.h:575
double * lzd
Definition: msm_defn.h:677
int approx
Definition: msm_defn.h:642
#define GRID_INDEX_CHECK(a, _i, _j, _k)
Definition: msm_defn.h:114
int ny
Definition: msm_defn.h:640
double cellvec1[3]
Definition: msm_defn.h:574
double hy
Definition: msm_defn.h:634
double cellvec3[3]
Definition: msm_defn.h:576
#define GRID_INIT(_p)
Definition: msm_defn.h:55
float * lyzd_f
Definition: msm_defn.h:681
NL_Msmgrid_float * gc_f
Definition: msm_defn.h:672
#define GRID_INDEX(_p, _i, _j, _k)
Definition: msm_defn.h:66
double * lyzd
Definition: msm_defn.h:678
double cellcenter[3]
Definition: msm_defn.h:577
#define ASSERT(E)
int nx
Definition: msm_defn.h:640
int max_lyzd
Definition: msm_defn.h:676
double a
Definition: msm_defn.h:635
NL_Msmgrid_float * eh_f
Definition: msm_defn.h:671
NL_Msmgrid_double * qh
Definition: msm_defn.h:666
double gy
Definition: msm_defn.h:647
double gx
Definition: msm_defn.h:647
std::vector< std::string > split(const std::string &text, std::string delimiter)
Definition: MoleculeQM.C:73
double hz
Definition: msm_defn.h:634
int split
Definition: msm_defn.h:643
static int setup_hgrid_1d(NL_Msm *pm, double len, double *hh, int *nn, int *aindex, int *bindex, int isperiodic)
Definition: msm_setup.c:658
double hx
Definition: msm_defn.h:634
NL_Msmgrid_float * qh_f
Definition: msm_defn.h:670
#define SPOLY(pg, pdg, ra, split)
Definition: msm_defn.h:140
int maxlevels
Definition: msm_defn.h:674
int nlevels
Definition: msm_defn.h:645
int nz
Definition: msm_defn.h:640
float gzero_f
Definition: msm_defn.h:651
float * lzd_f
Definition: msm_defn.h:680
int max_lzd
Definition: msm_defn.h:676
int setup_hgrid_1d ( NL_Msm pm,
double  len,
double *  hh,
int *  nn,
int *  aindex,
int *  bindex,
int  isperiodic 
)
static

Definition at line 658 of file msm_setup.c.

References NL_Msm_t::approx, ASSERT, NL_Msm_t::gridspacing, NL_MSM_ERROR_PARAM, NL_MSM_SUCCESS, and InterpParams_t::nu.

Referenced by setup_grids().

666  {
667 
668  const int nu = INTERP_PARAMS[pm->approx].nu; /* interp stencil radius */
669 
670  ASSERT(hmax > 0);
671  if (isperiodic) {
672  const double hmin = (4./5) * pm->gridspacing; /* minimum bound on h */
673  const double hmax = 1.5 * hmin;
674  double h = len;
675  int n = 1; /* start with one grid point across domain */
676  while (h >= hmax) {
677  h *= 0.5; /* halve h */
678  n <<= 1; /* double grid points */
679  }
680  if (h < hmin) {
681  if (n < 4) { /* either len is too small or hmin is too large */
682  return NL_MSM_ERROR_PARAM;
683  }
684  h *= (4./3); /* scale h by 4/3 */
685  n >>= 2; /* scale n by 3/4 */
686  n *= 3;
687  }
688  /* now we have: hmin <= h < hmax */
689  /* now we have: n is power of two times no more than one power of 3 */
690  *hh = h;
691  *nn = n;
692  *aindex = 0;
693  *bindex = n-1;
694  }
695  else { /* non-periodic */
696  double h = pm->gridspacing;
697  int n = (int) floorf(len / h) + 1;
698  *hh = h;
699  *nn = n;
700  *aindex = -nu;
701  *bindex = n + nu;
702  }
703  return NL_MSM_SUCCESS;
704 }
static InterpParams INTERP_PARAMS[]
Definition: msm_setup.c:206
int approx
Definition: msm_defn.h:642
#define ASSERT(E)
double gridspacing
Definition: msm_defn.h:633

Variable Documentation

InterpParams INTERP_PARAMS[]
static
Initial value:
= {
{ 1, 4, 6 },
{ 2, 6, 10 },
{ 2, 6, 10 },
{ 3, 8, 14 },
{ 3, 8, 14 },
{ 4, 10, 18 },
{ 4, 10, 18 },
{ 1, 4, 6 },
}

Definition at line 206 of file msm_setup.c.