NAMD
ComputeLCPO.C
Go to the documentation of this file.
1 
6 #include "ComputeLCPO.h"
7 #include "ReductionMgr.h"
8 #include "PatchMap.h"
9 #include "ComputeMgr.h"
10 #include "Molecule.h"
11 #include "Node.h"
12 #include "SimParameters.h"
13 #include "Debug.h"
14 #include "WorkDistrib.decl.h"
15 #include "Node.h"
16 #include "ComputeLCPO.h"
17 #include "Priorities.h"
18 #include "PatchMap.inl"
19 #include "Patch.h"
20 #include "ComputeMap.h"
21 #include "LdbCoordinator.h"
22 #include "common.h"
23 #include "time.h"
24 
25 //#define COUNT_FLOPS
26 
27 #ifdef COUNT_FLOPS
28 #define FLOPS(X) flops += X;
29 #else
30 #define FLOPS(X)
31 #endif
32 
33 //#define MIN_DEBUG_LEVEL 4
34 // #define DEBUGM
35 
37  ComputeNonbondedWorkArrays* _workArrays,
38  int minPartition, int maxPartition, int numPartitions, int numPatches)
39  : Compute(c), workArrays(_workArrays),
40  minPart(minPartition), maxPart(maxPartition),
41  strideIg(numPartitions), numParts(numPartitions),
42  maxAtomRadius(1.9+1.4)
43  {
44 
46 
47  setNumPatches(8);
49  surfTen = simParams->surface_tension;
50 
51  for (int i=0; i<getNumPatches(); i++) {
52  patchID[i] = p[i];
53  trans[i] = t[i];
54  patch[i] = NULL;
55  positionBox[i] = NULL;
56  forceBox[i] = NULL;
57  lcpoTypeBox[i] = NULL;
58  } // for all patches
59 } // constructor
60 
62  DebugM(4, "~ComputeLCPO("<<cid<<") numAtoms("<<patchID[0]<<") = "
63  << numAtoms[0]
64  << " numAtoms("<<patchID[1]<<") = " << numAtoms[1] << "\n" );
65  DebugM(4, "~ComputeLCPO("<<cid<<") addr("<<patchID[0]<<") = "
66  << PatchMap::Object()->patch(patchID[0]) << " addr("<<patchID[1]<<") = "
67  << PatchMap::Object()->patch(patchID[1]) << "\n");
68 
69  for (int i=0; i<getNumPatches(); i++) {
70  if (positionBox[i] != NULL) {
72  &positionBox[i]);
73  }
74  if (forceBox[i] != NULL) {
76  &forceBox[i]);
77  }
78  if (lcpoTypeBox[i] != NULL) {
80  &lcpoTypeBox[i]);
81  }
82  }
83  delete reduction;
84 } // destructor
85 
88  // How can we tell if BoxOwner has packed up and left? Need a mechanism
89  // to handle this or do we assume the Boxes have been dumped?
90  PatchMap *patchMap = PatchMap::Object();
91  //Check out Boxes
92  for (int i=0; i<8; i++) {
93  //invalid patch so don't even checkout boxes
94  if (positionBox[i] == NULL) { // We have yet to get boxes
95  patch[i] = PatchMap::Object()->patch(patchID[i]);
96  if (!(patch[i] = PatchMap::Object()->patch(patchID[i]))) {
97  DebugM(5,"invalid patch(" << patchID[i]
98  << ") pointer!\n");
99  }
101  forceBox[i] = patch[i]->registerForceDeposit(this);
103  // will need to open a box full of lcpo parameters
104  }
105  numAtoms[i] = patch[i]->getNumAtoms();
106  } // for all patches
107 
108  // set priority
110 
111  //get bounds of inner rectangular prism in octet
112  bounds[0][0] = 0.5*(patchMap->min_a(patchID[0])+patchMap->max_a(patchID[0]));
113  bounds[1][0] = 0.5*(patchMap->min_b(patchID[0])+patchMap->max_b(patchID[0]));
114  bounds[2][0] = 0.5*(patchMap->min_c(patchID[0])+patchMap->max_c(patchID[0]));
115  bounds[0][1] = 0.5*(patchMap->min_a(patchID[7])+patchMap->max_a(patchID[7]));
116  bounds[1][1] = 0.5*(patchMap->min_b(patchID[7])+patchMap->max_b(patchID[7]));
117  bounds[2][1] = 0.5*(patchMap->min_c(patchID[7])+patchMap->max_c(patchID[7]));
118 
119  //if only 1 patch in a dimenion, invalidate those patches
120  int gsa = patchMap->gridsize_a();
121  int gsb = patchMap->gridsize_b();
122  int gsc = patchMap->gridsize_c();
123  invalidPatch[0] = 0;
124  invalidPatch[1] = 0;
125  invalidPatch[2] = 0;
126  invalidPatch[3] = 0;
127  invalidPatch[4] = 0;
128  invalidPatch[5] = 0;
129  invalidPatch[6] = 0;
130  invalidPatch[7] = 0;
131 
132  if (gsa==1) {
133  //CkPrintf("ONLY 1 PATCH in A DIMENSION!\n");
134  invalidPatch[1] = 1;
135  invalidPatch[3] = 1;
136  invalidPatch[5] = 1;
137  invalidPatch[7] = 1;
138  }
139  if (gsb==1) {
140  //CkPrintf("ONLY 1 PATCH in B DIMENSION!\n");
141  invalidPatch[2] = 1;
142  invalidPatch[3] = 1;
143  invalidPatch[6] = 1;
144  invalidPatch[7] = 1;
145  }
146  if (gsc==1) {
147  //CkPrintf("ONLY 1 PATCH in C DIMENSION!\n");
148  invalidPatch[4] = 1;
149  invalidPatch[5] = 1;
150  invalidPatch[6] = 1;
151  invalidPatch[7] = 1;
152  }
153  //relative a,b,c index for 8 patches in ComputeLCPO
154  int idx[8][3] = {
155  { 0, 0, 0},
156  { 1, 0, 0},
157  { 0, 1, 0},
158  { 1, 1, 0},
159  { 0, 0, 1},
160  { 1, 0, 1},
161  { 0, 1, 1},
162  { 1, 1, 1} };
163 /*
164  int i_a = patchMap->index_a(patchID[0]);
165  int i_b = patchMap->index_b(patchID[0]);
166  int i_c = patchMap->index_c(patchID[0]);
167  CkPrintf("VALID[%d,%d,%d]=\n",i_a,i_b,i_c);
168 */
169  for (int pI = 0; pI < 8; pI++) {
170  int iia = patchMap->index_a(patchID[pI]);
171  int iib = patchMap->index_b(patchID[pI]);
172  int iic = patchMap->index_c(patchID[pI]);
173  for (int pJ = 0; pJ < 8; pJ++) {
174  int jia = patchMap->index_a(patchID[pJ]);
175  int jib = patchMap->index_b(patchID[pJ]);
176  int jic = patchMap->index_c(patchID[pJ]);
177  if ( ( gsa==1 && (jia>iia) != (idx[pJ][0]>idx[pI][0]) ) ||
178  ( gsb==1 && (jib>iib) != (idx[pJ][1]>idx[pI][1]) ) ||
179  ( gsc==1 && (jic>iic) != (idx[pJ][2]>idx[pI][2]) ) ||
180  ( invalidPatch[pI] ) ||
181  ( invalidPatch[pJ] ) )
182  valid[pI][pJ] = 0;
183  else
184  valid[pI][pJ] = 1;
185  //CkPrintf("%d ",valid[pI][pJ]);
186  }
187  //CkPrintf("\n");
188  }
189  //CkPrintf("\n");
190 
191 } // initialize
192 
194  for (int i=0; i<8; i++) {
195  numAtoms[i] = patch[i]->getNumAtoms();
196  }
197 }
198 
199 //---------------------------------------------------------------------
200 // doWork
201 //---------------------------------------------------------------------
204  for (int i=0; i<8; i++) {
205  pos[i] = positionBox[i]->open();
206  force[i] = forceBox[i]->open();
207  posExt[i] = patch[i]->getCompAtomExtInfo();
208  lcpoType[i] = lcpoTypeBox[i]->open();
209  }
210 
211  doForce();
212 
213  // Inform load balancer
215 
216  // Close up boxes
217  for (int i=0; i<getNumPatches(); i++) {
218  positionBox[i]->close(&pos[i]);
219  forceBox[i]->close(&force[i]);
220  lcpoTypeBox[i]->close(&lcpoType[i]);
221  }
222 } // doWork
223 
224 
226 
227  if ( patch[0]->flags.doNonbonded) {
228  return 0; // work to do, enqueue as usual
229  } else {
230 
231  // skip all boxes
232  for (int i=0; i<8; i++) {
233  positionBox[i]->skip();
234  forceBox[i]->skip();
235  lcpoTypeBox[i]->skip();
236  }
237 
239  reduction->submit();
241 
242  return 1; // no work to do, do not enqueue
243  }
244  return 0;
245 } // noWork
246 
247 // 1 - yes in bounds, 0 - not in bounds
248 // this does uniquely assign atoms to ComputeLCPO octets
249 int ComputeLCPO::isInBounds(Real x, Real y, Real z ) {
250 
251  //check x dimension
252  if ( bounds[0][0] < bounds[0][1] ) { // internal
253  if (x < bounds[0][0] || x >= bounds[0][1] )
254  return 0;
255  } else { // edge
256  if (x < bounds[0][0] && x >= bounds[0][1] )
257  return 0;
258  }
259 
260  //check y dimension
261  if ( bounds[1][0] < bounds[1][1] ) { // internal
262  if (y < bounds[1][0] || y >= bounds[1][1] )
263  return 0;
264  } else { // edge
265  if (y < bounds[1][0] && y >= bounds[1][1] )
266  return 0;
267  }
268 
269  //check z dimension
270  if ( bounds[2][0] < bounds[2][1] ) { // internal
271  if (z < bounds[2][0] || z >= bounds[2][1] )
272  return 0;
273  } else { // edge
274  if (z < bounds[2][0] && z >= bounds[2][1] )
275  return 0;
276  }
277 
278  return 1;
279 } // isInBounds
280 
281 inline BigReal calcOverlap( BigReal r, Real ri, Real rj ) {
282  return PI*ri*(2*ri-r-(ri*ri-rj*rj)/r);
283 }
284 
291  //CkPrintf("ComputeLCPO::doForce\n");
292  step = patch[0]->flags.sequence;
293 
294  Real probeRadius = 1.4f;
295  Real cutMargin = 0.f;//regenerating pairlists every step
296 
297  Position ngir, ngjr, ngkr;
298  Real ri, rj, rk;
299  BigReal dxij, dyij, dzij, r2ij;
300  BigReal dxik, dyik, dzik, r2ik;
301  BigReal dxjk, dyjk, dzjk, r2jk;
302 
303 #ifdef COUNT_FLOPS
304  int flops = 0;
305 #endif
306 
308 // Build Pairlists
310 //generate pairlists every step since contain coordinates
311 if ( true ) {
312  double t_start = 1.0*clock()/CLOCKS_PER_SEC;
313 
314  inAtomsPl.reset();
315  lcpoNeighborList.reset();
316  FLOPS(8);
317  cut2 = 2*maxAtomRadius+cutMargin; cut2 *= cut2;
318  maxAtomRadius = 0;
319  //find in-bounds atoms in each patch
320  for (int pI = 0; pI < 8; pI++) {
321  if (invalidPatch[pI]) continue;
322  if (numAtoms[pI] == 0) continue;
323 
324  int minIg = 0;
325  for (int s = 0; s < minPart; s++) {
326  minIg += pos[pI][minIg].hydrogenGroupSize;
327  FLOPS(1)
328  }
329  strideIg = numParts;//stride through partitions
330  plint *inAtoms = inAtomsPl.newlist(numAtoms[pI]);
331  int numAtomsInBounds = 0;
332 
333  //iterate over heavy atoms only
334  for ( int ngi = minIg; ngi < numAtoms[pI]; /* ngi */) {
335  ngir = pos[pI][ngi].position;
336  if ( isInBounds(ngir.x, ngir.y, ngir.z) && lcpoType[pI][ngi] > 0 ) {
337  inAtoms[numAtomsInBounds++] = ngi;
338  ri = probeRadius+lcpoParams[ lcpoType[pI][ngi] ][0];
339  maxAtomRadius = (ri > maxAtomRadius) ? ri : maxAtomRadius;
340  FLOPS(1);
341 
342  int maxAtoms = 0;
343  for (int pJ = 0; pJ < 8; pJ++) {
344  if (numAtoms[pJ] > 0) {
345  maxAtoms += numAtoms[pJ];
346  }
347  }
348  LCPOAtom *lcpoNeighbors = lcpoNeighborList.newlist(maxAtoms);
349  int numLcpoNeighbors = 0;
350 
351  //find pairs of this inAtom from all 8 patches
352  for (int pJ = 0; pJ < 8; pJ++) {
353  if (invalidPatch[pJ]) continue;
354  if (!valid[pI][pJ]) continue;
355 
356  // j atom pairs
357  for ( int ngj = 0; ngj < numAtoms[pJ]; /* ngj */) {
358  FLOPS(1)
359  ngjr = pos[pJ][ngj].position;
360  dxij = ngir.x - ngjr.x;
361  dyij = ngir.y - ngjr.y;
362  dzij = ngir.z - ngjr.z;
363 
364  // i-j coarse check if too far apart
365  r2ij = dxij*dxij + dyij*dyij + dzij*dzij;
366  FLOPS(8)
367  if (r2ij < cut2 && r2ij > 0.01) {
368 
369  // i-j precise check if too far apart
370  rj = probeRadius+lcpoParams[ lcpoType[pJ][ngj] ][0];
371  FLOPS(5)
372  BigReal rirjcutMargin2 = ri+rj+cutMargin;
373  rirjcutMargin2 *= rirjcutMargin2;
374  if (r2ij < rirjcutMargin2 && r2ij > 0.0001 &&
375  lcpoType[pJ][ngj] > 0) {
376  lcpoNeighbors[numLcpoNeighbors].x = ngjr.x;
377  lcpoNeighbors[numLcpoNeighbors].y = ngjr.y;
378  lcpoNeighbors[numLcpoNeighbors].z = ngjr.z;
379  lcpoNeighbors[numLcpoNeighbors].r = rj;
380  lcpoNeighbors[numLcpoNeighbors].f =
381  &force[pJ]->f[Results::nbond][ngj];
382  numLcpoNeighbors++;
383  FLOPS(2)
384  maxAtomRadius = (rj > maxAtomRadius) ? rj : maxAtomRadius;
385  } // precise cutoff
386  } // coarse cutoff
387  //jump to next nonbonded group
388  ngj += pos[pJ][ngj].hydrogenGroupSize;
389  FLOPS(1)
390  } // for j atoms
391  } // for patches J
392  lcpoNeighborList.newsize(numLcpoNeighbors);
393  } // in bounds
394  //jump to next nonbonded group for round-robin
395  for (int s = 0; s < strideIg; s++) {
396  ngi += pos[pI][ngi].hydrogenGroupSize;
397  FLOPS(1)
398  }
399  } // for i atoms
400  inAtomsPl.newsize(numAtomsInBounds);
401  } // for patches I
402 #ifdef COUNT_FLOPS
403  double t_stop = 1.0*clock()/CLOCKS_PER_SEC;
404  CkPrintf("LCPO_TIME_P %7.3f Gflops %9d @ %f\n", flops*1e-9/(t_stop-t_start),flops,(t_stop-t_start));
405 #endif
406 }
407 #ifdef COUNT_FLOPS
408  double t_start = 1.0*clock()/CLOCKS_PER_SEC;
409  flops = 0;
410 #endif
411 
412  //reset pairlists
413  inAtomsPl.reset();
414  lcpoNeighborList.reset();
415  cut2 = maxAtomRadius*2; cut2 *= cut2;
416 
417  //init values
418  BigReal totalSurfaceArea = 0;
419 
427  //for each patch in octet
428  for (int pI = 0; pI < 8; pI++) {
429  if (invalidPatch[pI]) continue;
430  if (numAtoms[pI] == 0) continue;
431  plint *inAtoms;
432  int numInAtoms;
433  inAtomsPl.nextlist( &inAtoms, &numInAtoms );
434  //for each inAtom in each patch
435  for (int i = 0; i < numInAtoms; i++) {
436  int iIndex = inAtoms[i];
437  int idi = posExt[pI][iIndex].id;
438  Real xi = pos[pI][iIndex].position.x;
439  Real yi = pos[pI][iIndex].position.y;
440  Real zi = pos[pI][iIndex].position.z;
441  const Real *lcpoParamI = lcpoParams[ lcpoType[pI][iIndex] ];
442  ri = probeRadius+lcpoParamI[0];
443  FLOPS(1)
444 
445  Real P1 = lcpoParamI[1];
446  Real P2 = lcpoParamI[2];
447  Real P3 = lcpoParamI[3];
448  Real P4 = lcpoParamI[4];
449 
451 // S1
453  BigReal S1 = 4.0*PI*ri*ri; // a
454  FLOPS(3)
455 
456  //for surface area calculation
457  BigReal AijSum = 0; // b
458  BigReal AjkSum = 0; // c
459  BigReal AjkjSum = 0; // d'
460  BigReal AijAjkSum = 0; // d
461 
462  //for force calculation
463  BigReal dAijdrijdxiSum = 0.0;
464  BigReal dAijdrijdyiSum = 0.0;
465  BigReal dAijdrijdziSum = 0.0;
466  BigReal dAijdrijdxiAjkSum = 0.0;
467  BigReal dAijdrijdyiAjkSum = 0.0;
468  BigReal dAijdrijdziAjkSum = 0.0;
469 
471 // for J Atoms
473  LCPOAtom *lcpoNeighbors;
474  int numLcpoNeighbors;
475  lcpoNeighborList.nextlist( &lcpoNeighbors, &numLcpoNeighbors );
476 
477  for (int j = 0; j < numLcpoNeighbors; j++) {
478  Real xj = lcpoNeighbors[j].x;
479  Real yj = lcpoNeighbors[j].y;
480  Real zj = lcpoNeighbors[j].z;
481  Real rj = lcpoNeighbors[j].r;
482 
483  // i-j coarse check if too far away
484  dxij = xj-xi;
485  dyij = yj-yi;
486  dzij = zj-zi;
487  r2ij = dxij*dxij + dyij*dyij + dzij*dzij;
488  FLOPS(7);
489  if (r2ij >= cut2 || r2ij < 0.01) { continue; }
490 
491  // i-j precise check if too far away
492  FLOPS(5)
493  BigReal rirj2 = ri+rj;
494  rirj2 *= rirj2;
495  if ( r2ij >= rirj2 ) { continue; }
496 
497  BigReal rij = sqrt(r2ij);
498  BigReal rij_1 = 1.f / rij;
499 
501 // S2
503  BigReal Aij = calcOverlap(rij, ri, rj);
504  AijSum += Aij;
505  FLOPS(12)
506 
507  //for dAi_drj force calculation
508  BigReal dAijdrij = PI*ri*(rij_1*rij_1*(ri*ri-rj*rj)-1);
509  BigReal dAijdrijdxj = dAijdrij*dxij*rij_1; // g k' i' l'
510  BigReal dAijdrijdyj = dAijdrij*dyij*rij_1;
511  BigReal dAijdrijdzj = dAijdrij*dzij*rij_1;
512  FLOPS(14)
513 
514  BigReal AjkjSum = 0; // i' l'
515  BigReal dAjkdrjkdxjSum = 0.0;
516  BigReal dAjkdrjkdyjSum = 0.0;
517  BigReal dAjkdrjkdzjSum = 0.0;
518 
520 // for K Atoms
522  for (int k = 0; k < numLcpoNeighbors; k++) {
523  Real xk = lcpoNeighbors[k].x;
524  Real yk = lcpoNeighbors[k].y;
525  Real zk = lcpoNeighbors[k].z;
526  Real rk = lcpoNeighbors[k].r;
527 
528  // i-k coarse check if too far away
529  dxik = xk-xi;
530  dyik = yk-yi;
531  dzik = zk-zi;
532  r2ik = dxik*dxik + dyik*dyik + dzik*dzik;
533  FLOPS(8)
534  if (r2ik >= cut2 || r2ik < 0.01) { continue; }
535 
536  // j-k coarse check if too far away
537  dxjk = xk-xj;
538  dyjk = yk-yj;
539  dzjk = zk-zj;
540  r2jk = dxjk*dxjk + dyjk*dyjk + dzjk*dzjk;
541  FLOPS(8)
542  if (r2jk >= cut2 || r2jk < 0.01) { continue; }
543 
544  // i-k precise check if too far away
545  FLOPS(3)
546  BigReal rirk2 = ri+rk;
547  rirk2 *= rirk2;
548  if ( r2ik >= rirk2 ) { continue; }
549 
550  // j-k precise check if too far away
551  FLOPS(2)
552  BigReal rjrk2 = rj+rk;
553  rjrk2 *= rjrk2;
554  if ( r2jk >= rjrk2 ) { continue; }
555  BigReal rjk = sqrt(r2jk);
556 
558 // S3
560  BigReal rjk_1 = 1.0/rjk;
561  BigReal Ajk = calcOverlap(rjk, rj, rk);
562  FLOPS(12)
563  AjkSum += Ajk;
564  AjkjSum += Ajk; // i' l'
565  FLOPS(5)
566 
567 // Force dAi_drk
570  BigReal dAjkdrjk = PI*rj*rjk_1*(rjk_1*rjk_1*(rj*rj-rk*rk) - 1.f);//ef'
571  BigReal dAjkdrjkdxj = -dAjkdrjk*dxjk; // e f h'
572  BigReal dAjkdrjkdyj = -dAjkdrjk*dyjk;
573  BigReal dAjkdrjkdzj = -dAjkdrjk*dzjk;
574  lcpoNeighbors[k].f->x -= -dAjkdrjkdxj*(P3+P4*Aij)*surfTen; // e f
575  lcpoNeighbors[k].f->y -= -dAjkdrjkdyj*(P3+P4*Aij)*surfTen;
576  lcpoNeighbors[k].f->z -= -dAjkdrjkdzj*(P3+P4*Aij)*surfTen;
577 
578  dAjkdrjkdxjSum += dAjkdrjkdxj; // h j'
579  dAjkdrjkdyjSum += dAjkdrjkdyj;
580  dAjkdrjkdzjSum += dAjkdrjkdzj;
581  FLOPS(34)
582 
583  } // k atoms
585 // S4
587  AijAjkSum += Aij*AjkjSum;
588 
590 // Force dAi_drj
592  BigReal lastxj = dAijdrijdxj*AjkjSum + Aij*dAjkdrjkdxjSum; // i j
593  BigReal lastyj = dAijdrijdyj*AjkjSum + Aij*dAjkdrjkdyjSum;
594  BigReal lastzj = dAijdrijdzj*AjkjSum + Aij*dAjkdrjkdzjSum;
595  BigReal dAidxj = (P2*dAijdrijdxj + P3*dAjkdrjkdxjSum + P4*lastxj);//ghij
596  BigReal dAidyj = (P2*dAijdrijdyj + P3*dAjkdrjkdyjSum + P4*lastyj);
597  BigReal dAidzj = (P2*dAijdrijdzj + P3*dAjkdrjkdzjSum + P4*lastzj);
598  lcpoNeighbors[j].f->x -= dAidxj*surfTen;
599  lcpoNeighbors[j].f->y -= dAidyj*surfTen;
600  lcpoNeighbors[j].f->z -= dAidzj*surfTen;
601 
602  //for dAi_dri force calculation
603  dAijdrijdxiSum -= dAijdrijdxj; // k
604  dAijdrijdyiSum -= dAijdrijdyj;
605  dAijdrijdziSum -= dAijdrijdzj;
606  dAijdrijdxiAjkSum -= dAijdrijdxj*AjkjSum; // l
607  dAijdrijdyiAjkSum -= dAijdrijdyj*AjkjSum;
608  dAijdrijdziAjkSum -= dAijdrijdzj*AjkjSum;
609  FLOPS(41)
610  } // j atoms
611 
613 // Force dAi_dri
615  BigReal dAidxi = (P2*dAijdrijdxiSum + P4*dAijdrijdxiAjkSum); // k l
616  BigReal dAidyi = (P2*dAijdrijdyiSum + P4*dAijdrijdyiAjkSum);
617  BigReal dAidzi = (P2*dAijdrijdziSum + P4*dAijdrijdziAjkSum);
618  force[pI]->f[Results::nbond][iIndex].x -= dAidxi*surfTen;
619  force[pI]->f[Results::nbond][iIndex].y -= dAidyi*surfTen;
620  force[pI]->f[Results::nbond][iIndex].z -= dAidzi*surfTen;
621 
623 // Atom I Surface Area
625  BigReal SAi = P1*S1 + P2*AijSum + P3*AjkSum + P4*AijAjkSum;
626  //CkPrintf("SurfArea[%05d] = % 7.3f\n",idi,SAi);
627  //SAi = (SAi > 0) ? SAi : 0;
628  totalSurfaceArea += SAi;
629  FLOPS(22)
630  } // for inAtoms
631  } // for patches I
632 #ifdef COUNT_FLOPS
633  double t_stop = 1.0*clock()/CLOCKS_PER_SEC;
634  CkPrintf("LCPO_TIME_F %7.3f Gflops %9d @ %f\n", 1e-9*flops/(t_stop-t_start),flops, (t_stop-t_start));
635 #endif
636 
638 // end calculation by submitting reduction
640 
642  reduction->item(REDUCTION_ELECT_ENERGY) += totalSurfaceArea * surfTen;
643  reduction->submit();
644 
645 }//end do Force
646 
647 
648 // Lookup table for lcpo paramters
649 //indices 0 -> 22 are determined in Molecule.C
650 const Real ComputeLCPO::lcpoParams[23][5] = { // neigh
651  { 0.00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00 }, // 0 H
652  { 1.70, 7.7887e-01, -2.8063e-01, -1.2968e-03, 3.9328e-04 }, // 1 C sp3 1
653  { 1.70, 5.6482e-01, -1.9608e-01, -1.0219e-03, 2.6580e-04 }, // 2 C sp3 2
654  { 1.70, 2.3348e-01, -7.2627e-02, -2.0079e-04, 7.9670e-05 }, // 3 C sp3 3
655  { 1.70, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00 }, // 4 C sp3 4
656  { 1.70, 5.1245e-01, -1.5966e-01, -1.9781e-04, 1.6392e-04 }, // 5 C sp2 2
657  { 1.70, 7.0344e-02, -1.9015e-02, -2.2009e-05, 1.6875e-05 }, // 6 C sp2 3
658  { 1.60, 7.7914e-01, -2.5262e-01, -1.6056e-03, 3.5071e-04 }, // 7 O sp3 1
659  { 1.60, 4.9392e-01, -1.6038e-01, -1.5512e-04, 1.6453e-04 }, // 8 O sp3 2
660  { 1.60, 6.8563e-01, -1.8680e-01, -1.3557e-03, 2.3743e-04 }, // 9 O sp2 1
661  { 1.60, 8.8857e-01, -3.3421e-01, -1.8683e-03, 4.9372e-04 }, // 10 O=C-O
662  { 1.65, 7.8602e-02, -2.9198e-01, -6.5370e-04, 3.6247e-04 }, // 11 N sp3 1
663  { 1.65, 2.2599e-01, -3.6648e-02, -1.2297e-03, 8.0038e-05 }, // 12 N sp3 2
664  { 1.65, 5.1481e-02, -1.2603e-02, -3.2006e-04, 2.4774e-05 }, // 13 N sp3 3
665  { 1.65, 7.3511e-01, -2.2116e-01, -8.9148e-04, 2.5230e-04 }, // 14 N sp2 1
666  { 1.65, 4.1102e-01, -1.2254e-01, -7.5448e-05, 1.1804e-04 }, // 15 N sp2 2
667  { 1.65, 6.2577e-02, -1.7874e-02, -8.3120e-05, 1.9849e-05 }, // 16 N sp2 3
668  { 1.90, 7.7220e-01, -2.6393e-01, 1.0629e-03, 2.1790e-04 }, // 17 S 1
669  { 1.90, 5.4581e-01, -1.9477e-01, -1.2873e-03, 2.9247e-04 }, // 18 S 2
670  { 1.90, 3.8650e-01, -1.8249e-01, -3.6598e-03, 4.2640e-04 }, // 19 P 3
671  { 1.90, 3.8730e-02, -8.9339e-03, 8.3582e-06, 3.0381e-06 }, // 20 P 4
672  { 1.80, 9.8318e-01, -4.0437e-01, 1.1249e-04, 4.9901e-04 }, // 21 Cl
673  { 1.18, 4.9392e-01, -1.6038e-01, -1.5512e-04, 1.6453e-04 } // 22 Mg
674 };
static Node * Object()
Definition: Node.h:86
void setNumPatches(int n)
Definition: Compute.h:52
#define PROXY_RESULTS_PRIORITY
Definition: Priorities.h:73
float r
Definition: ComputeLCPO.h:23
virtual int noWork()
Definition: ComputeLCPO.C:225
unsigned int hydrogenGroupSize
Definition: NamdTypes.h:58
int ComputeID
Definition: NamdTypes.h:183
Box< Patch, int > * lcpoTypeBox[8]
Definition: ComputeLCPO.h:132
int gridsize_c(void) const
Definition: PatchMap.h:66
static PatchMap * Object()
Definition: PatchMap.h:27
BigReal max_c(int pid) const
Definition: PatchMap.h:96
__global__ void const int const TileList *__restrict__ TileExcl *__restrict__ const int *__restrict__ const int const float2 *__restrict__ const int *__restrict__ const float3 const float3 const float3 const float4 *__restrict__ const float cudaTextureObject_t cudaTextureObject_t cudaTextureObject_t const int const float const PatchPairRecord *__restrict__ const int *__restrict__ const int2 *__restrict__ const unsigned int *__restrict__ unsigned int *__restrict__ int *__restrict__ int *__restrict__ TileListStat *__restrict__ const BoundingBox *__restrict__ float4 *__restrict__ float4 *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ const int numPatches
Definition: Vector.h:64
BigReal min_a(int pid) const
Definition: PatchMap.h:91
Vector * f
Definition: ComputeLCPO.h:24
SimParameters * simParameters
Definition: Node.h:178
int index_a(int pid) const
Definition: PatchMap.h:86
float Real
Definition: common.h:107
BigReal & item(int i)
Definition: ReductionMgr.h:312
#define DebugM(x, y)
Definition: Debug.h:59
void unregisterForceDeposit(Compute *cid, Box< Patch, Results > **const box)
Definition: Patch.C:239
void startWork(const LDObjHandle &handle)
BigReal z
Definition: Vector.h:66
Results * force[8]
Definition: ComputeLCPO.h:121
virtual void initialize()
Definition: Compute.h:56
Position position
Definition: NamdTypes.h:53
CompAtomExt * posExt[8]
Definition: ComputeLCPO.h:119
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:365
#define FLOPS(X)
Definition: ComputeLCPO.C:30
LDObjHandle ldObjHandle
Definition: Compute.h:44
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:278
Box< Patch, int > * registerLcpoTypePickup(Compute *cid)
Definition: Patch.C:188
int * lcpoType[8]
Definition: ComputeLCPO.h:122
Patch * patch(PatchID pid)
Definition: PatchMap.h:235
BigReal min_b(int pid) const
Definition: PatchMap.h:93
Box< Patch, Results > * forceBox[8]
Definition: ComputeLCPO.h:131
Flags flags
Definition: Patch.h:127
float x
Definition: ComputeLCPO.h:23
virtual void atomUpdate()
Definition: ComputeLCPO.C:193
virtual ~ComputeLCPO()
Definition: ComputeLCPO.C:61
#define PI
Definition: common.h:81
float z
Definition: ComputeLCPO.h:23
unsigned short plint
void nextlist(plint **list, int *list_size)
int gridsize_a(void) const
Definition: PatchMap.h:64
int valid[8][8]
Definition: ComputeLCPO.h:116
virtual void initialize()
Definition: ComputeLCPO.C:86
float y
Definition: ComputeLCPO.h:23
BigReal surface_tension
virtual void doWork()
Definition: ComputeLCPO.C:202
void skipWork(const LDObjHandle &handle)
gridSize z
int numAtoms[8]
Definition: ComputeLCPO.h:115
Force * f[maxNumForces]
Definition: PatchTypes.h:67
int index_b(int pid) const
Definition: PatchMap.h:87
BigReal x
Definition: Vector.h:66
PatchID patchID[8]
Definition: ComputeLCPO.h:128
void newsize(int list_size)
Definition: ComputeLCPO.h:75
Box< Patch, CompAtom > * positionBox[8]
Definition: ComputeLCPO.h:130
int sequence
Definition: PatchTypes.h:18
int PatchID
Definition: NamdTypes.h:182
static LdbCoordinator * Object()
LCPOAtom * newlist(int max_size)
Definition: ComputeLCPO.h:48
BigReal max_b(int pid) const
Definition: PatchMap.h:94
void skip(void)
Definition: Box.h:63
int index_c(int pid) const
Definition: PatchMap.h:88
ComputeLCPO(ComputeID c, PatchID pid[], int t[], ComputeNonbondedWorkArrays *_workArrays, int minPartition, int maxPartition, int numPartitions, int numPatches)
Definition: ComputeLCPO.C:36
Patch * patch[8]
Definition: ComputeLCPO.h:126
void unregisterPositionPickup(Compute *cid, Box< Patch, CompAtom > **const box)
Definition: Patch.C:122
BigReal max_a(int pid) const
Definition: PatchMap.h:92
void endWork(const LDObjHandle &handle)
#define simParams
Definition: Output.C:127
CompAtom * pos[8]
Definition: ComputeLCPO.h:120
BigReal calcOverlap(BigReal r, Real ri, Real rj)
Definition: ComputeLCPO.C:281
int getNumPatches()
Definition: Compute.h:53
BigReal y
Definition: Vector.h:66
int getNumAtoms()
Definition: Patch.h:105
virtual void doForce()
Definition: ComputeLCPO.C:290
int trans[8]
Definition: ComputeLCPO.h:129
void newsize(int list_size)
gridSize y
void submit(void)
Definition: ReductionMgr.h:323
int basePriority
Definition: Compute.h:37
BigReal min_c(int pid) const
Definition: PatchMap.h:95
void nextlist(LCPOAtom **list, int *list_size)
Definition: ComputeLCPO.h:84
int invalidPatch[8]
Definition: ComputeLCPO.h:118
gridSize x
plint * newlist(int max_size)
SubmitReduction * reduction
Definition: ComputeLCPO.h:137
Data * open(void)
Definition: Box.h:39
int gridsize_b(void) const
Definition: PatchMap.h:65
const ComputeID cid
Definition: Compute.h:43
void unregisterLcpoTypePickup(Compute *cid, Box< Patch, int > **const box)
Definition: Patch.C:191
void close(Data **const t)
Definition: Box.h:49
Box< Patch, CompAtom > * registerPositionPickup(Compute *cid)
Definition: Patch.C:107
double BigReal
Definition: common.h:112
#define PATCH_PRIORITY(PID)
Definition: Priorities.h:25
CompAtomExt * getCompAtomExtInfo()
Definition: Patch.h:117
Box< Patch, Results > * registerForceDeposit(Compute *cid)
Definition: Patch.C:228