NAMD
GridForceGrid.C
Go to the documentation of this file.
1 
7 #include <iostream>
8 #include <typeinfo>
9 
10 #include "GridForceGrid.h"
11 #include "Vector.h"
12 #include "Node.h"
13 #include "SimParameters.h"
14 #include "Molecule.h"
15 #include "InfoStream.h"
16 #include "common.h"
17 #include "ComputeGridForce.h"
18 
19 #include "MGridforceParams.h"
20 
21 #define MIN_DEBUG_LEVEL 4
22 //#define DEBUGM
23 #include "Debug.h"
24 
25 #include "GridForceGrid.inl"
26 
27 
28 /*****************/
29 /* GRIDFORCEGRID */
30 /*****************/
31 
32 // Class methods
33 
34 GridforceGrid * GridforceGrid::new_grid(int gridnum, char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams)
35 {
36  GridforceGrid *grid = NULL;
37  if (mgridParams->gridforceLite) {
38  grid = new GridforceLiteGrid(gridnum);
39  } else {
40  grid = new GridforceFullMainGrid(gridnum);
41  }
42 
43  grid->initialize(potfilename, simParams, mgridParams);
44 
45  return grid;
46 }
47 
49 
51 {
52  // Abstract interface for packing a grid into a message. This
53  // could easily be a non-static function as it was before, but is
54  // static to make it similar to unpack_grid below, which does need
55  // to be static since it creates a GridforceGrid object.
56  msg->put(grid->get_grid_type());
57  grid->pack(msg);
58 }
59 
61 {
62  // Abstract interface for unpacking a grid from a message.
63  GridforceGrid *grid = NULL;
64  int type;
65 
66  msg->get(type);
67 
68  switch (type) {
70  grid = new GridforceFullMainGrid(gridnum);
71  break;
73  grid = new GridforceLiteGrid(gridnum);
74  break;
75  default:
76  NAMD_bug("GridforceGrid::unpack_grid called with unknown grid type!");
77  }
78 
79  grid->unpack(msg);
80 
81  return grid;
82 }
83 
85 {
86  // Idea: Go through each grid corner and wrap it to the grid center;
87  // if the position moves, then the grid is too big and we return false
88  DebugM(4, "Checking whether grid fits periodic cell\n");
89  Position center = get_center();
90  for (int i = 0; i < 8; i++) {
91  Position pos = get_corner(i);
92  Position pos_wrapped = wrap_position(pos, lattice);
93  if ((pos - pos_wrapped).length() > 1.) {
94  DebugM(5, "(" << pos << ") != (" << pos_wrapped << ")\n" << endi);
95  return false;
96  }
97  }
98  return true;
99 }
100 
102 {
103  // idx -> (x,y,z) (cell basis coordinates)
104  // 0 -> (0,0,0)
105  // 1 -> (1,0,0)
106  // 2 -> (0,1,0)
107  // 3 -> (1,1,0)
108  // 4 -> (0,0,1)
109  // 5 -> (1,0,1)
110  // 6 -> (0,1,1)
111  // 7 -> (1,1,1)
112  Position pos;
113  if (idx >= 8 || idx < 0) {
114  // invalid index
115  pos = Vector(); // default value of Vector() is strange enough to be a decent "undefined" value (-99999, -99999, -999999)
116  } else if (corners[idx] != Vector()) {
117  // use cached value if possible
118  pos = corners[idx];
119  } else {
120  // must calculate
121  Tensor e = get_e();
122  pos = get_origin();
123  if (idx & (1 << 0)) pos += e * Vector(get_k0()-1, 0, 0);
124  if (idx & (1 << 1)) pos += e * Vector(0, get_k1()-1, 0);
125  if (idx & (1 << 2)) pos += e * Vector(0, 0, get_k2()-1);
126  corners[idx] = pos; // cache for future use
127  DebugM(4, "corner " << idx << " = " << pos << "\n" << endi);
128  }
129  return pos;
130 }
131 
132 
133 /*************************/
134 /* GRIDFORCEFULLBASEGRID */
135 /*************************/
136 
138 {
139  cont[0] = cont[1] = cont[2] = FALSE;
140  grid = NULL;
141  numSubgrids = 0;
142  subgrids = NULL;
143 }
144 
146 {
147  delete[] grid;
148  for (int i = 0; i < numSubgrids; i++) {
149  delete subgrids[i];
150  }
151  delete[] subgrids;
152 }
153 
154 
156 {
157  DebugM(2, "Packing message\n" << endi);
158 
159  msg->put(numSubgrids);
160  msg->put(generation);
161 
162  msg->put(3*sizeof(int), (char*)k);
163  msg->put(3*sizeof(int), (char*)k_nopad);
164  msg->put(size);
165  msg->put(size_nopad);
166  msg->put(3*sizeof(long int), (char*)dk);
167  msg->put(3*sizeof(long int), (char*)dk_nopad);
168  msg->put(factor);
169 
170  msg->put(sizeof(Vector), (char*)&origin);
171  msg->put(sizeof(Vector), (char*)&center);
172  msg->put(sizeof(Tensor), (char*)&e);
173  msg->put(sizeof(Tensor), (char*)&inv);
174 
175 // msg->put(3*sizeof(float), (char*)pad_p);
176 // msg->put(3*sizeof(float), (char*)pad_n);
177  msg->put(3*sizeof(Bool), (char*)cont);
178  msg->put(3*sizeof(float), (char*)offset);
179  msg->put(3*sizeof(float), (char*)gap);
180  msg->put(3*sizeof(float), (char*)gapinv);
181  msg->put(sizeof(Vector), (char*)&scale);
182  msg->put(sizeof(Bool), (char*)&checksize);
183 
184  DebugM(2, "Packing grid, size = " << size << "\n" << endi);
185 
186  msg->put(size*sizeof(float), (char*)grid);
187 
188  DebugM(2, "Packing subgrids\n" << endi);
189 
190  for (int i = 0; i < numSubgrids; i++) {
191  subgrids[i]->pack(msg);
192  }
193 }
194 
196 {
197  DebugM(3, "Unpacking message\n" << endi);
198 // iout << iINFO << CkMyPe() << " Unpacking message\n" << endi;
199 
200  delete[] grid;
201  grid = NULL;
202  for (int i = 0; i < numSubgrids; i++) {
203  delete subgrids[i];
204  }
205  numSubgrids = 0;
206  delete[] subgrids;
207  subgrids = NULL;
208 
209  msg->get(numSubgrids);
210  msg->get(generation);
211 
212  DebugM(3, "numSubgrids = " << numSubgrids << "\n");
213  DebugM(3, "generation = " << generation << "\n" << endi);
214 
215  msg->get(3*sizeof(int), (char*)k);
216  msg->get(3*sizeof(int), (char*)k_nopad);
217  msg->get(size);
218  msg->get(size_nopad);
219  msg->get(3*sizeof(long int), (char*)dk);
220  msg->get(3*sizeof(long int), (char*)dk_nopad);
221  msg->get(factor);
222 
223  DebugM(3, "size = " << size << "\n" << endi);
224 
225  msg->get(sizeof(Vector), (char*)&origin);
226  msg->get(sizeof(Vector), (char*)&center);
227  msg->get(sizeof(Tensor), (char*)&e);
228  msg->get(sizeof(Tensor), (char*)&inv);
229 
230 // msg->get(3*sizeof(float), (char*)pad_p);
231 // msg->get(3*sizeof(float), (char*)pad_n);
232  msg->get(3*sizeof(Bool), (char*)cont);
233  msg->get(3*sizeof(float), (char*)offset);
234  msg->get(3*sizeof(float), (char*)gap);
235  msg->get(3*sizeof(float), (char*)gapinv);
236  msg->get(sizeof(Vector), (char*)&scale);
237  msg->get(sizeof(Bool), (char*)&checksize);
238 
239  if (size) {
240  DebugM(3, "allocating grid, size = " << size << "\n" << endi);
241  grid = new float[size];
242  msg->get(size*sizeof(float), (char*)grid);
243  }
244 
245  if (numSubgrids) {
246  DebugM(3, "Creating subgrids array, size " << numSubgrids << "\n" << endi);
248  for (int i = 0; i < numSubgrids; i++) {
249  subgrids[i] = new GridforceFullSubGrid(this);
250  subgrids[i]->unpack(msg);
251  }
252  }
253 }
254 
255 
257 {
258  char line[256];
259  long int poten_offset;
260  do {
261  poten_offset = ftell(poten_fp);
262  fgets(line, 256, poten_fp); // Read comment lines
263  DebugM(4, "Read line: " << line << endi);
264  } while (line[0] == '#');
265  fseek(poten_fp, poten_offset, SEEK_SET);
266 
267  // read grid dimensions
268  fscanf(poten_fp, "object %*d class gridpositions counts %d %d %d\n",
269  &k_nopad[0], &k_nopad[1], &k_nopad[2]);
270  size_nopad = k_nopad[0] * k_nopad[1] * k_nopad[2];
271 
272  // Read origin
273  fscanf(poten_fp, "origin %lf %lf %lf\n",
274  &origin.x, &origin.y, &origin.z);
275 
276  // Read delta (unit vectors)
277  // These are column vectors, so must fill gridfrcE tensor to reflect this
278  fscanf(poten_fp, "delta %lf %lf %lf\n", &e.xx, &e.yx, &e.zx);
279  fscanf(poten_fp, "delta %lf %lf %lf\n", &e.xy, &e.yy, &e.zy);
280  fscanf(poten_fp, "delta %lf %lf %lf\n", &e.xz, &e.yz, &e.zz);
281 
282  center = origin + e * 0.5
283  * Position(k_nopad[0]-1, k_nopad[1]-1, k_nopad[2]-1);
284 
285  fscanf(poten_fp, "object %*d class gridconnections counts %*lf %*lf %*lf\n");
286  fscanf(poten_fp, "object %*d class array type double rank 0 items %*d data follows\n");
287 
288  // Calculate inverse tensor
289  BigReal det;
290  det = e.xx*(e.yy*e.zz - e.yz*e.zy) - e.xy*(e.yx*e.zz - e.yz*e.zx)
291  + e.xz*(e.yx*e.zy - e.yy*e.zx);
292  inv.xx = (e.yy*e.zz - e.yz*e.zy)/det;
293  inv.xy = -(e.xy*e.zz - e.xz*e.zy)/det;
294  inv.xz = (e.xy*e.yz - e.xz*e.yy)/det;
295  inv.yx = -(e.yx*e.zz - e.yz*e.zx)/det;
296  inv.yy = (e.xx*e.zz - e.xz*e.zx)/det;
297  inv.yz = -(e.xx*e.yz - e.xz*e.yx)/det;
298  inv.zx = (e.yx*e.zy - e.yy*e.zx)/det;
299  inv.zy = -(e.xx*e.zy - e.xy*e.zx)/det;
300  inv.zz = (e.xx*e.yy - e.xy*e.yx)/det;
301 
302  DebugM(4, "origin = " << origin << "\n");
303  DebugM(4, "e = " << e << "\n");
304  DebugM(4, "inv = " << inv << "\n" << endi);
305 }
306 
307 void GridforceFullBaseGrid::readSubgridHierarchy(FILE *poten, int &totalGrids)
308 {
309  DebugM(4, "Beginning of readSubgridHierarchy, generation = " << generation << ", totalGrids = " << totalGrids << "\n" << endi);
310 
311  int elems, generation_in;
312 
314 
315  for (int i = 0; i < numSubgrids; i++) {
316  subgrids[i] = new GridforceFullSubGrid(this);
317  elems = fscanf(poten_fp, "# namdnugrid subgrid %d generation %d min %d %d %d max %d %d %d subgrids count %d\n",
318  &subgrids[i]->subgridIdx, &generation_in,
319  &subgrids[i]->pmin[0], &subgrids[i]->pmin[1], &subgrids[i]->pmin[2],
320  &subgrids[i]->pmax[0], &subgrids[i]->pmax[1], &subgrids[i]->pmax[2],
321  &subgrids[i]->numSubgrids);
322  if (elems < 9) {
323  char msg[256];
324  sprintf(msg, "Problem reading Gridforce potential file! (%d < 9)", elems);
325  NAMD_die(msg);
326  }
327 
328  totalGrids++;
329 
330  if (subgrids[i]->subgridIdx != (totalGrids - 1)) {
331  char msg[256];
332  sprintf(msg, "Problem reading Gridforce potential file! (%d != %d)", subgrids[i]->subgridIdx, totalGrids - 1);
333  NAMD_die(msg);
334  }
335  if (subgrids[i]->generation != generation_in) {
336  char msg[256];
337  sprintf(msg, "Problem reading Gridforce potential file! (%d != %d)", subgrids[i]->generation, generation_in);
338  NAMD_die(msg);
339  }
340 
341 // DebugM(3, "setting maingrid\n");
342 // subgrids[i]->maingrid->subgrids_flat[subgrids[i]->subgridIdx] = subgrids[i];
343 // DebugM(3, "reading subgrid hierarchy\n");
344 
345  subgrids[i]->readSubgridHierarchy(poten, totalGrids);
346  }
347 }
348 
349 
350 /*************************/
351 /* GRIDFORCEFULLMAINGRID */
352 /*************************/
353 
355 {
356  mygridnum = gridnum;
357  generation = 0;
358  subgrids_flat = NULL;
360 }
361 
362 
364 {
365  delete[] subgrids_flat;
366 }
367 
368 
370 {
371  DebugM(4, "Packing maingrid\n" << endi);
372 
373 // msg->put(3*sizeof(float), (char*)pad_p);
374 // msg->put(3*sizeof(float), (char*)pad_n);
375  msg->put(totalGrids);
376  msg->put(mygridnum);
377  msg->put(129*sizeof(char), (char*)filename);
378 
379  DebugM(3, "calling GridforceFullBaseGrid::pack\n" << endi);
380 
382 }
383 
384 
386 {
387  DebugM(4, "Unpacking maingrid\n" << endi);
388 
389 // msg->get(3*sizeof(float), (char*)pad_p);
390 // msg->get(3*sizeof(float), (char*)pad_n);
391  msg->get(totalGrids);
392  msg->get(mygridnum);
393  msg->get(129*sizeof(char), (char*)filename);
394 
396 
397  DebugM(4, "size = " << size << "\n");
398  DebugM(4, "numSubgrids = " << numSubgrids << "\n");
399  DebugM(4, "gapinv = " << gapinv[0] << " " << gapinv[2] << " " << gapinv[2] << " " << "\n");
400  DebugM(4, "generation = " << generation << "\n" << endi);
401  DebugM(4, "filename = " << filename << "\n" << endi);
402 
404 }
405 
406 
408 {
409  DebugM(4, "buildSubgridsFlat() called, totalGrids-1 = " << totalGrids-1 << "\n" << endi);
410  delete[] subgrids_flat;
412  for (int i = 0; i < numSubgrids; i++) {
413  DebugM(3, "adding to subgridsFlat\n" << endi);
415  DebugM(3, "success!\n" << endi);
416  }
417  for (int i = 0; i < totalGrids-1; i++) {
418  DebugM(4, "subgrids_flat[" << i << "]->numSubgrids = " << subgrids_flat[i]->numSubgrids << "\n" << endi);
419  }
420  for (int i = 0; i < numSubgrids; i++) {
421  DebugM(4, "subgrids[" << i << "]->numSubgrids = " << subgrids[i]->numSubgrids << "\n" << endi);
422  }
423 }
424 
425 
426 void GridforceFullMainGrid::initialize(char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams, int brd)
427 {
428  if (brd >= 0) {
429  border = brd;
430  } else {
432  }
433 
434  // FROM init1
435  //FILE *poten = Fopen(potfilename, "r");
436  poten_fp = Fopen(potfilename, "r");
437  if (!poten_fp) {
438  NAMD_die("Problem reading grid force potential file");
439  }
440 
441  // save file name so that grid can be re-read via Tcl
442  strcpy(filename, potfilename);
443 
444  // Read special comment fields and create subgrid objects
445  totalGrids = 1;
446  char line[256];
447  Bool flag = FALSE;
448  numSubgrids = 0;
449  float version;
450  long int poten_offset;
451  do {
452  poten_offset = ftell(poten_fp);
453  fgets(line, 256, poten_fp); // Read comment lines
454  //flag = sscanf(line, "# maingrid subgrids count %d\n", &numSubgrids);
455  flag = sscanf(line, "# namdnugrid version %f\n", &version);
456  } while (line[0] == '#' && !flag);
457 
458  if (flag) {
459  if (version != 1.0) {
460  NAMD_die("Unsupported version of non-uniform grid file format!");
461  }
462  fscanf(poten_fp, "# namdnugrid maingrid subgrids count %d\n", &numSubgrids);
465  } else {
466  fseek(poten_fp, poten_offset, SEEK_SET);
467  }
468 
469  // Read header
470  readHeader(simParams, mgridParams);
471 
472  factor = 1.0;
473  if (mgridParams->gridforceVolts)
474  {
475  factor /= 0.0434; // convert V -> kcal/mol*e
476  }
477  scale = mgridParams->gridforceScale;
478  checksize = mgridParams->gridforceCheckSize;
479 
480  // Allocate storage for potential and read it
481  float *grid_nopad = new float[size_nopad];
482 
483  float tmp2;
484  for (long int count = 0; count < size_nopad; count++) {
485  int err = fscanf(poten_fp, "%f", &tmp2);
486  if (err == EOF || err == 0) {
487  NAMD_die("Grid force potential file incorrectly formatted");
488  }
489  grid_nopad[count] = tmp2 * factor; // temporary, so just store flat
490  }
491  fscanf(poten_fp, "\n");
492 
493  // Shortcuts for accessing 1-D array with four indices
494  dk_nopad[0] = k_nopad[1] * k_nopad[2];
495  dk_nopad[1] = k_nopad[2];
496  dk_nopad[2] = 1;
497 
498  Vector Kvec[3];
499  Kvec[0] = e * Position(k_nopad[0]-1, 0, 0);
500  Kvec[1] = e * Position(0, k_nopad[1]-1, 0);
501  Kvec[2] = e * Position(0, 0, k_nopad[2]-1);
502  Vector Avec[3];
503  Avec[0] = simParams->lattice.a();
504  Avec[1] = simParams->lattice.b();
505  Avec[2] = simParams->lattice.c();
506 
507  // Decide whether we're wrapping
508  for (int i0 = 0; i0 < 3; i0++) {
509  if (mgridParams->gridforceCont[i0])
510  {
511  Bool found = FALSE;
512  for (int i1 = 0; i1 < 3; i1++) {
513  if (cross(Avec[i0].unit(), Kvec[i1].unit()).length() < 1e-4) {
514  found = TRUE;
515  cont[i1] = TRUE;
516  offset[i1] = mgridParams->gridforceVOffset[i0] * factor;
517  // want in grid-point units (normal = 1)
518  gap[i1] = (inv * (Avec[i0] - Kvec[i1])).length();
519  gapinv[i1] = 1.0/gap[i1];
520 
521  if (gap[i1] < 0) {
522  NAMD_die("Gridforce Grid overlap!");
523  }
524 
525  DebugM(4, "cont[" << i1 << "] = " << cont[i1] << "\n");
526  DebugM(4, "gap[" << i1 << "] = " << gap[i1] << "\n");
527  DebugM(4, "gapinv[" << i1 << "] = " << gapinv[i1] << "\n" << endi);
528  }
529  }
530 
531  if (!found) {
532  NAMD_die("No Gridforce unit vector found parallel to requested continuous grid direction!");
533  }
534  } else {
535  // check for grid overlap in non-wrapping dimensions
536  // occurs below
537  }
538  }
539 
540  // Figure out size of true grid (padded on non-periodic sides)
541  Vector delta = 0;
542  for (int i = 0; i < 3; i++) {
543  if (cont[i]) {
544  k[i] = k_nopad[i];
545  } else {
546  k[i] = k_nopad[i] + 2*border;
547  delta[i] -= border;
548  }
549  }
550  DebugM(4, "delta = " << e * delta << " (" << delta << ")\n" << endi);
551  origin += e * delta;
552 
553  // Check for grid overlap
554  if (!fits_lattice(simParams->lattice)) {
555  char errmsg[512];
556  if (checksize) {
557  sprintf(errmsg, "Warning: Periodic cell basis too small for Gridforce grid %d. Set gridforcechecksize off in configuration file to ignore.\n", mygridnum);
558  NAMD_die(errmsg);
559  }
560  }
561 
562  size = k[0] * k[1] * k[2];
563  dk[0] = k[1] * k[2];
564  dk[1] = k[2];
565  dk[2] = 1;
566 
567  DebugM(3, "size = " << size << ", size_nopad = " << size_nopad << "\n" << endi);
568 
569  delete[] grid;
570  grid = new float[size];
571 
572  n_sum[0] = n_sum[1] = n_sum[2] = 0;
573  p_sum[0] = p_sum[1] = p_sum[2] = 0;
574  for (int i0 = 0; i0 < k_nopad[0]; i0++) {
575  for (int i1 = 0; i1 < k_nopad[1]; i1++) {
576  for (int i2 = 0; i2 < k_nopad[2]; i2++) {
577  // Edges are special cases -- take force there to be
578  // zero for smooth transition across potential
579  // boundary
580 
581  long int ind_nopad = i0*dk_nopad[0] + i1*dk_nopad[1] + i2*dk_nopad[2];
582  int j0 = (cont[0]) ? i0 : i0 + border;
583  int j1 = (cont[1]) ? i1 : i1 + border;
584  int j2 = (cont[2]) ? i2 : i2 + border;
585  long int ind = j0*dk[0] + j1*dk[1] + j2*dk[2];
586 
587  if (i0 == 0) n_sum[0] += grid_nopad[ind_nopad];
588  else if (i0 == k_nopad[0]-1) p_sum[0] += grid_nopad[ind_nopad];
589  if (i1 == 0) n_sum[1] += grid_nopad[ind_nopad];
590  else if (i1 == k_nopad[1]-1) p_sum[1] += grid_nopad[ind_nopad];
591  if (i2 == 0) n_sum[2] += grid_nopad[ind_nopad];
592  else if (i2 == k_nopad[2]-1) p_sum[2] += grid_nopad[ind_nopad];
593 
594  //grid[ind] = grid_nopad[ind_nopad];
595  set_grid(j0, j1, j2, grid_nopad[ind_nopad]);
596  }
597  }
598  }
599 
600  const BigReal modThresh = 1.0;
601 
602  BigReal n_avg[3], p_avg[3];
603  int i0;
604  for (int i0 = 0; i0 < 3; i0++) {
605  int i1 = (i0 + 1) % 3;
606  int i2 = (i0 + 2) % 3;
607  n_avg[i0] = n_sum[i0] / (k_nopad[i1] * k_nopad[i2]);
608  p_avg[i0] = p_sum[i0] / (k_nopad[i1] * k_nopad[i2]);
609 
610  if (cont[i0] && fabs(offset[i0] - (p_avg[i0]-n_avg[i0])) > modThresh)
611  {
612  iout << iWARN << "GRID FORCE POTENTIAL DIFFERENCE IN K" << i0
613  << " DIRECTION IS "
614  << offset[i0] - (p_avg[i0]-n_avg[i0])
615  << " KCAL/MOL*E\n" << endi;
616  }
617  }
618 
619  Bool twoPadVals = (cont[0] + cont[1] + cont[2] == 2);
620  double padVal = 0.0;
621  long int weight = 0;
622  if (!twoPadVals) {
623  // Determine pad value (must average)
624  if (!cont[0]) {
625  padVal += p_sum[0] + n_sum[0];
626  weight += 2 * k_nopad[1] * k_nopad[2];
627  }
628  if (!cont[1]) {
629  padVal += p_sum[1] + n_sum[1];
630  weight += 2 * k_nopad[0] * k_nopad[2];
631  }
632  if (!cont[2]) {
633  padVal += p_sum[2] + n_sum[2];
634  weight += 2 * k_nopad[0] * k_nopad[1];
635  }
636  padVal /= weight;
637  }
638 
639  for (int i = 0; i < 3; i++) {
640  pad_n[i] = (cont[i]) ? 0.0 : (twoPadVals) ? n_avg[i] : padVal;
641  pad_p[i] = (cont[i]) ? 0.0 : (twoPadVals) ? p_avg[i] : padVal;
642  DebugM(4, "pad_n[" << i << "] = " << pad_n[i] << "\n");
643  DebugM(4, "pad_p[" << i << "] = " << pad_p[i] << "\n" << endi);
644  }
645 
646  if (cont[0] && cont[1] && cont[2]) {
647  // Nothing to do
648  return;
649  }
650 
651  // Now fill in rest of new grid
652  for (int i0 = 0; i0 < k[0]; i0++) {
653  for (int i1 = 0; i1 < k[1]; i1++) {
654  for (int i2 = 0; i2 < k[2]; i2++) {
655  if ( (cont[0] || (i0 >= border && i0 < k[0]-border))
656  && (cont[1] || (i1 >= border && i1 < k[1]-border))
657  && (cont[2] || i2 == border) )
658  {
659  i2 += k_nopad[2]-1;
660  continue;
661  }
662 
663  long int ind = i0*dk[0] + i1*dk[1] + i2*dk[2];
664 
665  Position pos = e * Position(i0, i1, i2);
666  int var[3] = {i0, i1, i2};
667 
668  for (int dir = 0; dir < 3; dir++) {
669  if (cont[dir])
670  continue;
671 
672  if (var[dir] < border) {
673  //grid[ind] = pad_n[dir];
674  set_grid(i0, i1, i2, pad_n[dir]);
675  } else if (var[dir] >= k[dir]-border) {
676  //grid[ind] = pad_p[dir];
677  set_grid(i0, i1, i2, pad_p[dir]);
678  }
679  }
680 
681 // DebugM(2, "grid[" << ind << "; " << i0 << ", " << i1
682 // << ", " << i2 << "] = " << get_grid(ind)
683 // << "\n" << endi);
684  }
685  }
686  }
687 
688  for (int i0 = 0; i0 < k[0]; i0++) {
689  for (int i1 = 0; i1 < k[1]; i1++) {
690  for (int i2 = 0; i2 < k[2]; i2++) {
691  DebugM(1, "grid[" << i0 << ", " << i1 << ", " << i2 << "] = " << get_grid(i0,i1,i2) << "\n" << endi);
692  }
693  }
694  }
695 
696  // Clean up
697  DebugM(3, "clean up\n" << endi);
698  delete[] grid_nopad;
699 
700  // Call initialize for each subgrid
701  for (int i = 0; i < numSubgrids; i++) {
702  subgrids[i]->poten_fp = poten_fp;
703  subgrids[i]->initialize(simParams, mgridParams);
704  }
705 
706  // close file pointer
707  fclose(poten_fp);
708 }
709 
710 
712 {
713  DebugM(4, "reinitializing grid\n" << endi);
714  initialize(filename, simParams, mgridParams);
715 }
716 
717 
718 long int GridforceFullMainGrid::get_all_gridvals(float **all_gridvals) const
719 {
720  // Creates a flat array of all grid values, including subgrids,
721  // and puts it in the value pointed to by the 'grids'
722  // argument. Returns the resulting array size. Caller is
723  // responsible for destroying the array via 'delete[]'
724 
725  DebugM(4, "get_all_gridvals called\n" << endi);
726 
727  long int sz = 0;
728  sz += size;
729  for (int i = 0; i < totalGrids-1; i++) {
730  sz += subgrids_flat[i]->size;
731  }
732  DebugM(4, "size = " << sz << "\n" << endi);
733 
734  float *grid_vals = new float[sz];
735  long int idx = 0;
736  for (long int i = 0; i < size; i++) {
737  grid_vals[idx++] = grid[i];
738  }
739  for (int j = 0; j < totalGrids-1; j++) {
740  for (long int i = 0; i < subgrids_flat[j]->size; i++) {
741  grid_vals[idx++] = subgrids_flat[j]->grid[i];
742  }
743  }
744  CmiAssert(idx == sz);
745 
746  *all_gridvals = grid_vals;
747 
748  DebugM(4, "get_all_gridvals finished\n" << endi);
749 
750  return sz;
751 }
752 
753 
754 void GridforceFullMainGrid::set_all_gridvals(float *all_gridvals, long int sz)
755 {
756  DebugM(4, "set_all_gridvals called\n" << endi);
757 
758  long int sz_calc = 0;
759  sz_calc += size;
760  for (int i = 0; i < totalGrids-1; i++) {
761  sz_calc += subgrids_flat[i]->size;
762  }
763  CmiAssert(sz == sz_calc);
764 
765  long int idx = 0;
766  for (long int i = 0; i < size; i++) {
767  DebugM(1, "all_gridvals[" << idx << "] = " << all_gridvals[idx] << "\n" << endi);
768  grid[i] = all_gridvals[idx++];
769  }
770  for (int j = 0; j < totalGrids-1; j++) {
771  for (long int i = 0; i < subgrids_flat[j]->size; i++) {
772  DebugM(1, "all_gridvals[" << idx << "] = " << all_gridvals[idx] << "\n" << endi);
773  subgrids_flat[j]->grid[i] = all_gridvals[idx++];
774  }
775  }
776  CmiAssert(idx == sz);
777 
778  DebugM(4, "set_all_gridvals finished\n" << endi);
779 }
780 
781 
782 void GridforceFullMainGrid::compute_b(float *b, int *inds, Vector gapscale) const
783 {
784  for (int i0 = 0; i0 < 8; i0++) {
785  int inds2[3];
786  int zero_derivs = FALSE;
787 
788  float voff = 0.0;
789  int bit = 1; // bit = 2^i1 in the below loop
790  for (int i1 = 0; i1 < 3; i1++) {
791  inds2[i1] = (inds[i1] + ((i0 & bit) ? 1 : 0)) % k[i1];
792 
793  // Deal with voltage offsets
794  if (cont[i1] && inds[i1] == (k[i1]-1) && inds2[i1] == 0) {
795  voff += offset[i1];
796  DebugM(3, "offset[" << i1 << "] = " << offset[i1] << "\n" << endi);
797  }
798 
799  bit <<= 1; // i.e. multiply by 2
800  }
801 
802  DebugM(1, "inds2 = " << inds2[0] << " " << inds2[1] << " " << inds2[2] << "\n" << endi);
803 
804  // NOTE: leaving everything in terms of unit cell coordinates for now,
805  // eventually will multiply by inv tensor when applying the force
806 
807  // First set variables 'dk_{hi,lo}' (glob notation). The 'hi'
808  // ('lo') variable in a given dimension is the number added (subtracted)
809  // to go up (down) one grid point in that dimension; both are normally
810  // just the corresponding 'dk[i]'. However, if we are sitting on a
811  // boundary and we are using a continuous grid, then we want to map the
812  // next point off the grid back around to the other side. e.g. point
813  // (k[0], i1, k) maps to point (0, i1, k), which would be
814  // accomplished by changing 'dk1_hi' to -(k[0]-1)*dk1.
815 
816  int d_hi[3] = {1, 1, 1};
817  int d_lo[3] = {1, 1, 1};
818  float voffs[3];
819  float dscales[3] = {0.5, 0.5, 0.5};
820  for (int i1 = 0; i1 < 3; i1++) {
821  if (inds2[i1] == 0) {
822  if (cont[i1]) {
823  d_lo[i1] = -(k[i1]-1);
824  voffs[i1] = offset[i1];
825  dscales[i1] = 1.0/(1.0 + gap[i1]) * 1.0/gapscale[i1];
826  }
827  else zero_derivs = TRUE;
828  }
829  else if (inds2[i1] == k[i1]-1) {
830  if (cont[i1]) {
831  d_hi[i1] = -(k[i1]-1);
832  voffs[i1] = offset[i1];
833  dscales[i1] = 1.0/(1.0 + gap[i1]) * 1.0/gapscale[i1];
834  }
835  else zero_derivs = TRUE;
836  }
837  else {
838  voffs[i1] = 0.0;
839  }
840  }
841 
842 // DebugM(2, "cont = " << cont[0] << " " << cont[1] << " " << cont[2] << "\n" << endi);
843 // DebugM(2, "zero_derivs = " << zero_derivs << "\n" << endi);
844 // DebugM(2, "d_hi = " << d_hi[0] << " " << d_hi[1] << " " << d_hi[2] << "\n" << endi);
845 // DebugM(2, "d_lo = " << d_lo[0] << " " << d_lo[1] << " " << d_lo[2] << "\n" << endi);
846  DebugM(1, "dscales = " << dscales[0] << " " << dscales[1] << " " << dscales[2] << "\n" << endi);
847  DebugM(1, "voffs = " << voffs[0] << " " << voffs[1] << " " << voffs[2] << "\n" << endi);
848 
849  // V
850  b[i0] = get_grid(inds2[0],inds2[1],inds2[2]) + voff;
851 
852  if (zero_derivs) {
853  DebugM(2, "zero_derivs\n" << endi);
854  b[8+i0] = 0.0;
855  b[16+i0] = 0.0;
856  b[24+i0] = 0.0;
857  b[32+i0] = 0.0;
858  b[40+i0] = 0.0;
859  b[48+i0] = 0.0;
860  b[56+i0] = 0.0;
861  } else {
862  b[8+i0] = dscales[0] * (get_grid_d(inds2[0]+d_hi[0],inds2[1],inds2[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1],inds2[2]) + voffs[0]); // dV/dx
863  b[16+i0] = dscales[1] * (get_grid_d(inds2[0],inds2[1]+d_hi[1],inds2[2]) - get_grid_d(inds2[0],inds2[1]-d_lo[1],inds2[2]) + voffs[1]); // dV/dy
864  b[24+i0] = dscales[2] * (get_grid_d(inds2[0],inds2[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0],inds2[1],inds2[2]-d_lo[2]) + voffs[2]); // dV/dz
865  b[32+i0] = dscales[0] * dscales[1]
866  * (get_grid_d(inds2[0]+d_hi[0],inds2[1]+d_hi[1],inds2[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1]+d_hi[1],inds2[2]) -
867  get_grid_d(inds2[0]+d_hi[0],inds2[1]-d_lo[1],inds2[2]) + get_grid_d(inds2[0]-d_lo[0],inds2[1]-d_lo[1],inds2[2])); // d2V/dxdy
868  b[40+i0] = dscales[0] * dscales[2]
869  * (get_grid_d(inds2[0]+d_hi[0],inds2[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1],inds2[2]+d_hi[2]) -
870  get_grid_d(inds2[0]+d_hi[0],inds2[1],inds2[2]-d_lo[2]) + get_grid_d(inds2[0]-d_lo[0],inds2[1],inds2[2]-d_lo[2])); // d2V/dxdz
871  b[48+i0] = dscales[1] * dscales[2]
872  * (get_grid_d(inds2[0],inds2[1]+d_hi[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0],inds2[1]-d_lo[1],inds2[2]+d_hi[2]) -
873  get_grid_d(inds2[0],inds2[1]+d_hi[1],inds2[2]-d_lo[2]) + get_grid_d(inds2[0],inds2[1]-d_lo[1],inds2[2]-d_lo[2])); // d2V/dydz
874 
875  b[56+i0] = dscales[0] * dscales[1] * dscales[2] // d3V/dxdydz
876  * (get_grid_d(inds2[0]+d_hi[0],inds2[1]+d_hi[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0]+d_hi[0],inds2[1]+d_hi[1],inds2[2]-d_lo[2]) -
877  get_grid_d(inds2[0]+d_hi[0],inds2[1]-d_lo[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1]+d_hi[1],inds2[2]+d_hi[2]) +
878  get_grid_d(inds2[0]+d_hi[0],inds2[1]-d_lo[1],inds2[2]-d_lo[2]) + get_grid_d(inds2[0]-d_lo[0],inds2[1]+d_hi[1],inds2[2]-d_lo[2]) +
879  get_grid_d(inds2[0]-d_lo[0],inds2[1]-d_lo[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1]-d_lo[1],inds2[2]-d_lo[2]));
880  }
881 
882  DebugM(1, "V = " << b[i0] << "\n");
883 
884  DebugM(1, "dV/dx = " << b[8+i0] << "\n");
885  DebugM(1, "dV/dy = " << b[16+i0] << "\n");
886  DebugM(1, "dV/dz = " << b[24+i0] << "\n");
887 
888  DebugM(1, "d2V/dxdy = " << b[32+i0] << "\n");
889  DebugM(1, "d2V/dxdz = " << b[40+i0] << "\n");
890  DebugM(1, "d2V/dydz = " << b[48+i0] << "\n");
891 
892  DebugM(1, "d3V/dxdydz = " << b[56+i0] << "\n" << endi);
893  }
894 }
895 
896 
897 /************************/
898 /* GRIDFORCEFULLSUBGRID */
899 /************************/
900 
902  parent = parent_in;
906  while (tmp->generation > 0) {
907  tmp = ((GridforceFullSubGrid *)tmp)->parent;
908  }
910  DebugM(4, "generation = " << generation << "\n" << endi);
911 }
912 
913 
915 {
916  int tmp;
917  char line[256];
918  long int poten_offset;
919 
920  // Skip 'attribute's
921  DebugM(3, "Skipping 'attribute' keywords...\n" << endi);
922  char str[256];
923  do {
924  poten_offset = ftell(poten_fp);
925  fscanf(poten_fp, "%s", str);
926  fgets(line, 256, poten_fp);
927  DebugM(4, "Read line " << str << " " << line << endi);
928  } while (strcmp(str, "attribute") == 0);
929  fseek(poten_fp, poten_offset, SEEK_SET);
930 
931  // Skip 'field' object
932  DebugM(3, "Skipping 'field' object\n" << endi);
933  fscanf(poten_fp, "object");
934  int n;
935  n = fscanf(poten_fp, "\"%[^\"]\" class field\n", str);
936  if (n == 0) {
937  n = fscanf(poten_fp, "%d class field\n", &tmp);
938  }
939 
940  if (n == 0) {
941  NAMD_die("Error reading gridforce grid! Could not find field object!\n");
942  }
943 
944  // Skip 'component's
945  DebugM(3, "Skipping 'component' keywords\n" << endi);
946  do {
947  poten_offset = ftell(poten_fp);
948  fscanf(poten_fp, "%s", str);
949  fgets(line, 256, poten_fp);
950  } while (strcmp(str, "component") == 0);
951  fseek(poten_fp, poten_offset, SEEK_SET);
952 
953  // Read header
954  readHeader(simParams, mgridParams);
955 
956  factor = 1.0;
957  if (mgridParams->gridforceVolts)
958  {
959  factor /= 0.0434; // convert V -> kcal/mol*e
960  }
961  scale = mgridParams->gridforceScale;
962  checksize = mgridParams->gridforceCheckSize;
963 
964  for (int i = 0; i < 3; i++) {
965  k[i] = k_nopad[i]; // subgrids aren't padded
966  }
967 
968  // Make sure that each subgrid dimension is an integral
969  // number of spanned supergrid cells. This is to ensure that no
970  // supergrid nodes land in the middle of a subgrid, because in
971  // this case forces could not be matched properly.
972  for (int i = 0; i < 3; i++) {
973  if ((k[i] - 1) % (pmax[i] - pmin[i] + 1) != 0) {
974  iout << (k[i] - 1) << " % " << (pmax[i] - pmin[i] + 1) << " != 0\n" << endi;
975  NAMD_die("Error reading gridforce grid! Subgrid dimensions must be an integral number spanned parent cells!");
976  }
977  }
978 
979  for (int i = 0; i < 3; i++) {
980  if (parent->cont[i]) {
981  cont[i] = (pmin[i] == 0 && pmax[i] == parent->k[i]-2) ? TRUE : FALSE;
982  DebugM(3, "pmin[" << i << "] = " << pmin[i] << " pmax[" << i << "] = " << pmax[i] << " parent->k[" << i << "] = " << parent->k[i] << " cont[" << i << "] = " << cont[i] << "\n" << endi);
983  } else {
984  cont[i] = false;
985  if (parent->generation == 0) {
986  // Add to pmin, pmax since parent got extra gridpoint layer(s) (maybe)
987  int brd = parent->get_border();
988  pmin[i] += brd;
989  pmax[i] += brd;
990  }
991  }
992  }
993 
994  DebugM(4, "pmin = " << pmin[0] << " " << pmin[1] << " " << pmin[2] << "\n");
995  DebugM(4, "pmax = " << pmax[0] << " " << pmax[1] << " " << pmax[2] << "\n" << endi);
996 
997  Vector origin2 = parent->origin + parent->e * Position(pmin[0], pmin[1], pmin[2]);
998  Vector escale, invscale;
999  for (int i = 0; i < 3; i++) {
1000  escale[i] = double(pmax[i] - pmin[i] + 1)/(k[i]-1);
1001  invscale[i] = 1.0/escale[i];
1002  if (cont[i]) { pmax[i]++; }
1003  }
1004  Tensor e2 = tensorMult(parent->e, Tensor::diagonal(escale));
1005 
1006  // Check that lattice parameters agree with min and max numbers
1007  // from subgrid hierarchy.
1008  double TOL2 = 1e-4; // Totally arbitrary
1009  if (pow(origin2.x-origin.x, 2) > TOL2 ||
1010  pow(origin2.y-origin.y, 2) > TOL2 ||
1011  pow(origin2.z-origin.z, 2) > TOL2 ||
1012  pow(e2.xx-e.xx, 2) > TOL2 ||
1013  pow(e2.xy-e.xy, 2) > TOL2 ||
1014  pow(e2.xz-e.xz, 2) > TOL2 ||
1015  pow(e2.yx-e.yx, 2) > TOL2 ||
1016  pow(e2.yy-e.yy, 2) > TOL2 ||
1017  pow(e2.yz-e.yz, 2) > TOL2 ||
1018  pow(e2.zx-e.zx, 2) > TOL2 ||
1019  pow(e2.zy-e.zy, 2) > TOL2 ||
1020  pow(e2.zz-e.zz, 2) > TOL2)
1021  {
1022  NAMD_die("Error reading gridforce grid! Subgrid lattice does not match!");
1023  }
1024 
1025  // Overwrite what was read from the header
1026  origin = origin2;
1027  e = e2;
1028 
1029  inv = tensorMult(Tensor::diagonal(invscale), parent->inv);
1030  for (int i = 0; i < 3; i++) {
1031  gap[i] = escale[i] * parent->gap[i];
1032  gapinv[i] = invscale[i] * parent->gapinv[i];
1033  offset[i] = parent->offset[i];
1034  }
1035  center = origin + e * 0.5 * Position(k[0], k[1], k[2]);
1036 
1037  DebugM(4, "origin = " << origin << "\n");
1038  DebugM(4, "e = " << e << "\n");
1039  DebugM(4, "inv = " << inv << "\n");
1040  DebugM(4, "gap = " << gap[0] << " " << gap[2] << " " << gap[2] << " " << "\n");
1041  DebugM(4, "gapinv = " << gapinv[0] << " " << gapinv[2] << " " << gapinv[2] << " " << "\n");
1042  DebugM(4, "numSubgrids = " << numSubgrids << "\n");
1043  DebugM(4, "k = " << k[0] << " " << k[1] << " " << k[2] << "\n");
1044  DebugM(4, "escale = " << escale << "\n");
1045  DebugM(4, "invscale = " << invscale << "\n" << endi);
1046 
1047  /*** Set members ***/
1048  size = k[0] * k[1] * k[2];
1049  dk[0] = k[1] * k[2];
1050  dk[1] = k[2];
1051  dk[2] = 1;
1052 
1053  scale_dV = Tensor::diagonal(escale);
1054  scale_d2V = Tensor::diagonal(Vector(escale.x*escale.y, escale.x*escale.z, escale.y*escale.z));
1055  scale_d3V = escale.x * escale.y * escale.z;
1056 
1057  DebugM(4, "scale_dV = " << scale_dV << "\n");
1058  DebugM(4, "scale_d2V = " << scale_d2V << "\n");
1059  DebugM(4, "scale_d3V = " << scale_d3V << "\n" << endi);
1060 
1061  // Allocate storage for potential and read it
1062  float *grid_tmp = new float[size];
1063 
1064  float tmp2;
1065  DebugM(3, "size_nopad = " << size_nopad << "\n");
1066  for (long int count = 0; count < size_nopad; count++) {
1067 // poten_offset = ftell(poten_fp);
1068 // fscanf(poten_fp, "%s", str);
1069 // fgets(line, 256, poten_fp);
1070 // DebugM(4, "Read line " << str << " " << line << endi);
1071 // fseek(poten_fp, poten_offset, SEEK_SET);
1072 
1073  int err = fscanf(poten_fp, "%f", &tmp2);
1074  if (err == EOF || err == 0) {
1075  NAMD_die("Grid force potential file incorrectly formatted");
1076  }
1077  grid_tmp[count] = tmp2 * factor;
1078  }
1079  fscanf(poten_fp, "\n");
1080 
1081  // Set real grid
1082  DebugM(3, "allocating grid\n" << endi);
1083  delete[] grid;
1084  grid = new float[size];
1085  for (int i0 = 0; i0 < k_nopad[0]; i0++) {
1086  for (int i1 = 0; i1 < k_nopad[1]; i1++) {
1087  for (int i2 = 0; i2 < k_nopad[2]; i2++) {
1088  long int ind = i0*dk[0] + i1*dk[1] + i2*dk[2];
1089  set_grid(i0, i1, i2, grid_tmp[ind]);
1090  }
1091  }
1092  }
1093 
1094  for (int i0 = 0; i0 < k[0]; i0++) {
1095  for (int i1 = 0; i1 < k[1]; i1++) {
1096  for (int i2 = 0; i2 < k[2]; i2++) {
1097  DebugM(1, "grid[" << i0 << ", " << i1 << ", " << i2 << "] = " << get_grid(i0,i1,i2) << "\n" << endi);
1098  }
1099  }
1100  }
1101 
1102  // Clean up
1103  delete[] grid_tmp;
1104 
1105  // Call initialize for each subgrid
1106  for (int i = 0; i < numSubgrids; i++) {
1107  subgrids[i]->initialize(simParams, mgridParams);
1108  }
1109 }
1110 
1111 
1113 {
1114  DebugM(4, "Packing subgrid\n" << endi);
1115 
1116  msg->put(sizeof(Tensor), (char*)&scale_dV);
1117  msg->put(sizeof(Tensor), (char*)&scale_d2V);
1118  msg->put(sizeof(float), (char*)&scale_d3V);
1119 
1120  msg->put(3*sizeof(int), (char*)pmin);
1121  msg->put(3*sizeof(int), (char*)pmax);
1122  msg->put(subgridIdx);
1123 
1124  DebugM(3, "calling GridforceFullBaseGrid::pack\n" << endi);
1125 
1127 }
1128 
1129 
1131 {
1132  DebugM(4, "Unpacking subgrid\n" << endi);
1133 
1134  msg->get(sizeof(Tensor), (char*)&scale_dV);
1135  msg->get(sizeof(Tensor), (char*)&scale_d2V);
1136  msg->get(sizeof(float), (char*)&scale_d3V);
1137 
1138  msg->get(3*sizeof(int), (char*)pmin);
1139  msg->get(3*sizeof(int), (char*)pmax);
1140  msg->get(subgridIdx);
1141 
1143 
1144  DebugM(4, "size = " << size << "\n");
1145  DebugM(4, "numSubgrids = " << numSubgrids << "\n");
1146  DebugM(4, "gapinv = " << gapinv[0] << " " << gapinv[2] << " " << gapinv[2] << " " << "\n");
1147  DebugM(4, "generation = " << generation << "\n" << endi);
1148 }
1149 
1150 
1152 {
1153  DebugM(4, "addToSubgridsFlat() called, subgridIdx = " << subgridIdx << ", maingrid->numSubgrids = " << maingrid->numSubgrids << "\n" << endi);
1154  maingrid->subgrids_flat[subgridIdx-1] = this;
1155  for (int i = 0; i < numSubgrids; i++) {
1157  }
1158 }
1159 
1160 
1161 void GridforceFullSubGrid::compute_b(float *b, int *inds, Vector gapscale) const
1162 {
1163  for (int i0 = 0; i0 < 8; i0++) {
1164  int inds2[3];
1165 
1166  float voff = 0.0;
1167  int bit = 1; // bit = 2^i1 in the below loop
1168  for (int i1 = 0; i1 < 3; i1++) {
1169  inds2[i1] = (inds[i1] + ((i0 & bit) ? 1 : 0)) % k[i1];
1170 
1171  // Deal with voltage offsets
1172  if (cont[i1] && inds[i1] == (k[i1]-1) && inds2[i1] == 0) {
1173  voff += offset[i1];
1174  DebugM(3, "offset[" << i1 << "] = " << offset[i1] << "\n" << endi);
1175  }
1176 
1177  bit <<= 1; // i.e. multiply by 2
1178  }
1179 
1180  DebugM(3, "inds2 = " << inds2[0] << " " << inds2[1] << " " << inds2[2] << "\n" << endi);
1181 
1182  int d_hi[3] = {1, 1, 1};
1183  int d_lo[3] = {1, 1, 1};
1184  float voffs[3];
1185  float dscales[3] = {0.5, 0.5, 0.5};
1186  for (int i1 = 0; i1 < 3; i1++) {
1187  if (inds2[i1] == 0 && cont[i1]) {
1188  d_lo[i1] = -(k[i1]-1);
1189  voffs[i1] = offset[i1];
1190  dscales[i1] = 1.0/(1.0 + gap[i1]) * 1.0/gapscale[i1];
1191  }
1192  else if (inds2[i1] == k[i1]-1 && cont[i1]) {
1193  d_hi[i1] = -(k[i1]-1);
1194  voffs[i1] = offset[i1];
1195  dscales[i1] = 1.0/(1.0 + gap[i1]) * 1.0/gapscale[i1];
1196  }
1197  else {
1198  voffs[i1] = 0.0;
1199  }
1200  }
1201 
1202  bool edge = false;
1203  for (int i1 = 0; i1 < 3; i1++) {
1204  if (!cont[i1] && (inds2[i1] == 0 || inds2[i1] == k[i1]-1)) {
1205  edge = true;
1206  }
1207  }
1208 
1209  if (inds2[2] == 0) {
1210 // DebugM(3, "cont = " << cont[0] << " " << cont[1] << " " << cont[2] << " d_hi = " << d_hi[0] << " " << d_hi[1] << " " << d_hi[2] << " d_lo = " << d_lo[0] << " " << d_lo[1] << " " << d_lo[2] << " dscales = " << dscales[0] << " " << dscales[1] << " " << dscales[2] << "\n" << endi);
1211  DebugM(3, "cont = " << cont[0] << " " << cont[1] << " " << cont[2] << "\n" << endi);
1212  }
1213 
1214  if (edge) {
1215  DebugM(2, "Edge!\n" << endi);
1216 
1217  // Must get derivatives from parent
1218  Position pos = e * Vector(inds2[0], inds2[1], inds2[2]) + origin; // Gridpoint position in realspace
1219  Vector g = parent->inv * (pos - parent->origin); // Gridpoint position in parent's gridspace
1220  Vector dg;
1221  int inds3[3];
1222 
1223  DebugM(2, "g = " << g << "\n" << endi);
1224 
1225  for (int i = 0; i < 3; i++) {
1226  inds3[i] = (int)floor(g[i]);
1227  dg[i] = g[i] - inds3[i];
1228  }
1229 
1230  float x[4], y[4], z[4];
1231  x[0] = 1; y[0] = 1; z[0] = 1;
1232  for (int j = 1; j < 4; j++) {
1233  x[j] = x[j-1] * dg.x;
1234  y[j] = y[j-1] * dg.y;
1235  z[j] = z[j-1] * dg.z;
1236  DebugM(1, "x[" << j << "] = " << x[j] << "\n");
1237  DebugM(1, "y[" << j << "] = " << y[j] << "\n");
1238  DebugM(1, "z[" << j << "] = " << z[j] << "\n" << endi);
1239  }
1240 
1241  // Compute parent matrices
1242  float b_parent[64];
1243  parent->compute_b(b_parent, inds3, gapscale);
1244 
1245  float a_parent[64];
1246  parent->compute_a(a_parent, b_parent);
1247 
1248  // Compute parent derivatives
1249  float V = parent->compute_V(a_parent, x, y, z);
1250  Vector dV = scale_dV * parent->compute_dV(a_parent, x, y, z);
1251  Vector d2V = scale_d2V * parent->compute_d2V(a_parent, x, y, z);
1252  float d3V = scale_d3V * parent->compute_d3V(a_parent, x, y, z);
1253 
1254  b[i0] = V;
1255  b[8+i0] = dV[0];
1256  b[16+i0] = dV[1];
1257  b[24+i0] = dV[2];
1258  b[32+i0] = d2V[0];
1259  b[40+i0] = d2V[1];
1260  b[48+i0] = d2V[2];
1261  b[56+i0] = d3V;
1262  } else {
1263  b[i0] = get_grid(inds2[0],inds2[1],inds2[2]) + voff; // V
1264 
1265  b[8+i0] = dscales[0] * (get_grid_d(inds2[0]+d_hi[0],inds2[1],inds2[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1],inds2[2]) + voffs[0]); // dV/dx
1266  b[16+i0] = dscales[1] * (get_grid_d(inds2[0],inds2[1]+d_hi[1],inds2[2]) - get_grid_d(inds2[0],inds2[1]-d_lo[1],inds2[2]) + voffs[1]); // dV/dy
1267  b[24+i0] = dscales[2] * (get_grid_d(inds2[0],inds2[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0],inds2[1],inds2[2]-d_lo[2]) + voffs[2]); // dV/dz
1268  b[32+i0] = dscales[0] * dscales[1]
1269  * (get_grid_d(inds2[0]+d_hi[0],inds2[1]+d_hi[1],inds2[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1]+d_hi[1],inds2[2]) -
1270  get_grid_d(inds2[0]+d_hi[0],inds2[1]-d_lo[1],inds2[2]) + get_grid_d(inds2[0]-d_lo[0],inds2[1]-d_lo[1],inds2[2])); // d2V/dxdy
1271  b[40+i0] = dscales[0] * dscales[2]
1272  * (get_grid_d(inds2[0]+d_hi[0],inds2[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1],inds2[2]+d_hi[2]) -
1273  get_grid_d(inds2[0]+d_hi[0],inds2[1],inds2[2]-d_lo[2]) + get_grid_d(inds2[0]-d_lo[0],inds2[1],inds2[2]-d_lo[2])); // d2V/dxdz
1274  b[48+i0] = dscales[1] * dscales[2]
1275  * (get_grid_d(inds2[0],inds2[1]+d_hi[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0],inds2[1]-d_lo[1],inds2[2]+d_hi[2]) -
1276  get_grid_d(inds2[0],inds2[1]+d_hi[1],inds2[2]-d_lo[2]) + get_grid_d(inds2[0],inds2[1]-d_lo[1],inds2[2]-d_lo[2])); // d2V/dydz
1277 
1278  b[56+i0] = dscales[0] * dscales[1] * dscales[2] // d3V/dxdydz
1279  * (get_grid_d(inds2[0]+d_hi[0],inds2[1]+d_hi[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0]+d_hi[0],inds2[1]+d_hi[1],inds2[2]-d_lo[2]) -
1280  get_grid_d(inds2[0]+d_hi[0],inds2[1]-d_lo[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1]+d_hi[1],inds2[2]+d_hi[2]) +
1281  get_grid_d(inds2[0]+d_hi[0],inds2[1]-d_lo[1],inds2[2]-d_lo[2]) + get_grid_d(inds2[0]-d_lo[0],inds2[1]+d_hi[1],inds2[2]-d_lo[2]) +
1282  get_grid_d(inds2[0]-d_lo[0],inds2[1]-d_lo[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1]-d_lo[1],inds2[2]-d_lo[2]));
1283  }
1284 
1285  if (inds2[0] == 1 && inds2[1] == 1 && inds2[2] == 0) {
1286  DebugM(1, "Sub V = " << b[i0] << "\n");
1287 
1288  DebugM(1, "Sub dV/dx = " << b[8+i0] << "\n");
1289  DebugM(1, "Sub dV/dy = " << b[16+i0] << "\n");
1290  DebugM(1, "Sub dV/dz = " << b[24+i0] << "\n");
1291 
1292  DebugM(1, "Sub d2V/dxdy = " << b[32+i0] << "\n");
1293  DebugM(1, "Sub d2V/dxdz = " << b[40+i0] << "\n");
1294  DebugM(1, "Sub d2V/dydz = " << b[48+i0] << "\n");
1295 
1296  DebugM(1, "Sub d3V/dxdydz = " << b[56+i0] << "\n" << endi);
1297  }
1298  }
1299 }
1300 
1301 
1302 /*********************/
1303 /* GRIDFORCELITEGRID */
1304 /*********************/
1305 
1307 {
1308  mygridnum = gridnum;
1309  grid = NULL;
1311 }
1312 
1313 
1315 {
1316  delete[] grid;
1317 }
1318 
1319 
1321 {
1322  // cheat and use GridforceFullMainGrid to read the file
1324  tmp_grid->initialize(potfilename, simParams, mgridParams, 1);
1325 
1326  if (tmp_grid->get_total_grids() != 1) {
1327  NAMD_die("Cannot use gridforcelite option with multi-resolution grid!");
1328  }
1329 
1330  // save file name so that grid can be re-read via Tcl
1331  strcpy(filename, potfilename);
1332 
1333  // copy parameters
1334  k[0] = tmp_grid->get_k0();
1335  k[1] = tmp_grid->get_k1();
1336  k[2] = tmp_grid->get_k2();
1337  k[3] = 4; // for V, dV/dx, dV/dy, dV/dz grids
1338  origin = tmp_grid->get_origin();
1339  center = tmp_grid->get_center();
1340  e = tmp_grid->get_e();
1341  inv = tmp_grid->get_inv();
1342  scale = tmp_grid->get_scale();
1343 
1344  // calculate rest of parameters
1345  size = k[0] * k[1] * k[2] * k[3];
1346  dk[0] = k[1] * k[2] * k[3];
1347  dk[1] = k[2] * k[3];
1348  dk[2] = k[3];
1349  dk[3] = 1;
1350 
1351  // copy the potential grid
1352  delete[] grid;
1353  grid = new float[size];
1354  for (int i0 = 0; i0 < k[0]; i0++) {
1355  for (int i1 = 0; i1 < k[1]; i1++) {
1356  for (int i2 = 0; i2 < k[2]; i2++) {
1357  float V = tmp_grid->get_grid(i0, i1, i2);
1358  set_grid(i0, i1, i2, 0, V);
1359  DebugM(1, "V[" << i0 << "," << i1 << "," << i2 << "] = " << get_grid(i0, i1, i2, 0) << "(" << V << ")\n" << endi);
1360  }
1361  }
1362  }
1363 
1364  delete tmp_grid;
1365 
1367 }
1368 
1369 
1371 {
1372  // calculate derivative grids
1373  // separate loop so all V values have been set already
1374  for (int i0 = 0; i0 < k[0]; i0++) {
1375  for (int i1 = 0; i1 < k[1]; i1++) {
1376  for (int i2 = 0; i2 < k[2]; i2++) {
1377  float dx, dy, dz;
1378  if (i0 == 0 || i0 == k[0]-1 || i1 == 0 || i1 == k[1]-1 || i2 == 0 || i2 == k[2]-1) {
1379  // on edge, set ALL derivatives to zero (make up for lack of padding)
1380  dx = 0;
1381  dy = 0;
1382  dz = 0;
1383  } else {
1384  dx = 0.5 * (get_grid_d(i0+1,i1,i2,0) - get_grid_d(i0-1,i1,i2,0));
1385  dy = 0.5 * (get_grid_d(i0,i1+1,i2,0) - get_grid_d(i0,i1-1,i2,0));
1386  dz = 0.5 * (get_grid_d(i0,i1,i2+1,0) - get_grid_d(i0,i1,i2-1,0));
1387  }
1388  set_grid(i0, i1, i2, 1, dx);
1389  set_grid(i0, i1, i2, 2, dy);
1390  set_grid(i0, i1, i2, 3, dz);
1391  DebugM(1, "dx[" << i0 << "," << i1 << "," << i2 << "] = " << get_grid(i0, i1, i2, 1) << "(" << dx << ")\n" << endi);
1392  DebugM(1, "dy[" << i0 << "," << i1 << "," << i2 << "] = " << get_grid(i0, i1, i2, 2) << "(" << dy << ")\n" << endi);
1393  DebugM(1, "dz[" << i0 << "," << i1 << "," << i2 << "] = " << get_grid(i0, i1, i2, 3) << "(" << dz << ")\n" << endi);
1394  }
1395  }
1396  }
1397 }
1398 
1399 
1401 {
1402  initialize(filename, simParams, mgridParams);
1403 }
1404 
1405 
1407 {
1408  msg->put(4*sizeof(int), (char*)k);
1409  msg->put(size);
1410  msg->put(4*sizeof(long int), (char*)dk);
1411 
1412  msg->put(sizeof(Vector), (char*)&origin);
1413  msg->put(sizeof(Vector), (char*)&center);
1414  msg->put(sizeof(Tensor), (char*)&e);
1415  msg->put(sizeof(Tensor), (char*)&inv);
1416  msg->put(sizeof(Vector), (char*)&scale);
1417  msg->put(sizeof(Bool), (char*)&checksize);
1418 
1419  msg->put(129*sizeof(char), (char*)filename);
1420 
1421  msg->put(size*sizeof(float), (char*)grid);
1422 }
1423 
1424 
1426 {
1427  delete[] grid;
1428  grid = NULL;
1429 
1430  msg->get(4*sizeof(int), (char*)k);
1431  msg->get(size);
1432  msg->get(4*sizeof(long int), (char*)dk);
1433 
1434  msg->get(sizeof(Vector), (char*)&origin);
1435  msg->get(sizeof(Vector), (char*)&center);
1436  msg->get(sizeof(Tensor), (char*)&e);
1437  msg->get(sizeof(Tensor), (char*)&inv);
1438  msg->get(sizeof(Vector), (char*)&scale);
1439  msg->get(sizeof(Bool), (char*)&checksize);
1440 
1441  msg->get(129*sizeof(char), (char*)filename);
1442 
1443  if (size) {
1444  grid = new float[size];
1445  msg->get(size*sizeof(float), (char*)grid);
1446  }
1447 }
1448 
1449 
1450 long int GridforceLiteGrid::get_all_gridvals(float** all_gridvals) const
1451 {
1452  // Creates a flat array of all grid values and puts it in the
1453  // value pointed to by the 'all_gridvals' argument. Returns the
1454  // resulting array size. Caller is responsible for destroying the
1455  // array via 'delete[]'
1456 
1457  DebugM(4, "GridforceLiteGrid::get_all_gridvals called\n" << endi);
1458 
1459  long int sz = size;
1460  DebugM(4, "size = " << sz << "\n" << endi);
1461 
1462  float *grid_vals = new float[sz];
1463  long int idx = 0;
1464  for (long int i = 0; i < size; i++) {
1465  grid_vals[idx++] = grid[i];
1466  }
1467  CmiAssert(idx == sz);
1468 
1469  *all_gridvals = grid_vals;
1470 
1471  DebugM(4, "GridforceLiteGrid::get_all_gridvals finished\n" << endi);
1472 
1473  return sz;
1474 }
1475 
1476 
1477 void GridforceLiteGrid::set_all_gridvals(float* all_gridvals, long int sz)
1478 {
1479  DebugM(4, "GridforceLiteGrid::set_all_gridvals called\n" << endi);
1480 
1481  long int sz_calc = size;
1482  CmiAssert(sz == sz_calc);
1483 
1484  long int idx = 0;
1485  for (long int i = 0; i < size; i++) {
1486  grid[i] = all_gridvals[idx++];
1487  }
1488  CmiAssert(idx == sz);
1489 
1490  //compute_derivative_grids(); // not needed if we're sending all 4 grids
1491 
1492  DebugM(4, "GridforceLiteGrid::set_all_gridvals finished\n" << endi);
1493 }
void reinitialize(SimParameters *simParams, MGridforceParams *mgridParams)
void compute_b(float *b, int *inds, Vector gapscale) const
void pack(MOStream *msg) const
BigReal zy
Definition: Tensor.h:19
float compute_V(float *a, float *x, float *y, float *z) const
BigReal xz
Definition: Tensor.h:17
double get_grid_d(int i0, int i1, int i2, int i3) const
friend class GridforceFullSubGrid
Definition: GridForceGrid.h:84
virtual ~GridforceGrid()
Definition: GridForceGrid.C:48
static Tensor diagonal(const Vector &v1)
Definition: Tensor.h:37
float get_grid(int i0, int i1, int i2) const
static GridforceGrid * new_grid(int gridnum, char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams)
Definition: GridForceGrid.C:34
Definition: Vector.h:64
virtual void compute_b(float *b, int *inds, Vector gapscale) const =0
void initialize(char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams, int border)
Vector compute_dV(float *a, float *x, float *y, float *z) const
static void pack_grid(GridforceGrid *grid, MOStream *msg)
Definition: GridForceGrid.C:50
virtual Position get_origin(void) const =0
#define DebugM(x, y)
Definition: Debug.h:59
__device__ __forceinline__ float3 cross(const float3 v1, const float3 v2)
Tensor get_e(void) const
BigReal z
Definition: Vector.h:66
#define FALSE
Definition: common.h:116
BigReal yz
Definition: Tensor.h:18
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:108
MIStream * get(char &data)
Definition: MStream.h:29
virtual ~GridforceLiteGrid()
void pack(MOStream *msg) const
#define iout
Definition: InfoStream.h:87
virtual int get_k2(void) const =0
void compute_a(float *a, float *b) const
int get_k1(void) const
int get_k0(void) const
float get_grid(int i0, int i1, int i2, int i3) const
void compute_b(float *b, int *inds, Vector gapscale) const
float compute_d3V(float *a, float *x, float *y, float *z) const
virtual void pack(MOStream *msg) const
GridforceFullMainGrid * maingrid
void readHeader(SimParameters *simParams, MGridforceParams *mgridParams)
void unpack(MIStream *msg)
void compute_derivative_grids(void)
void set_grid(int i0, int i1, int i2, int i3, float V)
GridforceFullSubGrid ** subgrids
GridforceFullBaseGrid * parent
bool fits_lattice(const Lattice &lattice)
Definition: GridForceGrid.C:84
void set_grid(int i0, int i1, int i2, float V)
virtual void pack(MOStream *msg) const =0
void NAMD_bug(const char *err_msg)
Definition: common.C:123
BigReal yx
Definition: Tensor.h:18
GridforceFullSubGrid ** subgrids_flat
GridforceFullMainGrid(int gridnum)
virtual Position get_center(void) const =0
void reinitialize(SimParameters *simParams, MGridforceParams *mgridParams)
virtual ~GridforceFullBaseGrid()
gridSize z
int Bool
Definition: common.h:131
void set_all_gridvals(float *all_gridvals, long int sz)
FILE * Fopen(const char *filename, const char *mode)
Definition: common.C:265
Vector Position
Definition: NamdTypes.h:18
GridforceFullSubGrid(GridforceFullBaseGrid *parent_in)
virtual void initialize(char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams)=0
BigReal x
Definition: Vector.h:66
void readSubgridHierarchy(FILE *poten, int &totalGrids)
GridforceLiteGrid(int gridnum)
void addToSubgridsFlat(void)
void NAMD_die(const char *err_msg)
Definition: common.C:83
virtual int get_border(void) const =0
virtual void unpack(MIStream *msg)
int get_k2(void) const
long int get_all_gridvals(float **all_gridvals) const
void set_all_gridvals(float *all_gridvals, long int sz)
BigReal xx
Definition: Tensor.h:17
Position get_origin(void) const
BigReal zz
Definition: Tensor.h:19
GridforceGridType type
Definition: GridForceGrid.h:74
#define simParams
Definition: Output.C:127
void initialize(char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams)
int get_total_grids(void) const
long int get_all_gridvals(float **all_gridvals) const
static GridforceGrid * unpack_grid(int gridnum, MIStream *msg)
Definition: GridForceGrid.C:60
void buildSubgridsFlat(void)
Definition: Tensor.h:15
Tensor get_inv(void) const
BigReal xy
Definition: Tensor.h:17
Tensor tensorMult(const Tensor &t1, const Tensor &t2)
virtual int get_k0(void) const =0
BigReal y
Definition: Vector.h:66
Position get_corner(int idx)
Vector b() const
Definition: Lattice.h:253
virtual int get_k1(void) const =0
Position wrap_position(const Position &pos, const Lattice &lattice)
virtual void unpack(MIStream *msg)=0
virtual Tensor get_e(void) const =0
BigReal yy
Definition: Tensor.h:18
Vector get_scale(void) const
virtual ~GridforceFullMainGrid()
gridSize y
MOStream * put(char data)
Definition: MStream.h:112
GridforceGridType get_grid_type(void)
Definition: GridForceGrid.h:65
double get_grid_d(int i0, int i1, int i2) const
static const int default_border
infostream & endi(infostream &s)
Definition: InfoStream.C:38
gridSize x
void initialize(SimParameters *simParams, MGridforceParams *mgridParams)
void pack(MOStream *msg) const
BigReal zx
Definition: Tensor.h:19
Vector a() const
Definition: Lattice.h:252
#define TRUE
Definition: common.h:117
Vector compute_d2V(float *a, float *x, float *y, float *z) const
Vector c() const
Definition: Lattice.h:254
Position get_center(void) const
double BigReal
Definition: common.h:112
void unpack(MIStream *msg)
void unpack(MIStream *msg)