NAMD
ComputeNonbondedSelf.C
Go to the documentation of this file.
1 
7 #include "ComputeNonbondedSelf.h"
8 #include "ReductionMgr.h"
9 #include "Patch.h"
10 #include "LdbCoordinator.h"
11 #include "Molecule.h"
12 
13 #include "Node.h"
14 #include "SimParameters.h"
15 
16 #include "Parameters.h"
17 #include "LJTable.h"
18 
19 #define MIN_DEBUG_LEVEL 4
20 // #define DEBUGM
21 #include "Debug.h"
22 
24  ComputeNonbondedWorkArrays* _workArrays,
25  int minPartition, int maxPartition, int numPartitions)
26  : ComputePatch(c,pid), workArrays(_workArrays),
27  minPart(minPartition), maxPart(maxPartition), numParts(numPartitions)
28 {
30  if (pressureProfileOn) {
35  } else {
37  pressureProfileData = NULL;
38  }
39  pairlistsValid = 0;
40  pairlistTolerance = 0.;
44 }
45 
46 #ifdef NAMD_AVXTILES
49 
50  if (avxTilesMode) {
51  tileLists.setSimParams(scaling, scale14, c1, c3, switchOn2,
52  knl_fast_grad_table, knl_fast_ener_table,
53  knl_scor_grad_table, knl_scor_ener_table,
54  avx_tiles_eps4_sigma, avx_tiles_eps4_sigma_14,
55  (float *)ljTable->table_val(0,0),
56  ljTable->get_table_dim(), knl_slow_grad_table,
57  knl_slow_ener_table, knl_excl_grad_table,
58  knl_excl_ener_table, avxTilesMode);
59  tileLists.atomUpdate(patch->getTiles(), patch->getTiles());
60  }
61 }
62 #endif
63 
67  // BEGIN LA
69  // END LA
70 
76 
77 #ifdef NAMD_CUDA
79 #endif
80 }
81 
83 {
84  delete reduction;
86  delete [] pressureProfileData;
87  if (avgPositionBox != NULL) {
89  }
90  // BEGIN LA
91  if (velocityBox != NULL) {
93  }
94  // END LA
95 
96  if (psiSumBox != NULL)
98  if (intRadBox != NULL)
100  if (bornRadBox != NULL)
102  if (dEdaSumBox != NULL)
104  if (dHdrPrefixBox != NULL)
106 }
107 
109 
110  if (patch->flags.doGBIS) {
111  gbisPhase = 1 + (gbisPhase % 3);//1->2->3->1...
112  }
113 
114 #ifndef NAMD_CUDA
115  if ( patch->flags.doNonbonded && numAtoms ) {
116  return 0; // work to do, enqueue as usual
117  } else {
118 #else
119  {
120 #endif
121 
122  if (patch->flags.doGBIS) {
123  if (gbisPhase == 1) {
124  psiSumBox->skip();
125  intRadBox->skip();
126  if (patch->flags.doNonbonded) return 1;
127  else gbisPhase = 2;
128  }
129  if (gbisPhase == 2) {
130  bornRadBox->skip();
131  dEdaSumBox->skip();
132  if (patch->flags.doNonbonded) return 1;
133  else gbisPhase = 3;
134  }
135  if (gbisPhase == 3) {
136  dHdrPrefixBox->skip();
137  }
138  }
139 
140  // skip all boxes
141  positionBox->skip();
142  forceBox->skip();
144  // BEGIN LA
146  // END LA
147 
149  reduction->submit();
150 
151  if (pressureProfileOn)
153 
154 #ifndef NAMD_CUDA
155  // Inform load balancer
157 #endif
158 
159  return 1; // no work to do, do not enqueue
160  }
161 }
162 
164 {
165  #ifdef NAMD_AVXTILES
166  if (avxTilesMode) {
167  doForceTiles(p, pExt, r);
168  return;
169  }
170  #endif
171 
172  // Inform load balancer.
173  // I assume no threads will suspend until endWork is called
174  //single phase declarations
175  CompAtom* v;
176  int doEnergy = patch->flags.doEnergy;
177 
178 
179 /*******************************************************************************
180  * Prepare Parameters
181  ******************************************************************************/
182  if (!patch->flags.doGBIS || gbisPhase == 1) {
183 
184 #ifdef TRACE_COMPUTE_OBJECTS
185  double traceObjStartTime = CmiWallTimer();
186 #endif
187 
188  DebugM(2,"doForce() called.\n");
189  DebugM(1,numAtoms << " patch 1 atoms\n");
190  DebugM(3, "NUMATOMSxNUMATOMS = " << numAtoms*numAtoms << "\n");
191 
192  for ( int i = 0; i < reductionDataSize; ++i ) reductionData[i] = 0;
193  if (pressureProfileOn) {
194  int n = pressureProfileAtomTypes;
195  memset(pressureProfileData, 0, 3*n*n*pressureProfileSlabs*sizeof(BigReal));
196  // adjust lattice dimensions to allow constant pressure
197  const Lattice &lattice = patch->lattice;
199  pressureProfileMin = lattice.origin().z - 0.5*lattice.c().z;
200  }
201 
202  plint maxa = (plint)(-1);
203  if ( numAtoms > maxa ) {
204  char estr[1024];
205  sprintf(estr,"patch has %d atoms, maximum allowed is %d",numAtoms,maxa);
206  NAMD_die(estr);
207  }
208 
209  params.offset = 0.;
210  params.offset_f = 0.;
211  params.p[0] = p;
212  params.p[1] = p;
213  params.pExt[0] = pExt;
214  params.pExt[1] = pExt;
215 #ifdef NAMD_KNL
216  CompAtomFlt *pFlt = patch->getCompAtomFlt();
217  params.pFlt[0] = pFlt;
218  params.pFlt[1] = pFlt;
219 #endif
221  // BEGIN LA
223  if (params.doLoweAndersen) {
224  DebugM(4, "opening velocity box\n");
225  v = velocityBox->open();
226  params.v[0] = v;
227  params.v[1] = v;
228  }
229  // END LA
230 #ifndef NAMD_CUDA
231  params.ff[0] = r->f[Results::nbond_virial];
232  params.ff[1] = r->f[Results::nbond_virial];
233 #endif
234  params.numAtoms[0] = numAtoms;
235  params.numAtoms[1] = numAtoms;
236 
237  // DMK - Atom Separation (water vs. non-water)
238  #if NAMD_SeparateWaters != 0
239  params.numWaterAtoms[0] = numWaterAtoms;
240  params.numWaterAtoms[1] = numWaterAtoms;
241  #endif
242 
245 
249 
251 
253  params.savePairlists = 0;
254  params.usePairlists = 0;
255  if ( patch->flags.savePairlists ) {
256  params.savePairlists = 1;
257  params.usePairlists = 1;
258  } else if ( patch->flags.usePairlists ) {
259  if ( ! pairlistsValid ||
262  } else {
263  params.usePairlists = 1;
264  }
265  }
266  if ( ! params.usePairlists ) {
267  pairlistsValid = 0;
268  }
271  if ( params.savePairlists ) {
272  pairlistsValid = 1;
276  }
277 
278 
279 /*******************************************************************************
280  * Call Nonbonded Functions
281  ******************************************************************************/
282 
283  if (numAtoms) {
285  {
286 #ifndef NAMD_CUDA
287  params.fullf[0] = r->f[Results::slow_virial];
288  params.fullf[1] = r->f[Results::slow_virial];
289 #endif
290  if ( patch->flags.doMolly ) {
291  if ( doEnergy ) calcSelfEnergy(&params);
292  else calcSelf(&params);
293  CompAtom *p_avg = avgPositionBox->open();
294  params.p[0] = p_avg;
295  params.p[1] = p_avg;
296  if ( doEnergy ) calcSlowSelfEnergy(&params);
297  else calcSlowSelf(&params);
298  avgPositionBox->close(&p_avg);
299  } else if ( patch->flags.maxForceMerged == Results::slow ) {
300  if ( doEnergy ) calcMergeSelfEnergy(&params);
301  else calcMergeSelf(&params);
302  } else {
303  if ( doEnergy ) calcFullSelfEnergy(&params);
304  else calcFullSelf(&params);
305  }
306  }
307  else
308  if ( doEnergy ) calcSelfEnergy(&params);
309  else calcSelf(&params);
310  }//end if atoms
311 
312  // BEGIN LA
313  if (params.doLoweAndersen) {
314  DebugM(4, "closing velocity box\n");
315  velocityBox->close(&v);
316  }
317  // END LA
318  }//end if not gbis
319 
320 /*******************************************************************************
321  * gbis Loop
322 *******************************************************************************/
323 if (patch->flags.doGBIS) {
327  gbisParams.numPatches = 1;//self
331  gbisParams.epsilon_p = simParams->dielectric;
333  gbisParams.kappa = simParams->kappa;
334  gbisParams.cutoff = simParams->cutoff;
336  gbisParams.a_cut = simParams->alpha_cutoff;
337  gbisParams.delta = simParams->gbis_delta;
338  gbisParams.beta = simParams->gbis_beta;
339  gbisParams.gamma = simParams->gbis_gamma;
340  gbisParams.alpha_max = simParams->alpha_max;
341  gbisParams.cid = cid;
345  gbisParams.doEnergy = doEnergy;
346  gbisParams.fsMax = simParams->fsMax;
347  for (int i = 0; i < numGBISPairlists; i++)
349 
350  //open boxes
351  if (gbisPhase == 1) {
358  } else if (gbisPhase == 2) {
363  } else if (gbisPhase == 3) {
366  }
367 
368  //make call to calculate GBIS
371  }
372 
373  //close boxes
374  if (gbisPhase == 1) {
376  } else if (gbisPhase == 2) {
378  } else if (gbisPhase == 3) {
384  }
385 }// end if doGBIS
386 
387 /*******************************************************************************
388  * Reduction
389 *******************************************************************************/
390  if (!patch->flags.doGBIS || gbisPhase == 3) {
392  if (pressureProfileOn)
394 
395 
396 #ifdef TRACE_COMPUTE_OBJECTS
397  traceUserBracketEvent(TRACE_COMPOBJ_IDOFFSET+cid, traceObjStartTime, CmiWallTimer());
398 #endif
399 
400  reduction->submit();
401  if (pressureProfileOn)
403  }// end not gbis
404 }
405 
406 #ifdef NAMD_AVXTILES
407 void ComputeNonbondedSelf::doForceTiles(CompAtom* p, CompAtomExt* pExt,
408  Results* r)
409 {
410  const int doEnergy = patch->flags.doEnergy;
411  const int doVirial = patch->flags.doVirial;
412  const int doFull = patch->flags.doFullElectrostatics;
413 
414 #ifdef TRACE_COMPUTE_OBJECTS
415  double traceObjStartTime = CmiWallTimer();
416 #endif
417 
418  DebugM(2,"doForceTiles() called.\n");
419  DebugM(1,numAtoms << " patch 1 atoms\n");
420  DebugM(3, "NUMATOMSxNUMATOMS = " << numAtoms*numAtoms << "\n");
421 
422  const Lattice &lattice = patch->lattice;
423  for ( int i = 0; i < reductionDataSize; ++i ) reductionData[i] = 0;
424 
425  plint maxa = (plint)(-1);
426  if ( numAtoms > maxa ) {
427  char estr[1024];
428  sprintf(estr,"patch has %d atoms, maximum allowed is %d",numAtoms,maxa);
429  NAMD_die(estr);
430  }
431 
432  int doList = false;
433  if ( patch->flags.savePairlists ) {
434  doList = 1;
435  pairlistsValid = 1;
437  } else if ( patch->flags.usePairlists ) {
438  if (!pairlistsValid ||
441  pairlistsValid = 0;
442  doList = 1;
443  }
444  } else
445  doList = 1;
446 
447  if (doList) {
448  double plcutoff = cutoff;
450  plcutoff += patch->flags.pairlistTolerance * 2.0;
451  tileLists.updateBuildInfo(patch->flags.step, minPart, maxPart,
452  numParts, plcutoff);
453  }
454 
455  Vector offset(0., 0., 0.);
456  tileLists.updateParams(lattice, offset, cutoff);
457 
458  tileLists.nbForceAVX512(doEnergy, doVirial, doFull, doList);
459 
460  reductionData[exclChecksumIndex] = tileLists.exclusionChecksum();
461  if (doEnergy) {
462  reductionData[vdwEnergyIndex] = tileLists.energyVdw();
463  reductionData[electEnergyIndex] = tileLists.energyElec();
464  if (doFull)
466  }
467  if (doVirial) {
468  reductionData[virialIndex_XX] = tileLists.virialXX();
469  reductionData[virialIndex_XY] = tileLists.virialXY();
470  reductionData[virialIndex_XZ] = tileLists.virialXZ();
471  reductionData[virialIndex_YX] = tileLists.virialXY();
472  reductionData[virialIndex_YY] = tileLists.virialYY();
473  reductionData[virialIndex_YZ] = tileLists.virialYZ();
474  reductionData[virialIndex_ZX] = tileLists.virialXZ();
475  reductionData[virialIndex_ZY] = tileLists.virialYZ();
476  reductionData[virialIndex_ZZ] = tileLists.virialZZ();
477  if (doFull) {
478  reductionData[fullElectVirialIndex_XX] = tileLists.virialSlowXX();
479  reductionData[fullElectVirialIndex_XY] = tileLists.virialSlowXY();
480  reductionData[fullElectVirialIndex_XZ] = tileLists.virialSlowXZ();
481  reductionData[fullElectVirialIndex_YX] = tileLists.virialSlowXY();
482  reductionData[fullElectVirialIndex_YY] = tileLists.virialSlowYY();
483  reductionData[fullElectVirialIndex_YZ] = tileLists.virialSlowYZ();
484  reductionData[fullElectVirialIndex_ZX] = tileLists.virialSlowXZ();
485  reductionData[fullElectVirialIndex_ZY] = tileLists.virialSlowYZ();
486  reductionData[fullElectVirialIndex_ZZ] = tileLists.virialSlowZZ();
487  }
488  }
489 
490  if (!patch->flags.doGBIS || gbisPhase == 3) {
492 
493 #ifdef TRACE_COMPUTE_OBJECTS
494  traceUserBracketEvent(TRACE_COMPOBJ_IDOFFSET+cid, traceObjStartTime,
495  CmiWallTimer());
496 #endif
497 
498  reduction->submit();
499  } // end not gbis
500 }
501 #endif
static Node * Object()
Definition: Node.h:86
Box< Patch, GBReal > * registerDEdaSumDeposit(Compute *cid)
Definition: Patch.C:204
Pairlists * pairlists
Box< Patch, CompAtom > * avgPositionBox
Box< Patch, CompAtom > * registerAvgPositionPickup(Compute *cid)
Definition: Patch.C:134
CompAtom * p
Definition: ComputePatch.h:37
int sequence(void)
Definition: Compute.h:64
CompAtomExt * pExt
Definition: ComputePatch.h:36
void unregisterAvgPositionPickup(Compute *cid, Box< Patch, CompAtom > **const box)
Definition: Patch.C:140
BigReal solvent_dielectric
static void submitReductionData(BigReal *, SubmitReduction *)
Parameters * parameters
int ComputeID
Definition: NamdTypes.h:183
#define TRACE_COMPOBJ_IDOFFSET
Definition: Compute.h:77
void unregisterPsiSumDeposit(Compute *cid, Box< Patch, GBReal > **const box)
Definition: Patch.C:174
static void(* calcSelf)(nonbonded *)
Lattice & lattice
Definition: Patch.h:126
CompAtom * p[2]
ComputeNonbondedWorkArrays *const workArrays
Patch * patch
Definition: ComputePatch.h:46
Box< Patch, Real > * dHdrPrefixBox
void calcGBIS(nonbonded *params, GBISParamStruct *gbisParams)
Definition: ComputeGBIS.C:261
Definition: Vector.h:64
SimParameters * simParameters
Definition: Node.h:178
Box< Patch, Results > * forceBox
Definition: ComputePatch.h:51
int savePairlists
Definition: PatchTypes.h:39
BigReal & item(int i)
Definition: ReductionMgr.h:312
GBISParamStruct gbisParams
#define DebugM(x, y)
Definition: Debug.h:59
CompAtom * v[2]
BigReal gbis_gamma
BigReal z
Definition: Vector.h:66
void unregisterDHdrPrefixPickup(Compute *cid, Box< Patch, Real > **const box)
Definition: Patch.C:222
int usePairlists
Definition: PatchTypes.h:38
BigReal * reduction
static void submitPressureProfileData(BigReal *, SubmitReduction *)
static BigReal pressureProfileThickness
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:365
Box< Patch, Real > * registerBornRadPickup(Compute *cid)
Definition: Patch.C:196
int get_table_dim() const
Definition: LJTable.h:44
if(ComputeNonbondedUtil::goMethod==2)
LDObjHandle ldObjHandle
Definition: Compute.h:44
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:278
int doLoweAndersen
Definition: PatchTypes.h:26
Box< Patch, GBReal > * dEdaSumBox
Vector origin() const
Definition: Lattice.h:262
SimParameters * simParameters
BigReal coulomb_radius_offset
virtual void atomUpdate()
Definition: ComputePatch.C:78
Flags flags
Definition: Patch.h:127
void register_cuda_compute_self(ComputeID c, PatchID pid)
__global__ void const int const TileList *__restrict__ tileLists
BigReal alpha_max
Box< Patch, GBReal > * registerPsiSumDeposit(Compute *cid)
Definition: Patch.C:164
unsigned short plint
PatchID patchID
Definition: ComputePatch.h:49
static void(* calcMergeSelfEnergy)(nonbonded *)
const TableEntry * table_val(unsigned int i, unsigned int j) const
Definition: LJTable.h:35
Pairlists * gbisStepPairlists[4]
static BigReal pressureProfileMin
SubmitReduction * pressureProfileReduction
Results * r
Definition: ComputePatch.h:38
int doEnergy
Definition: PatchTypes.h:20
int doFullElectrostatics
Definition: PatchTypes.h:23
Box< Patch, Real > * registerIntRadPickup(Compute *cid)
Definition: Patch.C:179
void skipWork(const LDObjHandle &handle)
Force * f[maxNumForces]
Definition: PatchTypes.h:67
BigReal gbis_beta
SubmitReduction * reduction
CompAtomExt * pExt[2]
ComputeNonbondedSelf(ComputeID c, PatchID pid, ComputeNonbondedWorkArrays *_workArrays, int minPartition=0, int maxPartition=1, int numPartitions=1)
BigReal alpha_cutoff
int PatchID
Definition: NamdTypes.h:182
static void(* calcFullSelf)(nonbonded *)
int doNonbonded
Definition: PatchTypes.h:22
void NAMD_die(const char *err_msg)
Definition: common.C:83
static LdbCoordinator * Object()
static void(* calcSlowSelf)(nonbonded *)
BigReal maxAtomMovement
Definition: PatchTypes.h:41
void skip(void)
Definition: Box.h:63
void unregisterBornRadPickup(Compute *cid, Box< Patch, Real > **const box)
Definition: Patch.C:199
static void(* calcSlowSelfEnergy)(nonbonded *)
Box< Patch, Real > * bornRadBox
int gbisPhase
Definition: Compute.h:39
Parameters * parameters
Definition: Node.h:177
Box< Patch, GBReal > * psiSumBox
static void(* calcSelfEnergy)(nonbonded *)
PatchID getPatchID()
Definition: Patch.h:114
virtual void doForce(CompAtom *p, CompAtomExt *pExt, Results *r)
Box< Patch, Real > * intRadBox
#define simParams
Definition: Output.C:127
virtual void initialize()
Definition: ComputePatch.C:46
Random * rand
Definition: Node.h:172
Box< Patch, CompAtom > * positionBox
Definition: ComputePatch.h:50
void unregisterIntRadPickup(Compute *cid, Box< Patch, Real > **const box)
Definition: Patch.C:182
Box< Patch, Real > * registerDHdrPrefixPickup(Compute *cid)
Definition: Patch.C:218
int doVirial
Definition: PatchTypes.h:21
static const LJTable * ljTable
Box< Patch, CompAtom > * velocityBox
BigReal pairlistTolerance
Definition: PatchTypes.h:40
BigReal gbis_delta
void unregisterVelocityPickup(Compute *cid, Box< Patch, CompAtom > **const box)
Definition: Patch.C:154
BigReal dielectric
int doGBIS
Definition: PatchTypes.h:28
Pairlists gbisStepPairlists[numGBISPairlists]
void submit(void)
Definition: ReductionMgr.h:323
Force * fullf[2]
int maxForceMerged
Definition: PatchTypes.h:32
BigReal * pressureProfileReduction
Data * open(void)
Definition: Box.h:39
static const int numGBISPairlists
BigReal maxGroupRadius
Definition: PatchTypes.h:42
const ComputeID cid
Definition: Compute.h:43
void unregisterDEdaSumDeposit(Compute *cid, Box< Patch, GBReal > **const box)
Definition: Patch.C:212
BigReal reductionData[reductionDataSize]
Box< Patch, CompAtom > * registerVelocityPickup(Compute *cid)
Definition: Patch.C:148
void close(Data **const t)
Definition: Box.h:49
int doMolly
Definition: PatchTypes.h:24
Vector c() const
Definition: Lattice.h:254
ComputeNonbondedWorkArrays * workArrays
static void(* calcMergeSelf)(nonbonded *)
double BigReal
Definition: common.h:114
int step
Definition: PatchTypes.h:16
static void(* calcFullSelfEnergy)(nonbonded *)