NAMD
ComputeNonbondedBase.h
Go to the documentation of this file.
1 
7 // Several special cases are defined:
8 // NBTYPE: exclusion method (NBPAIR, NBSELF -- mutually exclusive)
9 // FULLELECT full electrostatics calculation?
10 
11 // XXX Remove support for dead architectures and unused features:
12 //
13 // 1. SLOWONLY to generate _slow kernels is used only by Molly
14 // (mollified impulse method) so removal of Molly allows removal
15 // of this supporting machinery.
16 // 2. LES (locally enhanced sampling) is no longer in use;
17 // feature overlap with FEP.
18 
19 #ifdef ARCH_POWERPC
20 #include <builtins.h>
21 #endif
22 
23 #if defined(__SSE2__) && ! defined(NAMD_DISABLE_SSE)
24 #include <emmintrin.h> // We're using SSE2 intrinsics
25 #if defined(__INTEL_COMPILER)
26 #define __align(X) __declspec(align(X) )
27 #elif defined(__GNUC__) || defined(__PGI)
28 #define __align(X) __attribute__((aligned(X) ))
29 #else
30 #define __align(X) __declspec(align(X) )
31 #endif
32 #endif
33 
34 #if defined(__GNUC__)
35 #define UNLIKELY(X) __builtin_expect((X), 0)
36 #else
37 #define UNLIKELY(X) (X)
38 #endif
39 
40 #ifdef DEFINITION // (
41  #include "LJTable.h"
42  #include "Molecule.h"
43  #include "ComputeNonbondedUtil.h"
44 #endif // )
45  #include "Parameters.h"
46 #if NAMD_ComputeNonbonded_SortAtoms != 0
47  #include "PatchMap.h"
48 #endif
49 
50 // determining class name
51 #undef NAME
52 #undef CLASS
53 #undef CLASSNAME
54 #define NAME CLASSNAME(calc)
55 
56 #undef PAIR
57 #if NBTYPE == NBPAIR
58  #define PAIR(X) X
59  #define CLASS ComputeNonbondedPair
60  #define CLASSNAME(X) ENERGYNAME( X ## _pair )
61 #else
62  #define PAIR(X)
63 #endif
64 
65 #undef SELF
66 #if NBTYPE == NBSELF
67  #define SELF(X) X
68  #define CLASS ComputeNonbondedSelf
69  #define CLASSNAME(X) ENERGYNAME( X ## _self )
70 #else
71  #define SELF(X)
72 #endif
73 
74 #undef ENERGYNAME
75 #undef ENERGY
76 #undef NOENERGY
77 #ifdef CALCENERGY
78  #define ENERGY(X) X
79  #define NOENERGY(X)
80  #define ENERGYNAME(X) SLOWONLYNAME( X ## _energy )
81 #else
82  #define ENERGY(X)
83  #define NOENERGY(X) X
84  #define ENERGYNAME(X) SLOWONLYNAME( X )
85 #endif
86 
87 #undef SLOWONLYNAME
88 #undef FAST
89 #undef NOFAST
90 #ifdef SLOWONLY
91  #define FAST(X)
92  #define NOFAST(X) X
93  #define SLOWONLYNAME(X) MERGEELECTNAME( X ## _slow )
94 #else
95  #define FAST(X) X
96  #define NOFAST(X)
97  #define SLOWONLYNAME(X) MERGEELECTNAME( X )
98 #endif
99 
100 #undef MERGEELECTNAME
101 #undef SHORT
102 #undef NOSHORT
103 #ifdef MERGEELECT
104  #define SHORT(X)
105  #define NOSHORT(X) X
106  #define MERGEELECTNAME(X) FULLELECTNAME( X ## _merge )
107 #else
108  #define SHORT(X) X
109  #define NOSHORT(X)
110  #define MERGEELECTNAME(X) FULLELECTNAME( X )
111 #endif
112 
113 #undef FULLELECTNAME
114 #undef FULL
115 #undef NOFULL
116 #ifdef FULLELECT
117  #define FULLELECTNAME(X) TABENERGYNAME( X ## _fullelect )
118  #define FULL(X) X
119  #define NOFULL(X)
120 #else
121  #define FULLELECTNAME(X) TABENERGYNAME( X )
122  #define FULL(X)
123  #define NOFULL(X) X
124 #endif
125 
126 #undef TABENERGYNAME
127 #undef TABENERGY
128 #undef NOTABENERGY
129 #ifdef TABENERGYFLAG
130  #define TABENERGYNAME(X) FEPNAME( X ## _tabener )
131  #define TABENERGY(X) X
132  #define NOTABENERGY(X)
133 #else
134  #define TABENERGYNAME(X) FEPNAME( X )
135  #define TABENERGY(X)
136  #define NOTABENERGY(X) X
137 #endif
138 
139 // Here are what these macros stand for:
140 // FEP/NOT_FEP: FEP free energy perturbation is active/inactive
141 // (does NOT use LAM)
142 // LES: locally-enhanced sampling is active
143 // LAM: scale the potential by a factor lambda (except FEP)
144 // INT: measure interaction energies
145 // PPROF: pressure profiling
146 
147 #undef FEPNAME
148 #undef FEP
149 #undef LES
150 #undef INT
151 #undef PPROF
152 #undef LAM
153 #undef ALCH
154 #undef TI
155 #undef GO
156 #define FEPNAME(X) LAST( X )
157 #define FEP(X)
158 #define ALCHPAIR(X)
159 #define NOT_ALCHPAIR(X) X
160 #define LES(X)
161 #define INT(X)
162 #define PPROF(X)
163 #define LAM(X)
164 #define ALCH(X)
165 #define TI(X)
166 #define GO(X)
167 #ifdef FEPFLAG
168  #undef FEPNAME
169  #undef FEP
170  #undef ALCH
171  #define FEPNAME(X) LAST( X ## _fep )
172  #define FEP(X) X
173  #define ALCH(X) X
174 #endif
175 #ifdef TIFLAG
176  #undef FEPNAME
177  #undef TI
178  #undef ALCH
179  #define FEPNAME(X) LAST( X ## _ti )
180  #define TI(X) X
181  #define ALCH(X) X
182 #endif
183 #ifdef LESFLAG
184  #undef FEPNAME
185  #undef LES
186  #undef LAM
187  #define FEPNAME(X) LAST( X ## _les )
188  #define LES(X) X
189  #define LAM(X) X
190 #endif
191 #ifdef INTFLAG
192  #undef FEPNAME
193  #undef INT
194  #define FEPNAME(X) LAST( X ## _int )
195  #define INT(X) X
196 #endif
197 #ifdef PPROFFLAG
198  #undef FEPNAME
199  #undef INT
200  #undef PPROF
201  #define FEPNAME(X) LAST( X ## _pprof )
202  #define INT(X) X
203  #define PPROF(X) X
204 #endif
205 #ifdef GOFORCES
206  #undef FEPNAME
207  #undef GO
208  #define FEPNAME(X) LAST( X ## _go )
209  #define GO(X) X
210 #endif
211 
212 #define LAST(X) X
213 
214 // see if things are really messed up
215 SELF( PAIR( foo bar ) )
216 LES( FEP( foo bar ) )
217 LES( INT( foo bar ) )
218 FEP( INT( foo bar ) )
219 LAM( INT( foo bar ) )
220 FEP( NOENERGY( foo bar ) )
221 ENERGY( NOENERGY( foo bar ) )
222 TABENERGY(NOTABENERGY( foo bar ) )
223 
224 #define KNL_MAKE_DEPENDS_INCLUDE
225 #include "ComputeNonbondedBase2KNL.h"
226 #undef KNL_MAKE_DEPENDS_INCLUDE
227 
228 #undef KNL
229 #undef NOKNL
230 #ifdef NAMD_KNL
231  #if ( TABENERGY(1+) FEP(1+) TI(1+) INT(1+) LES(1+) GO(1+) PPROF(1+) NOFAST(1+) 0 )
232  #define KNL(X)
233  #define NOKNL(X) X
234  #else
235  #define KNL(X) X
236  #define NOKNL(X)
237  #endif
238 #else
239  #define KNL(X)
240  #define NOKNL(X) X
241 #endif
242 
243 #if NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR( + 1 ) )
244  #define COMPONENT_DOTPRODUCT(A,B) ((A##_x * B##_x) + (A##_y * B##_y) + (A##_z * B##_z))
245 #endif
246 
247 
248 // ************************************************************
249 // function header
250 void ComputeNonbondedUtil :: NAME
251  ( nonbonded *params )
252 
253 // function body
254 {
255  // int NAME; // to track errors via unused variable warnings
256 
257  if ( ComputeNonbondedUtil::commOnly ) return;
258 
259  // speedup variables
260  BigReal *reduction = params->reduction;
261  SimParameters *simParams = params->simParameters;
262 
263  PPROF(
264  BigReal *pressureProfileReduction = params->pressureProfileReduction;
265  const BigReal invThickness = 1.0 / pressureProfileThickness;
266  )
267 
268  Pairlists &pairlists = *(params->pairlists);
269  int savePairlists = params->savePairlists;
270  int usePairlists = params->usePairlists;
271  pairlists.reset();
272  // PAIR(iout << "--------\n" << endi;)
273 
274  // BEGIN LA
275  const int doLoweAndersen = params->doLoweAndersen;
276  // END LA
277 
278  // local variables
279  int exclChecksum = 0;
280  FAST
281  (
282  // JLai
283  ENERGY( BigReal vdwEnergy = 0;
284  GO( BigReal groLJEnergy = 0;
285  BigReal groGaussEnergy = 0;
286  BigReal goEnergyNative = 0;
287  BigReal goEnergyNonnative = 0; ) )
288  SHORT
289  (
290  ENERGY( BigReal electEnergy = 0; )
291  Vector f_net = 0.;
292  )
293 
294  FEP
295  (
296  ENERGY( BigReal vdwEnergy_s = 0; )
297  SHORT
298  (
299  ENERGY( BigReal electEnergy_s = 0; )
300  )
301  )
302  )
303 
304  FULL
305  (
306  ENERGY( BigReal fullElectEnergy = 0; )
307  Vector fullf_net = 0.;
308  FEP
309  (
310  ENERGY( BigReal fullElectEnergy_s = 0; )
311  )
312  )
313 
314  // Bringing stuff into local namespace for speed.
315  const BigReal offset_x = params->offset.x;
316  const BigReal offset_y = params->offset.y;
317  const BigReal offset_z = params->offset.z;
318 
319  // Note that offset_f assumes positions are relative to patch centers
320  const float offset_x_f = params->offset_f.x;
321  const float offset_y_f = params->offset_f.y;
322  const float offset_z_f = params->offset_f.z;
323 
324  register const BigReal plcutoff2 = \
325  params->plcutoff * params->plcutoff;
326  register const BigReal groupplcutoff2 = \
327  params->groupplcutoff * params->groupplcutoff;
328  const BigReal dielectric_1 = ComputeNonbondedUtil:: dielectric_1;
329  const LJTable* const ljTable = ComputeNonbondedUtil:: ljTable;
330  LJTable::TableEntry ljNull; ljNull.A = 0; ljNull.B = 0;
331  const LJTable::TableEntry* const lj_null_pars = &ljNull;
332  const Molecule* const mol = ComputeNonbondedUtil:: mol;
333  SHORT
334  (
335  const BigReal* const table_four = ComputeNonbondedUtil:: table_short;
336  )
337  FULL
338  (
339  SHORT
340  (
341  const BigReal* const slow_table = ComputeNonbondedUtil:: slow_table;
342  )
343  NOSHORT
344  (
345 //#if 1 ALCH(-1)
346  const BigReal* const table_four = ComputeNonbondedUtil:: table_noshort;
347 //#else // have to switch this for ALCH
348 // BigReal* table_four = ComputeNonbondedUtil:: table_noshort;
349 //#endif
350  )
351  )
352  BigReal scaling = ComputeNonbondedUtil:: scaling;
353  const float scaling_f = scaling;
354  const BigReal modf_mod = 1.0 - scale14;
355  FAST
356  (
357  const BigReal switchOn2 = ComputeNonbondedUtil:: switchOn2;
358  const float switchOn2_f = switchOn2;
359  const float cutoff2_f = ComputeNonbondedUtil::cutoff2;
360  const BigReal c1 = ComputeNonbondedUtil:: c1;
361  const BigReal c3 = ComputeNonbondedUtil:: c3;
362  const float c1_f = c1;
363  const float c3_f = c3;
364  )
365  const BigReal r2_delta = ComputeNonbondedUtil:: r2_delta;
366  const int r2_delta_exp = ComputeNonbondedUtil:: r2_delta_exp;
367  // const int r2_delta_expc = 64 * (r2_delta_exp - 127);
368  const int r2_delta_expc = 64 * (r2_delta_exp - 1023);
369 
370  ALCH(
371  const BigReal switchdist2 = ComputeNonbondedUtil::switchOn2;
372  const BigReal cutoff2 = ComputeNonbondedUtil::cutoff2;
373  //JM: switchfactor is set to inf if the switching flag is off
374  const BigReal diff = cutoff2 - switchdist2;
375  const BigReal switchfactor = 1./(diff*diff*diff);
376  const BigReal alchVdwShiftCoeff = ComputeNonbondedUtil::alchVdwShiftCoeff;
377  const Bool alchWCAOn = ComputeNonbondedUtil::alchWCAOn;
378  const Bool vdwForceSwitching = ComputeNonbondedUtil::vdwForceSwitching;
379 
380  /* BKR - Enable dynamic lambda by passing the timestep to simParams.
381  The SimParameters object now has all knowledge of the lambda schedule
382  and all operations should ask _it_ for lambda, rather than recomputing
383  it.
384  */
385  const BigReal alchLambda = simParams->getCurrentLambda(params->step);
386 
387  /*lambda values 'up' are for atoms scaled up with lambda (partition 1)*/
388  BigReal lambdaUp = alchLambda; // needed in ComputeNonbondedBase2.h!
389  BigReal elecLambdaUp = simParams->getElecLambda(lambdaUp);
390  BigReal vdwLambdaUp = simParams->getVdwLambda(lambdaUp);
391  BigReal repLambdaUp = simParams->getRepLambda(lambdaUp);
392  BigReal vdwShiftUp = alchVdwShiftCoeff*(1 - vdwLambdaUp);
393  FEP(BigReal lambda2Up = simParams->getCurrentLambda2(params->step);)
394  FEP(BigReal elecLambda2Up = simParams->getElecLambda(lambda2Up);)
395  FEP(BigReal vdwLambda2Up = simParams->getVdwLambda(lambda2Up);)
396  FEP(BigReal repLambda2Up = simParams->getRepLambda(lambda2Up);)
397  FEP(BigReal vdwShift2Up = alchVdwShiftCoeff*(1 - vdwLambda2Up);)
398 
399  /*lambda values 'down' are for atoms scaled down with lambda (partition 2)*/
400  BigReal lambdaDown = 1 - alchLambda; // needed in ComputeNonbondedBase2.h!
401  BigReal elecLambdaDown = simParams->getElecLambda(lambdaDown);
402  BigReal vdwLambdaDown = simParams->getVdwLambda(lambdaDown);
403  BigReal repLambdaDown = simParams->getRepLambda(lambdaDown);
404  BigReal vdwShiftDown = alchVdwShiftCoeff*(1 - vdwLambdaDown);
405  FEP(BigReal lambda2Down = 1 - simParams->getCurrentLambda2(params->step);)
406  FEP(BigReal elecLambda2Down = simParams->getElecLambda(lambda2Down);)
407  FEP(BigReal vdwLambda2Down = simParams->getVdwLambda(lambda2Down);)
408  FEP(BigReal repLambda2Down = simParams->getRepLambda(lambda2Down);)
409  FEP(BigReal vdwShift2Down = alchVdwShiftCoeff*(1 - vdwLambda2Down);)
410 
411  // Thermodynamic Integration Notes:
412  // Separation of atom pairs into different pairlist according to lambda
413  // coupling, for alchemical free energy calculations. Normal pairlists
414  // (pairlist[nxm]_save) are re-used for non-lambda-coupled pairs; new ones
415  // (pairlist[nxm][12]_save are created for atoms switched up or down with
416  // lambda respectively.
417  // This makes TI coding far easier and more readable, since it's necessary
418  // to store dU/dlambda in different variables depending on whether atoms are
419  // being switched up or down. Further, it allows the prior explicit coding of
420  // the separation-shifted vdW potential to stay in place without a big
421  // performance hit, since it's no longer necessary to calculate vdW potentials
422  // explicity for the bulk (non-alchemical) interactions - so that part of the
423  // free energy code stays readable and easy to modify. Finally there should
424  // be some performance gain over the old FEP implementation as we no
425  // longer have to check the partitions of each atom pair and switch
426  // parameters accordingly for every single NonbondedBase2.h loop - this is
427  // done at the pairlist level
428 
429  // XXX pswitchTable is repeatedly calculated each time
430  // since its value is fixed, it should be calculated once up front
431  int pswitchTable[5*5];
432  // determines which pairlist alchemical pairings get sent to
433  // 0: non-alchemical pairs, partition 0 <-> partition 0
434  // 1 and 3: atoms scaled up as lambda increases, p0<->p1, p0<->p3
435  // 2 and 4: atoms scaled down as lambda increases, p0<->p2, p0<->p4
436  // all p1<->p2, p1<->p4, p3<->p2, p3<->p4 interactions to be dropped (99)
437  // in general, 'up' refers to 1 and 3, 'down' refers to 2 and 4
438  for (int ip=0; ip<5; ++ip) {
439  for (int jp=0; jp<5; ++jp ) {
440  pswitchTable[ip+5*jp] = 0;
441  if ((ip==1 && jp==0) || (ip==0 && jp==1)) pswitchTable[ip+5*jp] = 1;
442  if ((ip==2 && jp==0) || (ip==0 && jp==2)) pswitchTable[ip+5*jp] = 2;
443  if ((ip==3 && jp==0) || (ip==0 && jp==3)) pswitchTable[ip+5*jp] = 3;
444  if ((ip==4 && jp==0) || (ip==0 && jp==4)) pswitchTable[ip+5*jp] = 4;
445 
446  if (ip && jp && (abs(ip - jp) != 2)) pswitchTable[ip+5*jp] = 99; // no interaction
447  if (! ComputeNonbondedUtil::alchDecouple) {
448  // no decoupling: interactions within an alchemical moiety are treated like
449  // normal alchemical pairs
450  if ((ip == 1 && jp == 1) || (ip == 1 && jp == 3) || (ip == 3 && jp == 1)) pswitchTable[ip+5*jp] = 1;
451  if ((ip == 2 && jp == 2) || (ip == 2 && jp == 4) || (ip == 4 && jp == 2)) pswitchTable[ip+5*jp] = 2;
452  if (ip == 3 && jp == 3) pswitchTable[ip+5*jp] = 3;
453  if (ip == 4 && jp == 4) pswitchTable[ip+5*jp] = 4;
454  }
455  if (ComputeNonbondedUtil::alchDecouple) {
456  // decoupling: PME calculates extra grids so that while PME
457  // interaction with the full system is switched off, a new PME grid
458  // containing only alchemical atoms is switched on. Full interactions
459  // between alchemical atoms are maintained; potentials within one
460  // partition need not be scaled here.
461  if (ip == 1 && jp == 1) pswitchTable[ip+5*jp] = 0; //Single topology only supports alchDecouple off!
462  if (ip == 2 && jp == 2) pswitchTable[ip+5*jp] = 0; //So 'alchDecouple on' only works for dual topology!
463  }
464  }
465  }
466 
467  // dU/dlambda variables for thermodynamic integration
468  TI(
469  BigReal vdwEnergy_ti_1 = 0;
470  BigReal vdwEnergy_ti_2 = 0;
471  SHORT(BigReal electEnergy_ti_1 = 0;
472  BigReal electEnergy_ti_2 = 0;)
473  FULL(BigReal fullElectEnergy_ti_1 = 0;
474  BigReal fullElectEnergy_ti_2 = 0;)
475  )
476  )
477 
478 
479  const int i_upper = params->numAtoms[0];
480  register const int j_upper = params->numAtoms[1];
481  register int j;
482  register int i;
483  const CompAtom *p_0 = params->p[0];
484  const CompAtom *p_1 = params->p[1];
485  KNL( const CompAtomFlt *pFlt_0 = params->pFlt[0]; )
486  KNL( const CompAtomFlt *pFlt_1 = params->pFlt[1]; )
487  const CompAtomExt *pExt_0 = params->pExt[0];
488  const CompAtomExt *pExt_1 = params->pExt[1];
489 
490  char * excl_flags_buff = 0;
491  const int32 * full_excl = 0;
492  const int32 * mod_excl = 0;
493 
494  plint *pairlistn_save; int npairn;
495  plint *pairlistx_save; int npairx;
496  plint *pairlistm_save; int npairm;
497 
498  ALCH(
499  // n = normal, x = excluded, m = modified
500  // A[1-4] = alchemical partition 1-4
501  plint *pairlistnA1_save; int npairnA1;
502  plint *pairlistxA1_save; int npairxA1;
503  plint *pairlistmA1_save; int npairmA1;
504  plint *pairlistnA2_save; int npairnA2;
505  plint *pairlistxA2_save; int npairxA2;
506  plint *pairlistmA2_save; int npairmA2;
507  plint *pairlistnA3_save; int npairnA3;
508  plint *pairlistxA3_save; int npairxA3;
509  plint *pairlistmA3_save; int npairmA3;
510  plint *pairlistnA4_save; int npairnA4;
511  plint *pairlistxA4_save; int npairxA4;
512  plint *pairlistmA4_save; int npairmA4;
513  )
514 
515  NBWORKARRAYSINIT(params->workArrays);
516 
517  int arraysize = j_upper+5;
518 
519  NBWORKARRAY(int,pairlisti,arraysize)
520  NBWORKARRAY(BigReal,r2list,arraysize)
521  KNL( NBWORKARRAY(float,r2list_f,arraysize) )
522  KNL( NBWORKARRAY(float,xlist,arraysize) )
523  KNL( NBWORKARRAY(float,ylist,arraysize) )
524  KNL( NBWORKARRAY(float,zlist,arraysize) )
525 
526  union { double f; int32 i[2]; } byte_order_test;
527  byte_order_test.f = 1.0; // should occupy high-order bits only
528  int32 *r2iilist = (int32*)r2list + ( byte_order_test.i[0] ? 0 : 1 );
529 
530  if ( ! ( savePairlists || ! usePairlists ) ) arraysize = 0;
531 
532 
533  // DMK - Atom Sort
534  //
535  // Basic Idea: Determine the line between the center of masses of
536  // the two patches. Project and then sort the lists of atoms
537  // along this line. Then, as the pairlists are being generated
538  // for the atoms in the first atom list, use the sorted
539  // list to only add atoms from the second list that are between
540  // +/- ~cutoff from the atoms position on the line.
541  #if NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR( + 1 ) )
542 
543  // NOTE: For the first atom list, just the sort values themselves need to be
544  // calculated (a BigReal vs. SortEntry data structure). However, a second
545  // SortEntry array is needed to perform the merge sort on the second list of
546  // atoms. Use the atomSort_[0/1]_sortValues__ arrays to sort the second
547  // list of atoms and then use the left over array to make a version of the
548  // list that only includes non-fixed groups (fixg).
549 
550  NBWORKARRAY(SortEntry, atomSort_0_sortValues__, arraysize);
551  NBWORKARRAY(SortEntry, atomSort_1_sortValues__, arraysize);
552  NBWORKARRAY(BigReal, p_0_sortValues, arraysize);
553 
554  register SortEntry* p_1_sortValues = atomSort_0_sortValues__;
555  register SortEntry* p_1_sortValues_fixg = atomSort_1_sortValues__;
556 
557  int p_0_sortValues_len = 0;
558  int p_1_sortValues_len = 0;
559  int p_1_sortValues_fixg_len = 0;
560 
561  // Calculate the distance between to projected points on the line that
562  // represents the cutoff distance.
563  BigReal atomSort_windowRadius = sqrt(groupplcutoff2);
564 
565  if (savePairlists || !usePairlists) {
566 
567  register const BigReal projLineVec_x = params->projLineVec.x;
568  register const BigReal projLineVec_y = params->projLineVec.y;
569  register const BigReal projLineVec_z = params->projLineVec.z;
570 
571  // Calculate the sort values for the atoms in patch 1
572  {
573  register int nbgs = p_1->nonbondedGroupSize;
574  register BigReal p_x = p_1->position.x;
575  register BigReal p_y = p_1->position.y;
576  register BigReal p_z = p_1->position.z;
577  register int index = 0;
578 
579  for (register int j = nbgs; j < j_upper; j += nbgs) {
580 
581  // Set p_j_next to point to the atom for the next iteration and begin
582  // loading the 'nbgs' value for that atom.
583  register const CompAtom* p_j_next = p_1 + j;
584  nbgs = p_j_next->nonbondedGroupSize;
585 
586  // Calculate the distance along the projection vector
587  // NOTE: If the vector from the origin to the point is 'A' and the vector
588  // between the patches is 'B' then to project 'A' onto 'B' we take the dot
589  // product of the two vectors 'A dot B' divided by the length of 'B'.
590  // So... projection of A onto B = (A dot B) / length(B), however, note
591  // that length(B) == 1, therefore the sort value is simply (A dot B)
592  register BigReal sortVal = COMPONENT_DOTPRODUCT(p,projLineVec);
593 
594  // Start loading the next iteration's atom's position
595  p_x = p_j_next->position.x;
596  p_y = p_j_next->position.y;
597  p_z = p_j_next->position.z;
598 
599  // Store the caclulated sort value into the array of sort values
600  register SortEntry* p_1_sortValStorePtr = p_1_sortValues + p_1_sortValues_len;
601  p_1_sortValStorePtr->index = index;
602  p_1_sortValStorePtr->sortValue = sortVal;
603  p_1_sortValues_len++;
604 
605  // Update index for the next iteration
606  index = j;
607 
608  } // end while (j < j_upper)
609 
610  register BigReal sortVal = COMPONENT_DOTPRODUCT(p,projLineVec);
611 
612  register SortEntry* p_1_sortValStorePtr = p_1_sortValues + p_1_sortValues_len;
613  p_1_sortValStorePtr->index = index;
614  p_1_sortValStorePtr->sortValue = sortVal;
615  p_1_sortValues_len++;
616  }
617 
618  // NOTE: This list and another version of it with only non-fixed
619  // atoms will be used in place of grouplist and fixglist.
620  #if 0 // Selection Sort
621  sortEntries_selectionSort(p_1_sortValues, p_1_sortValues_len);
622  #elif 0 // Bubble Sort
623  sortEntries_bubbleSort(p_1_sortValues, p_1_sortValues_len);
624  #else // Merge Sort
625  #if NAMD_ComputeNonbonded_SortAtoms_LessBranches == 0
626  sortEntries_mergeSort_v1(p_1_sortValues, p_1_sortValues_fixg, p_1_sortValues_len);
627  #else
628  sortEntries_mergeSort_v2(p_1_sortValues, p_1_sortValues_fixg, p_1_sortValues_len);
629  #endif
630  #endif
631 
632  // Calculate the sort values for the atoms in patch 0
633  {
634  register int nbgs = p_0->nonbondedGroupSize;
635  register BigReal p_x = p_0->position.x + offset_x;
636  register BigReal p_y = p_0->position.y + offset_y;
637  register BigReal p_z = p_0->position.z + offset_z;
638  register int index = 0;
639 
640  for (register int i = nbgs; i < i_upper; i += nbgs) {
641 
642  // Set p_i_next to point to the atom for the next iteration and begin
643  // loading the 'nbgs' value for that atom.
644  register const CompAtom* p_i_next = p_0 + i;
645  nbgs = p_i_next->nonbondedGroupSize;
646 
647  // Calculate the distance along the projection vector
648  register BigReal sortVal = COMPONENT_DOTPRODUCT(p,projLineVec);
649 
650  // Start loading the next iteration's atom's position
651  p_x = p_i_next->position.x + offset_x;
652  p_y = p_i_next->position.y + offset_y;
653  p_z = p_i_next->position.z + offset_z;
654 
655  // Store the calculated sort value into the array of sort values
656  register BigReal* p_0_sortValStorePtr = p_0_sortValues + index;
657  *p_0_sortValStorePtr = sortVal;
658 
659  // Update index for the next iteration
660  index = i;
661  }
662 
663  register BigReal sortVal = COMPONENT_DOTPRODUCT(p,projLineVec);
664 
665  register BigReal* p_0_sortValStorePtr = p_0_sortValues + index;
666  *p_0_sortValStorePtr = sortVal;
667 
668  p_0_sortValues_len = i_upper;
669  }
670 
671  } // end if (savePairlists || !usePairlists)
672 
673  #endif // NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR( + 1 ) )
674 
675  // Atom Sort : The grouplist and fixglist arrays are not needed when the
676  // the atom sorting code is in use.
677  #if ! (NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR( + 1 ) ) )
678  NBWORKARRAY(int,grouplist,arraysize);
679  NBWORKARRAY(int,fixglist,arraysize);
680  #endif
681 
682  NBWORKARRAY(int,goodglist,arraysize);
683  NBWORKARRAY(plint,pairlistx,arraysize);
684  NBWORKARRAY(plint,pairlistm,arraysize);
685  NBWORKARRAY(int,pairlist,arraysize);
686  NBWORKARRAY(int,pairlist2,arraysize);
687  ALCH(
688  NBWORKARRAY(plint,pairlistnA1,arraysize);
689  NBWORKARRAY(plint,pairlistxA1,arraysize);
690  NBWORKARRAY(plint,pairlistmA1,arraysize);
691  NBWORKARRAY(plint,pairlistnA2,arraysize);
692  NBWORKARRAY(plint,pairlistxA2,arraysize);
693  NBWORKARRAY(plint,pairlistmA2,arraysize);
694  NBWORKARRAY(plint,pairlistnA3,arraysize);
695  NBWORKARRAY(plint,pairlistxA3,arraysize);
696  NBWORKARRAY(plint,pairlistmA3,arraysize);
697  NBWORKARRAY(plint,pairlistnA4,arraysize);
698  NBWORKARRAY(plint,pairlistxA4,arraysize);
699  NBWORKARRAY(plint,pairlistmA4,arraysize);
700  )
701 
702  int fixg_upper = 0;
703  int g_upper = 0;
704 
705  if ( savePairlists || ! usePairlists ) {
706 
707  #if NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR ( + 1 ) )
708 
709  // Create a sorted list of non-fixed groups
710  register int fixg = 0;
711  for (int tmpI = 0; tmpI < p_1_sortValues_len; tmpI++) {
712  register SortEntry* p_1_sortEntry = p_1_sortValues + tmpI;
713  register int p_1_index = p_1_sortEntry->index;
714  if (!pExt_1[p_1_index].groupFixed) {
715  register SortEntry* p_1_sortEntry_fixg = p_1_sortValues_fixg + p_1_sortValues_fixg_len;
716  p_1_sortEntry_fixg->index = p_1_sortEntry->index;
717  p_1_sortEntry_fixg->sortValue = p_1_sortEntry->sortValue;
718  p_1_sortValues_fixg_len++;
719  }
720  }
721 
722  #else
723 
724  register int g = 0;
725  for ( j = 0; j < j_upper; ++j ) {
726  if ( p_1[j].nonbondedGroupSize ) {
727  grouplist[g++] = j;
728  }
729  }
730  g_upper = g;
731  if ( g_upper ) grouplist[g_upper] = grouplist[g_upper-1];
732  int fixg = 0;
733 
734  if ( fixedAtomsOn ) {
735  for ( g = 0; g < g_upper; ++g ) {
736  j = grouplist[g];
737  if ( ! pExt_1[j].groupFixed ) {
738  fixglist[fixg++] = j;
739  }
740  }
741  }
742 
743  fixg_upper = fixg;
744  if ( fixg_upper ) fixglist[fixg_upper] = fixglist[fixg_upper-1];
745 
746  #endif // NAMD_ComputeNonbonded_SortAtoms != 0
747 
748  pairlists.addIndex();
749  pairlists.setIndexValue(i_upper);
750 
751  } else { // if ( savePairlists || ! usePairlists )
752 
753  if ( pairlists.getIndexValue() != i_upper )
754  NAMD_bug("pairlist i_upper mismatch!");
755 
756  } // if ( savePairlists || ! usePairlists )
757 
758  SELF(
759  int j_hgroup = 0;
760  int g_lower = 0;
761  int fixg_lower = 0;
762  )
763  int pairlistindex=0;
764  int pairlistoffset=0;
765 
766 
767 
768 #if ( SHORT( FAST( 1+ ) ) 0 )
769  Force *f_0 = params->ff[0];
770 #if ( PAIR( 1+ ) 0 )
771  Force *f_1 = params->ff[1];
772 #else
773 #define f_1 f_0
774 #endif
775 #endif
776 #if ( FULL( 1+ ) 0 )
777  Force *fullf_0 = params->fullf[0];
778 #if ( PAIR( 1+ ) 0 )
779  Force *fullf_1 = params->fullf[1];
780 #else
781 #define fullf_1 fullf_0
782 #endif
783 #endif
784 
785 
786  int maxPart = params->numParts - 1;
787  int groupCount = params->minPart;
788  PAIR(
789  if ( savePairlists || ! usePairlists ) {
790  i = 0;
791  pairlists.addIndex();
792  } else {
793  i = pairlists.getIndexValue();
794  }
795  )
796  PAIR(for ( ; i < (i_upper);)) SELF(for ( i=0; i < (i_upper- 1);i++))
797  {
798  const CompAtom &p_i = p_0[i];
799  KNL( const CompAtomFlt &pFlt_i = pFlt_0[i]; )
800  const CompAtomExt &pExt_i = pExt_0[i];
801 
802  PAIR(if (savePairlists || ! usePairlists){)
803  if ( p_i.hydrogenGroupSize ) {
804  if ( groupCount ) { // skip this group
805  --groupCount;
806  i += p_i.hydrogenGroupSize;
807 #ifdef ARCH_POWERPC
808  __dcbt((void *) &(p_0[i]));
809 #endif
810  SELF(--i;)
811  continue;
812  } else { // compute this group
813  groupCount = maxPart;
814  }
815  }
816  PAIR(})
817 
818  register const BigReal p_i_x = p_i.position.x + offset_x;
819  register const BigReal p_i_y = p_i.position.y + offset_y;
820  register const BigReal p_i_z = p_i.position.z + offset_z;
821 #if KNL(1)+0
822  register const BigReal p_i_x_f = pFlt_i.position.x + offset_x_f;
823  register const BigReal p_i_y_f = pFlt_i.position.y + offset_y_f;
824  register const BigReal p_i_z_f = pFlt_i.position.z + offset_z_f;
825 #endif
826 
827  ALCH(const int p_i_partition = p_i.partition;)
828 
829  PPROF(
830  const int p_i_partition = p_i.partition;
831  int n1 = (int)floor((p_i.position.z-pressureProfileMin)*invThickness);
832  pp_clamp(n1, pressureProfileSlabs);
833  )
834 
835  SELF ( if ( p_i.hydrogenGroupSize ) j_hgroup = i + p_i.hydrogenGroupSize; )
836 
837  if ( savePairlists || ! usePairlists ) {
838 
839  #ifdef MEM_OPT_VERSION
840  const ExclusionCheck *exclcheck = mol->get_excl_check_for_idx(pExt_i.exclId);
841  const int excl_min = pExt_i.id + exclcheck->min;
842  const int excl_max = pExt_i.id + exclcheck->max;
843  #else
844  const ExclusionCheck *exclcheck = mol->get_excl_check_for_atom(pExt_i.id);
845  const int excl_min = exclcheck->min;
846  const int excl_max = exclcheck->max;
847  #endif
848  const char * excl_flags_var;
849  if ( exclcheck->flags ) excl_flags_var = exclcheck->flags - excl_min;
850  else { // need to build list on the fly
851 
852  //TODO: Should change later!!!!!!!!!! --Chao Mei
853  //Now just for the sake of passing compilation
854  #ifndef MEM_OPT_VERSION
855 #define REZERO_EXCL_FLAGS_BUFF if ( excl_flags_buff ) { \
856  int nl,l; \
857  nl = full_excl[0] + 1; \
858  for ( l=1; l<nl; ++l ) excl_flags_buff[full_excl[l]] = 0; \
859  nl = mod_excl[0] + 1; \
860  for ( l=1; l<nl; ++l ) excl_flags_buff[mod_excl[l]] = 0; \
861  }
862  REZERO_EXCL_FLAGS_BUFF else {
863  ResizeArray<char> &wa = computeNonbondedWorkArrays->excl_flags_buff;
864  if ( wa.size() < mol->numAtoms ) {
865  int oldsize = wa.size();
866  wa.resize(mol->numAtoms);
867  memset( (void*) &wa[oldsize], 0, wa.size() - oldsize);
868  }
869  excl_flags_buff = &wa[0];
870  }
871  int nl,l;
872  full_excl = mol->get_full_exclusions_for_atom(pExt_i.id);
873  nl = full_excl[0] + 1;
874  for ( l=1; l<nl; ++l ) excl_flags_buff[full_excl[l]] = EXCHCK_FULL;
875  mod_excl = mol->get_mod_exclusions_for_atom(pExt_i.id);
876  nl = mod_excl[0] + 1;
877  for ( l=1; l<nl; ++l ) excl_flags_buff[mod_excl[l]] = EXCHCK_MOD;
878  excl_flags_var = excl_flags_buff;
879  #endif
880 
881  }
882  const char * const excl_flags = excl_flags_var;
883 
884  if ( p_i.nonbondedGroupSize ) {
885 
886  pairlistindex = 0; // initialize with 0 elements
887  pairlistoffset=0;
888  const int groupfixed = ( fixedAtomsOn && pExt_i.groupFixed );
889 
890  // If patch divisions are not made by hydrogen groups, then
891  // hydrogenGroupSize is set to 1 for all atoms. Thus we can
892  // carry on as if we did have groups - only less efficiently.
893  // An optimization in this case is to not rebuild the temporary
894  // pairlist but to include every atom in it. This should be a
895  // a very minor expense.
896 
897  SELF
898  (
899  if ( p_i.hydrogenGroupSize ) {
900  // exclude child hydrogens of i
901  // j_hgroup = i + p_i.hydrogenGroupSize; (moved above)
902  while ( g_lower < g_upper &&
903  grouplist[g_lower] < j_hgroup ) ++g_lower;
904  while ( fixg_lower < fixg_upper &&
905  fixglist[fixg_lower] < j_hgroup ) ++fixg_lower;
906  }
907  // add all child or sister hydrogens of i
908  for ( j = i + 1; j < j_hgroup; ++j ) {
909  pairlist[pairlistindex++] = j;
910  }
911  )
912 
913  // add remaining atoms to pairlist via hydrogen groups
914  register int *pli = pairlist + pairlistindex;
915 
916  {
917  // Atom Sort : Modify the values of g and gu based on the added information
918  // of the linear projections (sort values) information.
919  #if NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR ( + 1 ) )
920 
921  register int g = 0;
922  register BigReal p_i_sortValue = p_0_sortValues[i];
923  const BigReal p_i_sortValue_plus_windowRadius = p_i_sortValue + atomSort_windowRadius;
924  register SortEntry* sortValues = ( groupfixed ? p_1_sortValues_fixg : p_1_sortValues );
925 
926  // Find the actual gu (upper bound in sorted list for this outer-loop atom) based on the sort values
927  register int lower = 0;
928  register int upper = (groupfixed ? p_1_sortValues_fixg_len : p_1_sortValues_len);
929  while ((upper - lower) > 1) {
930  register int j = ((lower + upper) >> 1);
931  register BigReal jSortVal = sortValues[j].sortValue;
932  if (jSortVal < p_i_sortValue_plus_windowRadius) {
933  lower = j;
934  } else {
935  upper = j;
936  }
937  }
938  const int gu = (sortValues[lower].sortValue >= p_i_sortValue_plus_windowRadius) ? lower : upper;
939 
940  #else
941 
942  register int *gli = goodglist;
943  const int *glist = ( groupfixed ? fixglist : grouplist );
944  SELF( const int gl = ( groupfixed ? fixg_lower : g_lower ); )
945  const int gu = ( groupfixed ? fixg_upper : g_upper );
946  register int g = PAIR(0) SELF(gl);
947 
948  #endif
949 
950 
951  if ( g < gu ) {
952  int hu = 0;
953 #ifndef NAMD_KNL
954 #if defined(__SSE2__) && ! defined(NAMD_DISABLE_SSE)
955  if ( gu - g > 6 ) {
956 
957  #if NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR ( + 1 ) )
958  register SortEntry* sortEntry0 = sortValues + g;
959  register SortEntry* sortEntry1 = sortValues + g + 1;
960  register int jprev0 = sortEntry0->index;
961  register int jprev1 = sortEntry1->index;
962  #else
963  register int jprev0 = glist[g ];
964  register int jprev1 = glist[g + 1];
965  #endif
966 
967  register int j0;
968  register int j1;
969 
970  __m128d PJ_X_01 = _mm_set_pd(p_1[jprev1].position.x, p_1[jprev0].position.x);
971  __m128d PJ_Y_01 = _mm_set_pd(p_1[jprev1].position.y, p_1[jprev0].position.y);
972  __m128d PJ_Z_01 = _mm_set_pd(p_1[jprev1].position.z, p_1[jprev0].position.z);
973 
974  // these don't change here, so we could move them into outer scope
975  const __m128d P_I_X = _mm_set1_pd(p_i_x);
976  const __m128d P_I_Y = _mm_set1_pd(p_i_y);
977  const __m128d P_I_Z = _mm_set1_pd(p_i_z);
978 
979  g += 2;
980  for ( ; g < gu - 2; g +=2 ) {
981  // compute 1d distance, 2-way parallel
982  j0 = jprev0;
983  j1 = jprev1;
984 
985  __m128d T_01 = _mm_sub_pd(P_I_X, PJ_X_01);
986  __m128d R2_01 = _mm_mul_pd(T_01, T_01);
987  T_01 = _mm_sub_pd(P_I_Y, PJ_Y_01);
988  R2_01 = _mm_add_pd(R2_01, _mm_mul_pd(T_01, T_01));
989  T_01 = _mm_sub_pd(P_I_Z, PJ_Z_01);
990  R2_01 = _mm_add_pd(R2_01, _mm_mul_pd(T_01, T_01));
991 
992  #if NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR ( + 1 ) )
993  sortEntry0 = sortValues + g;
994  sortEntry1 = sortValues + g + 1;
995  jprev0 = sortEntry0->index;
996  jprev1 = sortEntry1->index;
997  #else
998  jprev0 = glist[g ];
999  jprev1 = glist[g+1];
1000  #endif
1001 
1002  PJ_X_01 = _mm_set_pd(p_1[jprev1].position.x, p_1[jprev0].position.x);
1003  PJ_Y_01 = _mm_set_pd(p_1[jprev1].position.y, p_1[jprev0].position.y);
1004  PJ_Z_01 = _mm_set_pd(p_1[jprev1].position.z, p_1[jprev0].position.z);
1005 
1006  __align(16) double r2_01[2];
1007  _mm_store_pd(r2_01, R2_01); // 16-byte-aligned store
1008 
1009  // XXX these could possibly benefit from SSE-based conditionals
1010  bool test0 = ( r2_01[0] < groupplcutoff2 );
1011  bool test1 = ( r2_01[1] < groupplcutoff2 );
1012 
1013  //removing ifs benefits on many architectures
1014  //as the extra stores will only warm the cache up
1015  goodglist [ hu ] = j0;
1016  goodglist [ hu + test0 ] = j1;
1017 
1018  hu += test0 + test1;
1019  }
1020  g-=2;
1021  }
1022 #else
1023  if ( gu - g > 6 ) {
1024 
1025  #if NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR ( + 1 ) )
1026  register SortEntry* sortEntry0 = sortValues + g;
1027  register SortEntry* sortEntry1 = sortValues + g + 1;
1028  register int jprev0 = sortEntry0->index;
1029  register int jprev1 = sortEntry1->index;
1030  #else
1031  register int jprev0 = glist[g ];
1032  register int jprev1 = glist[g + 1];
1033  #endif
1034 
1035  register int j0;
1036  register int j1;
1037 
1038  register BigReal pj_x_0, pj_x_1;
1039  register BigReal pj_y_0, pj_y_1;
1040  register BigReal pj_z_0, pj_z_1;
1041  register BigReal t_0, t_1, r2_0, r2_1;
1042 
1043  pj_x_0 = p_1[jprev0].position.x;
1044  pj_x_1 = p_1[jprev1].position.x;
1045 
1046  pj_y_0 = p_1[jprev0].position.y;
1047  pj_y_1 = p_1[jprev1].position.y;
1048 
1049  pj_z_0 = p_1[jprev0].position.z;
1050  pj_z_1 = p_1[jprev1].position.z;
1051 
1052  g += 2;
1053  for ( ; g < gu - 2; g +=2 ) {
1054  // compute 1d distance, 2-way parallel
1055  j0 = jprev0;
1056  j1 = jprev1;
1057 
1058  t_0 = p_i_x - pj_x_0;
1059  t_1 = p_i_x - pj_x_1;
1060  r2_0 = t_0 * t_0;
1061  r2_1 = t_1 * t_1;
1062 
1063  t_0 = p_i_y - pj_y_0;
1064  t_1 = p_i_y - pj_y_1;
1065  r2_0 += t_0 * t_0;
1066  r2_1 += t_1 * t_1;
1067 
1068  t_0 = p_i_z - pj_z_0;
1069  t_1 = p_i_z - pj_z_1;
1070  r2_0 += t_0 * t_0;
1071  r2_1 += t_1 * t_1;
1072 
1073  #if NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR ( + 1 ) )
1074  sortEntry0 = sortValues + g;
1075  sortEntry1 = sortValues + g + 1;
1076  jprev0 = sortEntry0->index;
1077  jprev1 = sortEntry1->index;
1078  #else
1079  jprev0 = glist[g ];
1080  jprev1 = glist[g+1];
1081  #endif
1082 
1083  pj_x_0 = p_1[jprev0].position.x;
1084  pj_x_1 = p_1[jprev1].position.x;
1085  pj_y_0 = p_1[jprev0].position.y;
1086  pj_y_1 = p_1[jprev1].position.y;
1087  pj_z_0 = p_1[jprev0].position.z;
1088  pj_z_1 = p_1[jprev1].position.z;
1089 
1090  bool test0 = ( r2_0 < groupplcutoff2 );
1091  bool test1 = ( r2_1 < groupplcutoff2 );
1092 
1093  //removing ifs benefits on many architectures
1094  //as the extra stores will only warm the cache up
1095  goodglist [ hu ] = j0;
1096  goodglist [ hu + test0 ] = j1;
1097 
1098  hu += test0 + test1;
1099  }
1100  g-=2;
1101  }
1102 #endif
1103 #endif // NAMD_KNL
1104 
1105  const int gp = g;
1106 #pragma omp simd simdlen(16)
1107  for (g = gp; g < gu; g++) {
1108 
1109  #if NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR ( + 1 ) )
1110  register SortEntry* sortEntry = sortValues + g;
1111  register int j = sortEntry->index;
1112  #else
1113  int j = glist[g];
1114  #endif
1115 
1116  BigReal p_j_x = p_1[j].position.x;
1117  BigReal p_j_y = p_1[j].position.y;
1118  BigReal p_j_z = p_1[j].position.z;
1119 
1120  BigReal r2 = p_i_x - p_j_x;
1121  r2 *= r2;
1122  BigReal t2 = p_i_y - p_j_y;
1123  r2 += t2 * t2;
1124  t2 = p_i_z - p_j_z;
1125  r2 += t2 * t2;
1126 
1127 #pragma omp ordered simd monotonic(hu:1)
1128  if ( r2 <= groupplcutoff2 )
1129  goodglist[hu ++] = j;
1130  }
1131 
1132  for ( int h=0; h<hu; ++h ) {
1133  int j = goodglist[h];
1134  int nbgs = p_1[j].nonbondedGroupSize;
1135  pli[0] = j; // copy over the next three in any case
1136  pli[1] = j+1;
1137  pli[2] = j+2;
1138  if ( nbgs & 4 ) { // if nbgs > 3, assume nbgs <= 5
1139  pli[3] = j+3;
1140  pli[4] = j+4;
1141  }
1142  pli += nbgs;
1143  }
1144  }
1145  }
1146 
1147  pairlistindex = pli - pairlist;
1148  // make sure padded element on pairlist points to real data
1149  if ( pairlistindex ) {
1150  pairlist[pairlistindex] = pairlist[pairlistindex-1];
1151  } PAIR( else { // skip empty loops if no pairs were found
1152  i += p_i.nonbondedGroupSize;
1153  continue;
1154  } )
1155  } // if i is nonbonded group parent
1156  SELF
1157  (
1158  // self-comparisions require list to be incremented
1159  // pair-comparisions use entire list (pairlistoffset is 0)
1160  else pairlistoffset++;
1161  )
1162 
1163  const int atomfixed = ( fixedAtomsOn && pExt_i.atomFixed );
1164 
1165  register int *pli = pairlist2;
1166 
1167  plint *pairlistn = pairlists.newlist(j_upper + 5 + 5 ALCH( + 20 ));
1168  register plint *plin = pairlistn;
1169 
1170 
1171  if ( qmForcesOn ) {
1172 
1173  Real qmGroup_i = mol->get_qmAtomGroup(pExt_i.id);
1174 
1175  for (int k=pairlistoffset; k<pairlistindex; k++) {
1176  j = pairlist[k];
1177 
1178  Real qmGroup_j = mol->get_qmAtomGroup(pExt_1[j].id);
1179 
1180  // There can be no non-bonded interaction between an QM-QM pair of the same QM group
1181  if (qmGroup_i > 0) {
1182  if (qmGroup_i == qmGroup_j) {
1183  continue;
1184  } else {
1185  *(pli++) = j;
1186  }
1187  } else {
1188  *(pli++) = j;
1189  }
1190  }
1191 
1192  int npair2_int = pli - pairlist2;
1193  pli = pairlist2;
1194  for (int k=0; k<npair2_int; k++) {
1195 
1196  j = pairlist2[k];
1197 
1198  BigReal p_j_x = p_1[j].position.x;
1199  BigReal r2 = p_i_x - p_j_x;
1200  r2 *= r2;
1201  BigReal p_j_y = p_1[j].position.y;
1202  BigReal t2 = p_i_y - p_j_y;
1203  r2 += t2 * t2;
1204  BigReal p_j_z = p_1[j].position.z;
1205  t2 = p_i_z - p_j_z;
1206  r2 += t2 * t2;
1207 
1208  if ( ( ! (atomfixed && pExt_1[j].atomFixed) ) && (r2 <= plcutoff2) ) {
1209  int atom2 = pExt_1[j].id;
1210  if ( atom2 >= excl_min && atom2 <= excl_max ) *(pli++) = j;
1211  else *(plin++) = j;
1212  }
1213  }
1214  } else
1215 
1216 
1217 
1218  INT(
1219  if ( pairInteractionOn ) {
1220  const int ifep_type = p_i.partition;
1221  if (pairInteractionSelf) {
1222  if (ifep_type != 1) { PAIR(++i;) continue; }
1223  for (int k=pairlistoffset; k<pairlistindex; k++) {
1224  j = pairlist[k];
1225  const int jfep_type = p_1[j].partition;
1226  // for pair-self, both atoms must be in group 1.
1227  if (jfep_type == 1) {
1228  *(pli++) = j;
1229  }
1230  }
1231  } else {
1232  if (ifep_type != 1 && ifep_type != 2) { PAIR(++i;) continue; }
1233  for (int k=pairlistoffset; k<pairlistindex; k++) {
1234  j = pairlist[k];
1235  const int jfep_type = p_1[j].partition;
1236  // for pair, must have one from each group.
1237  if (ifep_type + jfep_type == 3) {
1238  *(pli++) = j;
1239  }
1240  }
1241  }
1242  int npair2_int = pli - pairlist2;
1243  pli = pairlist2;
1244  for (int k=0; k<npair2_int; k++) {
1245  j = pairlist2[k];
1246  BigReal p_j_x = p_1[j].position.x;
1247  BigReal r2 = p_i_x - p_j_x;
1248  r2 *= r2;
1249  BigReal p_j_y = p_1[j].position.y;
1250  BigReal t2 = p_i_y - p_j_y;
1251  r2 += t2 * t2;
1252  BigReal p_j_z = p_1[j].position.z;
1253  t2 = p_i_z - p_j_z;
1254  r2 += t2 * t2;
1255  if ( ( ! (atomfixed && pExt_1[j].atomFixed) ) && (r2 <= plcutoff2) ) {
1256  int atom2 = pExt_1[j].id;
1257  if ( atom2 >= excl_min && atom2 <= excl_max ) *(pli++) = j;
1258  else *(plin++) = j;
1259  }
1260  }
1261  } else
1262  )
1263  if ( atomfixed ) {
1264  for (int k=pairlistoffset; k<pairlistindex; k++) {
1265  j = pairlist[k];
1266  BigReal p_j_x = p_1[j].position.x;
1267  BigReal r2 = p_i_x - p_j_x;
1268  r2 *= r2;
1269  BigReal p_j_y = p_1[j].position.y;
1270  BigReal t2 = p_i_y - p_j_y;
1271  r2 += t2 * t2;
1272  BigReal p_j_z = p_1[j].position.z;
1273  t2 = p_i_z - p_j_z;
1274  r2 += t2 * t2;
1275  if ( (! pExt_1[j].atomFixed) && (r2 <= plcutoff2) ) {
1276  int atom2 = pExt_1[j].id;
1277  if ( atom2 >= excl_min && atom2 <= excl_max ) *(pli++) = j;
1278  else *(plin++) = j;
1279  }
1280  }
1281  } else {
1282  int k = pairlistoffset;
1283  int ku = pairlistindex;
1284  if ( k < ku ) {
1285 #ifndef NAMD_KNL
1286 #if defined(__SSE2__) && ! defined(NAMD_DISABLE_SSE)
1287  if ( ku - k > 6 ) {
1288  register int jprev0 = pairlist [k ];
1289  register int jprev1 = pairlist [k + 1];
1290 
1291  register int j0;
1292  register int j1;
1293 
1294  __m128d PJ_X_01 = _mm_set_pd(p_1[jprev1].position.x, p_1[jprev0].position.x);
1295  __m128d PJ_Y_01 = _mm_set_pd(p_1[jprev1].position.y, p_1[jprev0].position.y);
1296  __m128d PJ_Z_01 = _mm_set_pd(p_1[jprev1].position.z, p_1[jprev0].position.z);
1297 
1298  // these don't change here, so we could move them into outer scope
1299  const __m128d P_I_X = _mm_set1_pd(p_i_x);
1300  const __m128d P_I_Y = _mm_set1_pd(p_i_y);
1301  const __m128d P_I_Z = _mm_set1_pd(p_i_z);
1302 
1303  int atom2_0 = pExt_1[jprev0].id;
1304  int atom2_1 = pExt_1[jprev1].id;
1305 
1306  k += 2;
1307  for ( ; k < ku - 2; k +=2 ) {
1308  // compute 1d distance, 2-way parallel
1309  j0 = jprev0;
1310  j1 = jprev1;
1311 
1312  __m128d T_01 = _mm_sub_pd(P_I_X, PJ_X_01);
1313  __m128d R2_01 = _mm_mul_pd(T_01, T_01);
1314  T_01 = _mm_sub_pd(P_I_Y, PJ_Y_01);
1315  R2_01 = _mm_add_pd(R2_01, _mm_mul_pd(T_01, T_01));
1316  T_01 = _mm_sub_pd(P_I_Z, PJ_Z_01);
1317  R2_01 = _mm_add_pd(R2_01, _mm_mul_pd(T_01, T_01));
1318 
1319  jprev0 = pairlist[k];
1320  jprev1 = pairlist[k+1];
1321 
1322  PJ_X_01 = _mm_set_pd(p_1[jprev1].position.x, p_1[jprev0].position.x);
1323  PJ_Y_01 = _mm_set_pd(p_1[jprev1].position.y, p_1[jprev0].position.y);
1324  PJ_Z_01 = _mm_set_pd(p_1[jprev1].position.z, p_1[jprev0].position.z);
1325 
1326  __align(16) double r2_01[2];
1327  _mm_store_pd(r2_01, R2_01); // 16-byte-aligned store
1328 
1329  if (r2_01[0] <= plcutoff2) {
1330  if (UNLIKELY( atom2_0 >= excl_min && atom2_0 <= excl_max ))
1331  *(pli++) = j0;
1332  else
1333  *(plin++) = j0;
1334  }
1335  atom2_0 = pExt_1[jprev0].id;
1336 
1337  if (r2_01[1] <= plcutoff2) {
1338  if (UNLIKELY( atom2_1 >= excl_min && atom2_1 <= excl_max ))
1339  *(pli++) = j1;
1340  else
1341  *(plin++) = j1;
1342  }
1343  atom2_1 = pExt_1[jprev1].id;
1344  }
1345  k-=2;
1346  }
1347 #else
1348  if ( ku - k > 6 ) {
1349  register int jprev0 = pairlist [k];
1350  register int jprev1 = pairlist [k + 1];
1351 
1352  register int j0;
1353  register int j1;
1354 
1355  register BigReal pj_x_0, pj_x_1;
1356  register BigReal pj_y_0, pj_y_1;
1357  register BigReal pj_z_0, pj_z_1;
1358  register BigReal t_0, t_1, r2_0, r2_1;
1359 
1360  pj_x_0 = p_1[jprev0].position.x;
1361  pj_x_1 = p_1[jprev1].position.x;
1362 
1363  pj_y_0 = p_1[jprev0].position.y;
1364  pj_y_1 = p_1[jprev1].position.y;
1365 
1366  pj_z_0 = p_1[jprev0].position.z;
1367  pj_z_1 = p_1[jprev1].position.z;
1368 
1369  int atom2_0 = pExt_1[jprev0].id;
1370  int atom2_1 = pExt_1[jprev1].id;
1371 
1372  k += 2;
1373  for ( ; k < ku - 2; k +=2 ) {
1374  // compute 1d distance, 2-way parallel
1375  j0 = jprev0;
1376  j1 = jprev1;
1377 
1378  t_0 = p_i_x - pj_x_0;
1379  t_1 = p_i_x - pj_x_1;
1380  r2_0 = t_0 * t_0;
1381  r2_1 = t_1 * t_1;
1382 
1383  t_0 = p_i_y - pj_y_0;
1384  t_1 = p_i_y - pj_y_1;
1385  r2_0 += t_0 * t_0;
1386  r2_1 += t_1 * t_1;
1387 
1388  t_0 = p_i_z - pj_z_0;
1389  t_1 = p_i_z - pj_z_1;
1390  r2_0 += t_0 * t_0;
1391  r2_1 += t_1 * t_1;
1392 
1393  jprev0 = pairlist[k];
1394  jprev1 = pairlist[k+1];
1395 
1396  pj_x_0 = p_1[jprev0].position.x;
1397  pj_x_1 = p_1[jprev1].position.x;
1398  pj_y_0 = p_1[jprev0].position.y;
1399  pj_y_1 = p_1[jprev1].position.y;
1400  pj_z_0 = p_1[jprev0].position.z;
1401  pj_z_1 = p_1[jprev1].position.z;
1402 
1403  if (r2_0 <= plcutoff2) {
1404  if ( atom2_0 >= excl_min && atom2_0 <= excl_max )
1405  *(pli++) = j0;
1406  else
1407  *(plin++) = j0;
1408  }
1409  atom2_0 = pExt_1[jprev0].id;
1410 
1411  if (r2_1 <= plcutoff2) {
1412  if ( atom2_1 >= excl_min && atom2_1 <= excl_max )
1413  *(pli++) = j1;
1414  else
1415  *(plin++) = j1;
1416  }
1417  atom2_1 = pExt_1[jprev1].id;
1418  }
1419  k-=2;
1420  }
1421 #endif
1422 #endif // NAMD_KNL
1423 
1424  const int kp = k;
1425 #pragma omp simd simdlen(16)
1426  for (k = kp; k < ku; k++) {
1427  int j = pairlist[k];
1428  int atom2 = pExt_1[j].id;
1429 
1430  BigReal p_j_x = p_1[j].position.x;
1431  BigReal p_j_y = p_1[j].position.y;
1432  BigReal p_j_z = p_1[j].position.z;
1433 
1434  BigReal r2 = p_i_x - p_j_x;
1435  r2 *= r2;
1436  BigReal t2 = p_i_y - p_j_y;
1437  r2 += t2 * t2;
1438  t2 = p_i_z - p_j_z;
1439  r2 += t2 * t2;
1440 
1441 #pragma omp ordered simd monotonic(plin:1, pli:1)
1442  if (r2 <= plcutoff2) {
1443  if ( atom2 >= excl_min && atom2 <= excl_max )
1444  *(pli++) = j;
1445  else
1446  *(plin++) = j;
1447  }
1448  }
1449  }
1450  }
1451 
1452  PAIR(
1453  if ( plin == pairlistn && pli == pairlist2 ) {
1454  ++i;
1455  continue;
1456  }
1457  pairlists.setIndexValue(i);
1458  )
1459 
1460  plint *plix = pairlistx;
1461  plint *plim = pairlistm;
1462  ALCH(
1463  plint *plinA1 = pairlistnA1;
1464  plint *plixA1 = pairlistxA1;
1465  plint *plimA1 = pairlistmA1;
1466  plint *plinA2 = pairlistnA2;
1467  plint *plixA2 = pairlistxA2;
1468  plint *plimA2 = pairlistmA2;
1469  plint *plinA3 = pairlistnA3;
1470  plint *plixA3 = pairlistxA3;
1471  plint *plimA3 = pairlistmA3;
1472  plint *plinA4 = pairlistnA4;
1473  plint *plixA4 = pairlistxA4;
1474  plint *plimA4 = pairlistmA4;
1475  )
1476 
1477  int k;
1478 
1479 #if 0 ALCH(+1)
1480  int unsortedNpairn = plin - pairlistn;
1481  plin = pairlistn;
1482  for ( k=0; k<unsortedNpairn; ++k ) {
1483  int j = pairlistn[k];
1484  switch(pswitchTable[p_i_partition + 5*(p_1[j].partition)]) {
1485  case 0: *(plin++) = j; break;
1486  case 1: *(plinA1++) = j; break;
1487  case 2: *(plinA2++) = j; break;
1488  case 3: *(plinA3++) = j; break;
1489  case 4: *(plinA4++) = j; break;
1490  }
1491  }
1492 #endif
1493 
1494  int npair2 = pli - pairlist2;
1495  // if ( npair2 ) pairlist2[npair2] = pairlist2[npair2-1];
1496  // removed code for implicit exclusions within hydrogen groups -JCP
1497  for (k=0; k < npair2; ++k ) {
1498  int j = pairlist2[k];
1499  int atom2 = pExt_1[j].id;
1500  int excl_flag = excl_flags[atom2];
1501  ALCH(int pswitch = pswitchTable[p_i_partition + 5*(p_1[j].partition)];)
1502  switch ( excl_flag ALCH( + 3 * pswitch)) {
1503  case 0: *(plin++) = j; break;
1504  case 1: *(plix++) = j; break;
1505  case 2: *(plim++) = j; break;
1506  ALCH(
1507  case 3: *(plinA1++) = j; break;
1508  case 6: *(plinA2++) = j; break;
1509  case 9: *(plinA3++) = j; break;
1510  case 12: *(plinA4++) = j; break;
1511  case 5: *(plimA1++) = j; break;
1512  case 8: *(plimA2++) = j; break;
1513  case 11: *(plimA3++) = j; break;
1514  case 14: *(plimA4++) = j; break;
1515  case 4: *(plixA1++) = j; break;
1516  case 7: *(plixA2++) = j; break;
1517  case 10: *(plixA3++) = j; break;
1518  case 13: *(plixA4++) = j; break;
1519  )
1520  }
1521  }
1522 
1523  npairn = plin - pairlistn;
1524  pairlistn_save = pairlistn;
1525  pairlistn_save[npairn] = npairn ? pairlistn_save[npairn-1] : -1;
1526  pairlists.newsize(npairn + 1);
1527 
1528  npairx = plix - pairlistx;
1529  pairlistx_save = pairlists.newlist();
1530  for ( k=0; k<npairx; ++k ) {
1531  pairlistx_save[k] = pairlistx[k];
1532  }
1533  pairlistx_save[k] = k ? pairlistx_save[k-1] : -1;
1534  pairlists.newsize(npairx + 1);
1535 
1536  npairm = plim - pairlistm;
1537  pairlistm_save = pairlists.newlist();
1538  for ( k=0; k<npairm; ++k ) {
1539  pairlistm_save[k] = pairlistm[k];
1540  }
1541  pairlistm_save[k] = k ? pairlistm_save[k-1] : -1;
1542  pairlists.newsize(npairm + 1);
1543 
1544 
1545 #if 0 ALCH(+1)
1546 #define PAIRLISTFROMARRAY(NPAIRS,PL1,PL2,PLSAVE) \
1547 { \
1548  (NPAIRS) = (PL2) - (PL1); \
1549  (PLSAVE) = pairlists.newlist(); \
1550  for ( k=0; k<(NPAIRS); ++k ) { \
1551  (PLSAVE)[k] = (PL1)[k]; \
1552  } \
1553  (PLSAVE)[k] = k ? (PLSAVE)[k-1] : -1; \
1554  pairlists.newsize((NPAIRS) + 1); \
1555 }
1556 
1557  PAIRLISTFROMARRAY(npairnA1,pairlistnA1,plinA1,pairlistnA1_save);
1558  PAIRLISTFROMARRAY(npairxA1,pairlistxA1,plixA1,pairlistxA1_save);
1559  PAIRLISTFROMARRAY(npairmA1,pairlistmA1,plimA1,pairlistmA1_save);
1560  PAIRLISTFROMARRAY(npairnA2,pairlistnA2,plinA2,pairlistnA2_save);
1561  PAIRLISTFROMARRAY(npairxA2,pairlistxA2,plixA2,pairlistxA2_save);
1562  PAIRLISTFROMARRAY(npairmA2,pairlistmA2,plimA2,pairlistmA2_save);
1563  PAIRLISTFROMARRAY(npairnA3,pairlistnA3,plinA3,pairlistnA3_save);
1564  PAIRLISTFROMARRAY(npairxA3,pairlistxA3,plixA3,pairlistxA3_save);
1565  PAIRLISTFROMARRAY(npairmA3,pairlistmA3,plimA3,pairlistmA3_save);
1566  PAIRLISTFROMARRAY(npairnA4,pairlistnA4,plinA4,pairlistnA4_save);
1567  PAIRLISTFROMARRAY(npairxA4,pairlistxA4,plixA4,pairlistxA4_save);
1568  PAIRLISTFROMARRAY(npairmA4,pairlistmA4,plimA4,pairlistmA4_save);
1569 #undef PAIRLISTFROMARRAY
1570 
1571 #endif
1572 
1573  if ( ! savePairlists ) pairlists.reset(); // limit space usage
1574  PAIR( pairlists.addIndex(); )
1575 
1576  // PAIR( iout << i << " " << i_upper << " save\n" << endi;)
1577  } else { // if ( savePairlists || ! usePairlists )
1578  // PAIR( iout << i << " " << i_upper << " use\n" << endi;)
1579 
1580  pairlists.nextlist(&pairlistn_save,&npairn); --npairn;
1581  pairlists.nextlist(&pairlistx_save,&npairx); --npairx;
1582  pairlists.nextlist(&pairlistm_save,&npairm); --npairm;
1583  ALCH(
1584  pairlists.nextlist(&pairlistnA1_save,&npairnA1); --npairnA1;
1585  pairlists.nextlist(&pairlistxA1_save,&npairxA1); --npairxA1;
1586  pairlists.nextlist(&pairlistmA1_save,&npairmA1); --npairmA1;
1587  pairlists.nextlist(&pairlistnA2_save,&npairnA2); --npairnA2;
1588  pairlists.nextlist(&pairlistxA2_save,&npairxA2); --npairxA2;
1589  pairlists.nextlist(&pairlistmA2_save,&npairmA2); --npairmA2;
1590  pairlists.nextlist(&pairlistnA3_save,&npairnA3); --npairnA3;
1591  pairlists.nextlist(&pairlistxA3_save,&npairxA3); --npairxA3;
1592  pairlists.nextlist(&pairlistmA3_save,&npairmA3); --npairmA3;
1593  pairlists.nextlist(&pairlistnA4_save,&npairnA4); --npairnA4;
1594  pairlists.nextlist(&pairlistxA4_save,&npairxA4); --npairxA4;
1595  pairlists.nextlist(&pairlistmA4_save,&npairmA4); --npairmA4;
1596  )
1597 
1598  } // if ( savePairlists || ! usePairlists )
1599 
1600  LES( BigReal *lambda_table_i =
1601  lambda_table + (lesFactor+1) * p_i.partition; )
1602 
1603  INT(
1604  const BigReal force_sign = (
1605  ( pairInteractionOn && ! pairInteractionSelf ) ?
1606  ( ( p_i.partition == 1 ) ? 1. : -1. ) : 0.
1607  );
1608 
1609  )
1610 
1611  const BigReal kq_i = COULOMB * p_i.charge * scaling * dielectric_1;
1612  const LJTable::TableEntry * const lj_row =
1613  ljTable->table_row(p_i.vdwType);
1614 
1615  SHORT( FAST( BigReal f_i_x = 0.; ) )
1616  SHORT( FAST( BigReal f_i_y = 0.; ) )
1617  SHORT( FAST( BigReal f_i_z = 0.; ) )
1618  FULL( BigReal fullf_i_x = 0.; )
1619  FULL( BigReal fullf_i_y = 0.; )
1620  FULL( BigReal fullf_i_z = 0.; )
1621 
1622  int npairi;
1623  int k;
1624 
1625 #if KNL(1)+0
1626  const float kq_i_f = kq_i;
1627 
1628  npairi = pairlist_from_pairlist_knl(ComputeNonbondedUtil::cutoff2_f,
1629  p_i_x_f, p_i_y_f, p_i_z_f, pFlt_1, pairlistn_save, npairn, pairlisti,
1630  r2list_f, xlist, ylist, zlist);
1631 #else
1632  npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
1633  p_i_x, p_i_y, p_i_z, p_1, pairlistn_save, npairn, pairlisti,
1634  r2_delta, r2list);
1635 #endif
1636 
1637  if ( npairi ) {
1638 
1639 // BEGIN NBTHOLE OF DRUDE MODEL
1640 #if (FAST(1+)0)
1641  if (drudeNbthole && p_i.hydrogenGroupSize > 1) {
1642 
1643  Parameters *parameters = params->parameters;
1644  const NbtholePairValue * const nbthole_array = parameters->nbthole_array;
1645  const int NumNbtholePairParams = parameters->NumNbtholePairParams;
1646  BigReal drudeNbtholeCut = simParams -> drudeNbtholeCut;
1647  NOKNL( BigReal drudeNbtholeCut2 = (drudeNbtholeCut * drudeNbtholeCut) + r2_delta; )
1648  KNL( float drudeNbtholeCut2 = (drudeNbtholeCut * drudeNbtholeCut); )
1649  BigReal CC = COULOMB * scaling * dielectric_1;
1650  int kk;
1651 
1652  for (k = 0; k < npairi; k++) {
1653  NOKNL( if (r2list[k] > drudeNbtholeCut2) { continue; } )
1654  KNL( if (r2list_f[k] > drudeNbtholeCut2) { continue; } )
1655 
1656  const int j = pairlisti[k];
1657  const CompAtom& p_j = p_1[j];
1658 
1659  if ( p_j.hydrogenGroupSize < 2 ) { continue; }
1660 
1661  for (kk = 0;kk < NumNbtholePairParams; kk++) {
1662 
1663  if (((nbthole_array[kk].ind1 == p_i.vdwType) &&
1664  (nbthole_array[kk].ind2 == p_j.vdwType)) ||
1665  ((nbthole_array[kk].ind2 == p_i.vdwType) &&
1666  (nbthole_array[kk].ind1 == p_j.vdwType))) {
1667  break;
1668  }
1669  }
1670  if ( kk < NumNbtholePairParams ) {
1671 
1672  BigReal aprod = mol->GetAtomAlpha(pExt_0[i].id) * mol->GetAtomAlpha(pExt_1[j].id);
1673  const BigReal aa = nbthole_array[kk].tholeij * powf(aprod, -1.f/6);
1674  const BigReal qqaa = CC * p_0[i].charge * p_1[j].charge;
1675  const BigReal qqad = CC * p_0[i].charge * p_1[j+1].charge;
1676  const BigReal qqda = CC * p_0[i+1].charge * p_1[j].charge;
1677  const BigReal qqdd = CC * p_0[i+1].charge * p_1[j+1].charge;
1678 
1679  Vector raa = (p_0[i].position + params->offset) - p_1[j].position;
1680  Vector rad = (p_0[i].position + params->offset) - p_1[j+1].position;
1681  Vector rda = (p_0[i+1].position + params->offset) - p_1[j].position;
1682  Vector rdd = (p_0[i+1].position + params->offset) - p_1[j+1].position;
1683 
1684  BigReal raa_invlen = raa.rlength(); // reciprocal of length
1685  BigReal rad_invlen = rad.rlength();
1686  BigReal rda_invlen = rda.rlength();
1687  BigReal rdd_invlen = rdd.rlength();
1688 
1689  BigReal auaa = aa / raa_invlen;
1690  BigReal auad = aa / rad_invlen;
1691  BigReal auda = aa / rda_invlen;
1692  BigReal audd = aa / rdd_invlen;
1693 
1694  BigReal expauaa = exp(-auaa);
1695  BigReal expauad = exp(-auad);
1696  BigReal expauda = exp(-auda);
1697  BigReal expaudd = exp(-audd);
1698 
1699  BigReal polyauaa = 1 + 0.5*auaa;
1700  BigReal polyauad = 1 + 0.5*auad;
1701  BigReal polyauda = 1 + 0.5*auda;
1702  BigReal polyaudd = 1 + 0.5*audd;
1703 
1704  BigReal ethole = 0;
1705  ethole += qqaa * raa_invlen * (- polyauaa * expauaa);
1706  ethole += qqad * rad_invlen * (- polyauad * expauad);
1707  ethole += qqda * rda_invlen * (- polyauda * expauda);
1708  ethole += qqdd * rdd_invlen * (- polyaudd * expaudd);
1709 
1710  polyauaa = 1 + auaa*polyauaa;
1711  polyauad = 1 + auad*polyauad;
1712  polyauda = 1 + auda*polyauda;
1713  polyaudd = 1 + audd*polyaudd;
1714 
1715  BigReal raa_invlen3 = raa_invlen * raa_invlen * raa_invlen;
1716  BigReal rad_invlen3 = rad_invlen * rad_invlen * rad_invlen;
1717  BigReal rda_invlen3 = rda_invlen * rda_invlen * rda_invlen;
1718  BigReal rdd_invlen3 = rdd_invlen * rdd_invlen * rdd_invlen;
1719 
1720  BigReal dfaa = qqaa * raa_invlen3 * (polyauaa*expauaa);
1721  BigReal dfad = qqad * rad_invlen3 * (polyauad*expauad);
1722  BigReal dfda = qqda * rda_invlen3 * (polyauda*expauda);
1723  BigReal dfdd = qqdd * rdd_invlen3 * (polyaudd*expaudd);
1724 
1725  Vector faa = -dfaa * raa;
1726  Vector fad = -dfad * rad;
1727  Vector fda = -dfda * rda;
1728  Vector fdd = -dfdd * rdd;
1729 
1730  SHORT(f_net) NOSHORT(fullf_net) += (faa + fad) + (fda + fdd);
1731  params->ff[0][i] += faa + fad;
1732  params->ff[0][i+1] += fda + fdd;
1733  params->ff[1][j] -= faa + fda;
1734  params->ff[1][j+1] -= fad + fdd;
1735 
1736  reduction[electEnergyIndex] += ethole;
1737  }
1738  }
1739  }
1740 #endif
1741 // END NBTHOLE OF DRUDE MODEL
1742 
1743  // BEGIN LA
1744 #if (FAST(1+)0)
1745  if (doLoweAndersen && p_i.hydrogenGroupSize) {
1746  BigReal loweAndersenCutoff = simParams->loweAndersenCutoff;
1747  NOKNL( BigReal loweAndersenCutoff2 = (loweAndersenCutoff * loweAndersenCutoff) + r2_delta; )
1748  KNL( float loweAndersenCutoff2 = (loweAndersenCutoff * loweAndersenCutoff); )
1749  BigReal loweAndersenProb = simParams->loweAndersenRate * (simParams->dt * simParams->nonbondedFrequency) * 0.001; // loweAndersenRate is in 1/ps
1750  const bool loweAndersenUseCOMvelocity = (simParams->rigidBonds != RIGID_NONE);
1751  const BigReal kbT = BOLTZMANN * (simParams->loweAndersenTemp);
1752  const BigReal dt_inv = TIMEFACTOR / (simParams->dt * simParams->nonbondedFrequency);
1753  //const BigReal dt_inv = 1.0 / simParams->dt;
1754  //BigReal kbT = 8.3145e-7 * (simParams->loweAndersenTemp); // in A^2/fs^2 * K
1755 
1756  const CompAtom* v_0 = params->v[0];
1757  const CompAtom* v_1 = params->v[1];
1758  const CompAtom& v_i = v_0[i];
1759  Mass mass_i = v_i.charge;
1760  Velocity vel_i = v_i.position;
1761  Position pos_i = p_i.position;
1762  int atom_i = pExt_0[i].id;
1763 
1764  if (loweAndersenUseCOMvelocity) {
1765  Mass mass = 0;
1766  Velocity vel = 0;
1767  Position pos = 0;
1768  for (int l = 0; l < p_i.hydrogenGroupSize; l++) {
1769  vel += v_0[i+l].charge * v_0[i+l].position;
1770  pos += v_0[i+l].charge * p_0[i+l].position;
1771  mass += v_0[i+l].charge;
1772  }
1773  vel_i = vel / mass;
1774  pos_i = pos / mass;
1775  mass_i = mass;
1776  }
1777 
1778  // Note v[0].charge is actually mass
1779  //Random rand(CkMyPe()); // ?? OK ?? NO!!!!
1780  Random *rand = params->random;
1781  for (k = 0; k < npairi; k++) {
1782  NOKNL( if (r2list[k] > loweAndersenCutoff2) { continue; } )
1783  KNL( if (r2list_f[k] > loweAndersenCutoff2) { continue; } )
1784 
1785  const int j = pairlisti[k];
1786  const CompAtom& v_j = v_1[j];
1787  const CompAtom& p_j = p_1[j];
1788 
1789  if (!p_j.hydrogenGroupSize) { continue; }
1790  if (rand->uniform() > loweAndersenProb) { continue; }
1791 
1792  Mass mass_j = v_j.charge;
1793  Velocity vel_j = v_j.position;
1794  Position pos_j = p_j.position;
1795  int atom_j = pExt_1[j].id;
1796 
1797  if (loweAndersenUseCOMvelocity) {
1798  Mass mass = 0;
1799  Velocity vel = 0;
1800  Position pos = 0;
1801  for (int l = 0; l < p_j.hydrogenGroupSize; l++) {
1802  vel += v_1[j+l].charge * v_1[j+l].position;
1803  pos += v_1[j+l].charge * p_1[j+l].position;
1804  mass += v_1[j+l].charge;
1805  }
1806  vel_j = vel / mass;
1807  pos_j = pos / mass;
1808  mass_j = mass;
1809  }
1810 
1811  //Velocity deltaV = v_i.position - v_j.position;
1812  Velocity deltaV = vel_i - vel_j;
1813  Mass mu_ij = (mass_i * mass_j)/(mass_i + mass_j);
1814  //Vector sep = (p_i.position + params->offset) - p_j.position;
1815  Vector sep = (pos_i + params->offset) - pos_j;
1816  Vector sigma_ij = sep.unit();
1817  BigReal lambda = rand->gaussian() * sqrt(kbT / mu_ij);
1818  Force force = mu_ij * dt_inv * (lambda - (deltaV * sigma_ij)) * sigma_ij;
1819 
1820  //DebugM(1, "atom1 atom2 = " << atom1 << " " << atom2 << " lambda = " << lambda << " force = " << force << " deltaP = " << sep.length() << " sqrt(r2) = " << sqrt(r2list[k]) << "\n");
1821  //DebugM(1, "atom1 atom2 = " << atom_i << " " << atom_j << " mass1 mass2 = " << mass_i << " " << mass_j << " hgrp1 hgrp2 " << p_i.hydrogenGroupSize << " " << p_j.hydrogenGroupSize << "\n");
1822 
1823  if (loweAndersenUseCOMvelocity) {
1824  BigReal inv_mass_i = 1.0 / mass_i;
1825  BigReal inv_mass_j = 1.0 / mass_j;
1826  for (int l = 0; l < p_i.hydrogenGroupSize; l++) {
1827  params->ff[0][i+l] += (v_0[i+l].charge * inv_mass_i) * force;
1828  }
1829  for (int l = 0; l < p_j.hydrogenGroupSize; l++) {
1830  params->ff[1][j+l] -= (v_1[j+l].charge * inv_mass_j) * force;
1831  }
1832  } else {
1833  params->ff[0][i] += force;
1834  params->ff[1][j] -= force;
1835  }
1836  SHORT(f_net) NOSHORT(fullf_net) += force;
1837  }
1838  }
1839 #endif
1840  // END LA
1841 
1842 #define NORMAL(X) X
1843 #define EXCLUDED(X)
1844 #define MODIFIED(X)
1845 #define PRAGMA_SIMD
1846 #if KNL(1+)0
1847  switch ( vdw_switch_mode ) {
1848 
1849 #define VDW_SWITCH_MODE VDW_SWITCH_MODE_ENERGY
1850  case VDW_SWITCH_MODE:
1851 #include "ComputeNonbondedBase2KNL.h"
1852  break;
1853 #undef VDW_SWITCH_MODE
1854 
1855 #define VDW_SWITCH_MODE VDW_SWITCH_MODE_MARTINI
1856  case VDW_SWITCH_MODE:
1857 #include "ComputeNonbondedBase2KNL.h"
1858  break;
1859 #undef VDW_SWITCH_MODE
1860 
1861 #define VDW_SWITCH_MODE VDW_SWITCH_MODE_FORCE
1862  case VDW_SWITCH_MODE:
1863 #include "ComputeNonbondedBase2KNL.h"
1864  break;
1865 #undef VDW_SWITCH_MODE
1866 
1867  }
1868 #else
1869 #include "ComputeNonbondedBase2.h"
1870 #endif
1871 #undef PRAGMA_SIMD
1872 #undef NORMAL
1873 #undef EXCLUDED
1874 #undef MODIFIED
1875 
1876  } // if ( npairi )
1877 
1878  npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
1879  p_i_x, p_i_y, p_i_z, p_1, pairlistm_save, npairm, pairlisti,
1880  r2_delta, r2list);
1881  exclChecksum += npairi;
1882 
1883 #define NORMAL(X)
1884 #define EXCLUDED(X)
1885 #define MODIFIED(X) X
1886 #include "ComputeNonbondedBase2.h"
1887 #undef NORMAL
1888 #undef EXCLUDED
1889 #undef MODIFIED
1890 
1891 #ifdef FULLELECT
1892  npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
1893  p_i_x, p_i_y, p_i_z, p_1, pairlistx_save, npairx, pairlisti,
1894  r2_delta, r2list);
1895  exclChecksum += npairi;
1896 
1897 #undef FAST
1898 #define FAST(X)
1899 #define NORMAL(X)
1900 #define EXCLUDED(X) X
1901 #define MODIFIED(X)
1902 #include "ComputeNonbondedBase2.h"
1903 #undef FAST
1904 #ifdef SLOWONLY
1905  #define FAST(X)
1906 #else
1907  #define FAST(X) X
1908 #endif
1909 #undef NORMAL
1910 #undef EXCLUDED
1911 #undef MODIFIED
1912 #else
1913  exclChecksum += npairx;
1914 #endif
1915 
1916 #if 0 ALCH(+1) // nonbondedbase2 for alchemical-separeted pairlists
1917 
1918  #undef ALCHPAIR
1919  #define ALCHPAIR(X) X
1920  #undef NOT_ALCHPAIR
1921  #define NOT_ALCHPAIR(X)
1922  BigReal myLambda; FEP(BigReal myLambda2;)
1923  BigReal myElecLambda; FEP(BigReal myElecLambda2;)
1924  BigReal myVdwLambda; FEP(BigReal myVdwLambda2;)
1925  BigReal myRepLambda; FEP(BigReal myRepLambda2;)
1926  BigReal myVdwShift; FEP(BigReal myVdwShift2;)
1927  BigReal alch_vdw_energy; BigReal alch_vdw_force;
1928  FEP(BigReal alch_vdw_energy_2; ) TI(BigReal alch_vdw_dUdl;)
1929  BigReal shiftedElec; BigReal shiftedElecForce;
1930 
1931  /********************************************************************/
1932  /*******NONBONDEDBASE2 FOR NORMAL INTERACTIONS SCALED BY LAMBDA******/
1933  /********************************************************************/
1934  #define NORMAL(X) X
1935  #define EXCLUDED(X)
1936  #define MODIFIED(X)
1937 
1938  #define ALCH1(X) X
1939  #define ALCH2(X)
1940  #define ALCH3(X)
1941  #define ALCH4(X)
1942  npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
1943  p_i_x, p_i_y, p_i_z, p_1, pairlistnA1_save, npairnA1, pairlisti,
1944  r2_delta, r2list);
1945  #include "ComputeNonbondedBase2.h" // normal, direction 'up'
1946  #undef ALCH1
1947  #undef ALCH2
1948  #undef ALCH3
1949  #undef ALCH4
1950 
1951  #define ALCH1(X)
1952  #define ALCH2(X) X
1953  #define ALCH3(X)
1954  #define ALCH4(X)
1955  npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
1956  p_i_x, p_i_y, p_i_z, p_1, pairlistnA2_save, npairnA2, pairlisti,
1957  r2_delta, r2list);
1958  #include "ComputeNonbondedBase2.h" // normal, direction 'down'
1959  #undef ALCH1
1960  #undef ALCH2
1961  #undef ALCH3
1962  #undef ALCH4
1963 
1964  #define ALCH3(X) X
1965  #define ALCH1(X)
1966  #define ALCH2(X)
1967  #define ALCH4(X)
1968  npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
1969  p_i_x, p_i_y, p_i_z, p_1, pairlistnA3_save, npairnA3, pairlisti,
1970  r2_delta, r2list);
1971  #include "ComputeNonbondedBase2.h" // normal, direction 'up'
1972  #undef ALCH1
1973  #undef ALCH2
1974  #undef ALCH3
1975  #undef ALCH4
1976 
1977  #define ALCH4(X) X
1978  #define ALCH1(X)
1979  #define ALCH2(X)
1980  #define ALCH3(X)
1981  npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
1982  p_i_x, p_i_y, p_i_z, p_1, pairlistnA4_save, npairnA4, pairlisti,
1983  r2_delta, r2list);
1984  #include "ComputeNonbondedBase2.h" // normal, direction 'down'
1985  #undef ALCH1
1986  #undef ALCH2
1987  #undef ALCH3
1988  #undef ALCH4
1989 
1990  #undef NORMAL
1991  #undef EXCLUDED
1992  #undef MODIFIED
1993  /********************************************************************/
1994  /******NONBONDEDBASE2 FOR MODIFIED INTERACTIONS SCALED BY LAMBDA*****/
1995  /********************************************************************/
1996  #define NORMAL(X)
1997  #define EXCLUDED(X)
1998  #define MODIFIED(X) X
1999 
2000  #define ALCH1(X) X
2001  #define ALCH2(X)
2002  #define ALCH3(X)
2003  #define ALCH4(X)
2004  npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
2005  p_i_x, p_i_y, p_i_z, p_1, pairlistmA1_save, npairmA1, pairlisti,
2006  r2_delta, r2list);
2007  exclChecksum += npairi;
2008  #include "ComputeNonbondedBase2.h" // modified, direction 'up'
2009  #undef ALCH1
2010  #undef ALCH2
2011  #undef ALCH3
2012  #undef ALCH4
2013 
2014  #define ALCH1(X)
2015  #define ALCH2(X) X
2016  #define ALCH3(X)
2017  #define ALCH4(X)
2018  npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
2019  p_i_x, p_i_y, p_i_z, p_1, pairlistmA2_save, npairmA2, pairlisti,
2020  r2_delta, r2list);
2021  exclChecksum += npairi;
2022  #include "ComputeNonbondedBase2.h" // modified, direction 'down'
2023  #undef ALCH1
2024  #undef ALCH2
2025  #undef ALCH3
2026  #undef ALCH4
2027 
2028  #define ALCH1(X)
2029  #define ALCH2(X)
2030  #define ALCH3(X) X
2031  #define ALCH4(X)
2032  npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
2033  p_i_x, p_i_y, p_i_z, p_1, pairlistmA3_save, npairmA3, pairlisti,
2034  r2_delta, r2list);
2035  exclChecksum += npairi;
2036  #include "ComputeNonbondedBase2.h" // modified, direction 'up'
2037  #undef ALCH1
2038  #undef ALCH2
2039  #undef ALCH3
2040  #undef ALCH4
2041 
2042  #define ALCH1(X)
2043  #define ALCH2(X)
2044  #define ALCH3(X)
2045  #define ALCH4(X) X
2046  npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
2047  p_i_x, p_i_y, p_i_z, p_1, pairlistmA4_save, npairmA4, pairlisti,
2048  r2_delta, r2list);
2049  exclChecksum += npairi;
2050  #include "ComputeNonbondedBase2.h" // modified, direction 'down'
2051  #undef ALCH1
2052  #undef ALCH2
2053  #undef ALCH3
2054  #undef ALCH4
2055 
2056  #undef NORMAL
2057  #undef EXCLUDED
2058  #undef MODIFIED
2059 
2060  /********************************************************************/
2061  /******NONBONDEDBASE2 FOR EXCLUDED INTERACTIONS SCALED BY LAMBDA*****/
2062  /********************************************************************/
2063  #ifdef FULLELECT
2064  #undef FAST
2065  #define FAST(X)
2066  #define NORMAL(X)
2067  #define EXCLUDED(X) X
2068  #define MODIFIED(X)
2069 
2070  #define ALCH1(X) X
2071  #define ALCH2(X)
2072  #define ALCH3(X)
2073  #define ALCH4(X)
2074  npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
2075  p_i_x, p_i_y, p_i_z, p_1, pairlistxA1_save, npairxA1, pairlisti,
2076  r2_delta, r2list);
2077  exclChecksum += npairi;
2078  #include "ComputeNonbondedBase2.h" //excluded, direction 'up'
2079  #undef ALCH1
2080  #undef ALCH2
2081  #undef ALCH3
2082  #undef ALCH4
2083 
2084  #define ALCH1(X)
2085  #define ALCH2(X) X
2086  #define ALCH3(X)
2087  #define ALCH4(X)
2088  npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
2089  p_i_x, p_i_y, p_i_z, p_1, pairlistxA2_save, npairxA2, pairlisti,
2090  r2_delta, r2list);
2091  exclChecksum += npairi;
2092  #include "ComputeNonbondedBase2.h" //excluded, direction 'down'
2093  #undef ALCH1
2094  #undef ALCH2
2095  #undef ALCH3
2096  #undef ALCH4
2097 
2098  #define ALCH1(X)
2099  #define ALCH2(X)
2100  #define ALCH3(X) X
2101  #define ALCH4(X)
2102  npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
2103  p_i_x, p_i_y, p_i_z, p_1, pairlistxA3_save, npairxA3, pairlisti,
2104  r2_delta, r2list);
2105  exclChecksum += npairi;
2106  #include "ComputeNonbondedBase2.h"
2107  #undef ALCH1
2108  #undef ALCH2
2109  #undef ALCH3
2110  #undef ALCH4
2111 
2112  #define ALCH1(X)
2113  #define ALCH2(X)
2114  #define ALCH3(X)
2115  #define ALCH4(X) X
2116  npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
2117  p_i_x, p_i_y, p_i_z, p_1, pairlistxA4_save, npairxA4, pairlisti,
2118  r2_delta, r2list);
2119  exclChecksum += npairi;
2120  #include "ComputeNonbondedBase2.h"
2121  #undef ALCH1
2122  #undef ALCH2
2123  #undef ALCH3
2124  #undef ALCH4
2125 
2126  #undef FAST
2127  #ifdef SLOWONLY
2128  #define FAST(X)
2129  #else
2130  #define FAST(X) X
2131  #endif
2132  #undef NORMAL
2133  #undef EXCLUDED
2134  #undef MODIFIED
2135  #else
2136  exclChecksum += npairxA1 + npairxA2 + npairxA3 + npairxA4;
2137  #endif
2138 
2139  #undef ALCHPAIR
2140  #define ALCHPAIR(X)
2141  #undef NOT_ALCHPAIR
2142  #define NOT_ALCHPAIR(X) X
2143 
2144 #endif // end nonbondedbase2.h includes for alchemical pairlists
2145 
2146 #ifdef NETWORK_PROGRESS
2147  CkNetworkProgress();
2148 #endif
2149 
2150 #ifdef ARCH_POWERPC
2151  //data cache block touch the position structure
2152  __dcbt ((void *) &(p_0[i+1]));
2153  //__prefetch_by_load ((void *)&(groupCount));
2154 #endif
2155 
2156  SHORT( FAST( f_net.x += f_i_x; ) )
2157  SHORT( FAST( f_net.y += f_i_y; ) )
2158  SHORT( FAST( f_net.z += f_i_z; ) )
2159  FULL( fullf_net.x += fullf_i_x; )
2160  FULL( fullf_net.y += fullf_i_y; )
2161  FULL( fullf_net.z += fullf_i_z; )
2162  SHORT( FAST( f_0[i].x += f_i_x; ) )
2163  SHORT( FAST( f_0[i].y += f_i_y; ) )
2164  SHORT( FAST( f_0[i].z += f_i_z; ) )
2165  FULL( fullf_0[i].x += fullf_i_x; )
2166  FULL( fullf_0[i].y += fullf_i_y; )
2167  FULL( fullf_0[i].z += fullf_i_z; )
2168 PAIR(
2169  if ( savePairlists || ! usePairlists ) {
2170  i++;
2171  } else {
2172  i = pairlists.getIndexValue();
2173  }
2174 )
2175 
2176  // PAIR( iout << i << " " << i_upper << " end\n" << endi;)
2177  } // for i
2178 
2179  // PAIR(iout << "++++++++\n" << endi;)
2180  PAIR( if ( savePairlists ) { pairlists.setIndexValue(i); } )
2181 
2182 #ifdef f_1
2183 #undef f_1
2184 #endif
2185 #if ( SHORT( FAST( 1+ ) ) 0 )
2186 #if ( PAIR( 1+ ) 0 )
2187  {
2188  BigReal virial_xx = f_net.x * params->offset_f.x;
2189  BigReal virial_xy = f_net.x * params->offset_f.y;
2190  BigReal virial_xz = f_net.x * params->offset_f.z;
2191  BigReal virial_yy = f_net.y * params->offset_f.y;
2192  BigReal virial_yz = f_net.y * params->offset_f.z;
2193  BigReal virial_zz = f_net.z * params->offset_f.z;
2194 
2195  reduction[virialIndex_XX] += virial_xx;
2196  reduction[virialIndex_XY] += virial_xy;
2197  reduction[virialIndex_XZ] += virial_xz;
2198  reduction[virialIndex_YX] += virial_xy;
2199  reduction[virialIndex_YY] += virial_yy;
2200  reduction[virialIndex_YZ] += virial_yz;
2201  reduction[virialIndex_ZX] += virial_xz;
2202  reduction[virialIndex_ZY] += virial_yz;
2203  reduction[virialIndex_ZZ] += virial_zz;
2204  }
2205 #endif
2206 #endif
2207 
2208 #ifdef fullf_1
2209 #undef fullf_1
2210 #endif
2211 #if ( FULL( 1+ ) 0 )
2212 #if ( PAIR( 1+ ) 0 )
2213  {
2214  BigReal fullElectVirial_xx = fullf_net.x * params->offset_f.x;
2215  BigReal fullElectVirial_xy = fullf_net.x * params->offset_f.y;
2216  BigReal fullElectVirial_xz = fullf_net.x * params->offset_f.z;
2217  BigReal fullElectVirial_yy = fullf_net.y * params->offset_f.y;
2218  BigReal fullElectVirial_yz = fullf_net.y * params->offset_f.z;
2219  BigReal fullElectVirial_zz = fullf_net.z * params->offset_f.z;
2220 
2221  reduction[fullElectVirialIndex_XX] += fullElectVirial_xx;
2222  reduction[fullElectVirialIndex_XY] += fullElectVirial_xy;
2223  reduction[fullElectVirialIndex_XZ] += fullElectVirial_xz;
2224  reduction[fullElectVirialIndex_YX] += fullElectVirial_xy;
2225  reduction[fullElectVirialIndex_YY] += fullElectVirial_yy;
2226  reduction[fullElectVirialIndex_YZ] += fullElectVirial_yz;
2227  reduction[fullElectVirialIndex_ZX] += fullElectVirial_xz;
2228  reduction[fullElectVirialIndex_ZY] += fullElectVirial_yz;
2229  reduction[fullElectVirialIndex_ZZ] += fullElectVirial_zz;
2230  }
2231 #endif
2232 #endif
2233 
2234  reduction[exclChecksumIndex] += exclChecksum;
2235  FAST
2236  (
2237  // JLai
2238  ENERGY( reduction[vdwEnergyIndex] += vdwEnergy;
2239  GO( reduction[groLJEnergyIndex] += groLJEnergy;
2240  reduction[groGaussEnergyIndex] += groGaussEnergy;
2241  reduction[goNativeEnergyIndex] += goEnergyNative;
2242  reduction[goNonnativeEnergyIndex] += goEnergyNonnative; ) )
2243  SHORT
2244  (
2245  ENERGY( reduction[electEnergyIndex] += electEnergy; )
2246  )
2247  ALCH
2248  (
2249  TI(reduction[vdwEnergyIndex_ti_1] += vdwEnergy_ti_1;)
2250  TI(reduction[vdwEnergyIndex_ti_2] += vdwEnergy_ti_2;)
2251  FEP( reduction[vdwEnergyIndex_s] += vdwEnergy_s; )
2252  SHORT
2253  (
2254  FEP( reduction[electEnergyIndex_s] += electEnergy_s; )
2255  TI(reduction[electEnergyIndex_ti_1] += electEnergy_ti_1;)
2256  TI(reduction[electEnergyIndex_ti_2] += electEnergy_ti_2;)
2257  )
2258  )
2259  )
2260  FULL
2261  (
2262  ENERGY( reduction[fullElectEnergyIndex] += fullElectEnergy; )
2263  ALCH
2264  (
2265  FEP( reduction[fullElectEnergyIndex_s] += fullElectEnergy_s; )
2266  TI(reduction[fullElectEnergyIndex_ti_1] += fullElectEnergy_ti_1;)
2267  TI(reduction[fullElectEnergyIndex_ti_2] += fullElectEnergy_ti_2;)
2268  )
2269  )
2270 
2271 #ifndef MEM_OPT_VERSION
2272  REZERO_EXCL_FLAGS_BUFF
2273 #undef REZERO_EXCL_FLAGS_BUFF
2274 #endif
2275 
2276 }
2277 
virial_xy
register BigReal virial_xy
Definition: ComputeFullDirectBase.h:19
nonbonded::pairlists
Pairlists * pairlists
Definition: ComputeNonbondedUtil.h:220
sortEntries_mergeSort_v1
void sortEntries_mergeSort_v1(SortEntry *&se, SortEntry *&buf, int seLen)
Definition: ComputeNonbondedInl.h:425
virial_xz
register BigReal virial_xz
Definition: ComputeFullDirectBase.h:20
LJTable::table_row
const TableEntry * table_row(unsigned int i) const
Definition: LJTable.h:31
ResizeArray::size
int size(void) const
Definition: ResizeArray.h:131
PAIR
#define PAIR(X)
Definition: ComputeNonbondedBase.h:58
virial_yz
register BigReal virial_yz
Definition: ComputeFullDirectBase.h:22
Random::uniform
BigReal uniform(void)
Definition: Random.h:109
nonbonded::parameters
Parameters * parameters
Definition: ComputeNonbondedUtil.h:214
BOLTZMANN
#define BOLTZMANN
Definition: common.h:54
sortEntries_selectionSort
void sortEntries_selectionSort(SortEntry *const se, const int seLen)
Definition: ComputeNonbondedInl.h:354
LAM
#define LAM(X)
Definition: ComputeNonbondedBase.h:163
Random::gaussian
BigReal gaussian(void)
Definition: Random.h:116
NOKNL
#define NOKNL(X)
Definition: ComputeNonbondedBase.h:240
partition
static void partition(int *order, const FullAtom *atoms, int begin, int end)
Definition: SortAtoms.C:45
PPROF
#define PPROF(X)
Definition: ComputeNonbondedBase.h:162
nonbonded::p
CompAtom * p[2]
Definition: ComputeNonbondedUtil.h:189
nbthole_pair_value
Definition: Parameters.h:221
Vector
Definition: Vector.h:72
electEnergy
register BigReal electEnergy
Definition: ComputeFullDirectBase.h:17
ComputeNonbondedUtil::mol
static const Molecule * mol
Definition: ComputeNonbondedUtil.h:291
__sort_entry
Definition: NamdTypes.h:305
ExclusionCheck::flags
char * flags
Definition: Molecule.h:72
Real
float Real
Definition: common.h:118
COULOMB
#define COULOMB
Definition: common.h:53
int32
int32_t int32
Definition: common.h:38
nonbonded::v
CompAtom * v[2]
Definition: ComputeNonbondedUtil.h:195
pairlist_from_pairlist
int pairlist_from_pairlist(BigReal cutoff2, BigReal p_i_x, BigReal p_i_y, BigReal p_i_z, const CompAtom *p_j, const plint *list, int list_size, int *newlist, BigReal r2_delta, BigReal *r2list)
Definition: ComputeNonbondedInl.h:80
Vector::z
BigReal z
Definition: Vector.h:74
NBWORKARRAY
#define NBWORKARRAY(TYPE, NAME, SIZE)
Definition: ComputeNonbondedUtil.h:92
CompAtom::position
Position position
Definition: NamdTypes.h:77
NOSHORT
#define NOSHORT(X)
Definition: ComputeNonbondedBase.h:109
nonbonded::reduction
BigReal * reduction
Definition: ComputeNonbondedUtil.h:211
ComputeNonbondedUtil::pressureProfileThickness
static BigReal pressureProfileThickness
Definition: ComputeNonbondedUtil.h:401
Molecule::get_full_exclusions_for_atom
const int32 * get_full_exclusions_for_atom(int anum) const
Definition: Molecule.h:1225
Molecule::GetAtomAlpha
BigReal GetAtomAlpha(int i) const
Definition: Molecule.h:511
p_j
#define p_j
Molecule::get_mod_exclusions_for_atom
const int32 * get_mod_exclusions_for_atom(int anum) const
Definition: Molecule.h:1227
nonbonded::simParameters
SimParameters * simParameters
Definition: ComputeNonbondedUtil.h:215
Molecule::get_excl_check_for_atom
const ExclusionCheck * get_excl_check_for_atom(int anum) const
Definition: Molecule.h:1241
EXCHCK_MOD
#define EXCHCK_MOD
Definition: Molecule.h:87
Molecule
Molecule stores the structural information for the system.
Definition: Molecule.h:175
sortEntries_bubbleSort
void sortEntries_bubbleSort(SortEntry *const se, const int seLen)
Definition: ComputeNonbondedInl.h:388
TI
#define TI(X)
Definition: ComputeNonbondedBase.h:165
ResizeArray::resize
void resize(int i)
Definition: ResizeArray.h:84
CompAtomExt::id
uint32 id
Definition: NamdTypes.h:156
__sort_entry::index
int32 index
Definition: NamdTypes.h:306
INT
#define INT(X)
Definition: ComputeNonbondedBase.h:161
CompAtom::charge
Charge charge
Definition: NamdTypes.h:78
NBWORKARRAYSINIT
#define NBWORKARRAYSINIT(ARRAYS)
Definition: ComputeNonbondedUtil.h:83
NOTABENERGY
#define NOTABENERGY(X)
Definition: ComputeNonbondedBase.h:136
virial_yy
register BigReal virial_yy
Definition: ComputeFullDirectBase.h:21
plint
unsigned short plint
Definition: ComputeNonbondedUtil.h:20
pp_clamp
void pp_clamp(int &n, int nslabs)
Definition: PressureProfile.h:15
Random
Definition: Random.h:37
Pairlists::getIndexValue
plint getIndexValue()
Definition: ComputeNonbondedUtil.h:77
Pairlists::nextlist
void nextlist(plint **list, int *list_size)
Definition: ComputeNonbondedUtil.h:64
Parameters::NumNbtholePairParams
int NumNbtholePairParams
Definition: Parameters.h:343
ComputeNonbondedUtil::table_noshort
static BigReal * table_noshort
Definition: ComputeNonbondedUtil.h:299
nbthole_pair_value::tholeij
Real tholeij
Definition: Parameters.h:227
NAMD_bug
void NAMD_bug(const char *err_msg)
Definition: common.C:195
ComputeNonbondedUtil::pressureProfileMin
static BigReal pressureProfileMin
Definition: ComputeNonbondedUtil.h:402
UNLIKELY
#define UNLIKELY(X)
Definition: ComputeNonbondedBase.h:37
NAME
#define NAME
Definition: ComputeNonbondedBase.h:54
EXCHCK_FULL
#define EXCHCK_FULL
Definition: Molecule.h:86
CompAtom::vdwType
int16 vdwType
Definition: NamdTypes.h:79
SHORT
#define SHORT(X)
Definition: ComputeNonbondedBase.h:108
__sort_entry::sortValue
BigReal sortValue
Definition: NamdTypes.h:307
Bool
int Bool
Definition: common.h:142
SELF
#define SELF(X)
Definition: ComputeNonbondedBase.h:67
ALCH
#define ALCH(X)
Definition: ComputeNonbondedBase.h:164
nonbonded::pExt
CompAtomExt * pExt[2]
Definition: ComputeNonbondedUtil.h:193
CompAtom::partition
uint8 partition
Definition: NamdTypes.h:80
virial_zz
register BigReal virial_zz
Definition: ComputeFullDirectBase.h:23
Vector::x
BigReal x
Definition: Vector.h:74
CompAtom::hydrogenGroupSize
uint8 hydrogenGroupSize
Definition: NamdTypes.h:88
Molecule::get_qmAtomGroup
const Real * get_qmAtomGroup() const
Definition: Molecule.h:853
Molecule::numAtoms
int numAtoms
Definition: Molecule.h:585
FAST
#define FAST(X)
Definition: ComputeNonbondedBase.h:95
Parameters::nbthole_array
NbtholePairValue * nbthole_array
Definition: Parameters.h:321
FEP
#define FEP(X)
Definition: ComputeNonbondedBase.h:157
ComputeNonbondedUtil::lambda_table
static BigReal * lambda_table
Definition: ComputeNonbondedUtil.h:393
ComputeNonbondedUtil::slow_table
static BigReal * slow_table
Definition: ComputeNonbondedUtil.h:302
CompAtom::nonbondedGroupSize
uint8 nonbondedGroupSize
Definition: NamdTypes.h:81
REZERO_EXCL_FLAGS_BUFF
#define REZERO_EXCL_FLAGS_BUFF
simParams
#define simParams
Definition: Output.C:129
KNL
#define KNL(X)
Definition: ComputeNonbondedBase.h:239
nbthole_pair_value::ind1
Index ind1
Definition: Parameters.h:223
Vector::y
BigReal y
Definition: Vector.h:74
ComputeNonbondedUtil::ljTable
static const LJTable * ljTable
Definition: ComputeNonbondedUtil.h:290
ENERGY
#define ENERGY(X)
Definition: ComputeNonbondedBase.h:82
virial_xx
register BigReal virial_xx
Definition: ComputeFullDirectBase.h:18
FULL
#define FULL(X)
Definition: ComputeNonbondedBase.h:122
Pairlists::newsize
void newsize(int list_size)
Definition: ComputeNonbondedUtil.h:59
ComputeNonbondedUtil::table_short
static BigReal * table_short
Definition: ComputeNonbondedUtil.h:298
TIMEFACTOR
#define TIMEFACTOR
Definition: common.h:55
Pairlists::setIndexValue
void setIndexValue(plint i)
Definition: ComputeNonbondedUtil.h:73
Mass
float Mass
Definition: ComputeGBIS.inl:20
COMPONENT_DOTPRODUCT
#define COMPONENT_DOTPRODUCT(A, B)
Definition: ComputeNonbondedBase.h:244
LES
#define LES(X)
Definition: ComputeNonbondedBase.h:160
nonbonded::fullf
Force * fullf[2]
Definition: ComputeNonbondedUtil.h:199
ComputeNonbondedUtil::alchVdwShiftCoeff
static BigReal alchVdwShiftCoeff
Definition: ComputeNonbondedUtil.h:385
Vector::rlength
NAMD_HOST_DEVICE BigReal rlength(void) const
Definition: Vector.h:210
nonbonded::pressureProfileReduction
BigReal * pressureProfileReduction
Definition: ComputeNonbondedUtil.h:212
Pairlists::newlist
plint * newlist(int max_size)
Definition: ComputeNonbondedUtil.h:40
NOENERGY
#define NOENERGY(X)
Definition: ComputeNonbondedBase.h:83
RIGID_NONE
#define RIGID_NONE
Definition: SimParameters.h:79
Vector::unit
NAMD_HOST_DEVICE Vector unit(void) const
Definition: Vector.h:215
nbthole_pair_value::ind2
Index ind2
Definition: Parameters.h:224
sortEntries_mergeSort_v2
void sortEntries_mergeSort_v2(SortEntry *&se, SortEntry *&buf, int seLen)
Definition: ComputeNonbondedInl.h:556
GO
#define GO(X)
Definition: ComputeNonbondedBase.h:166
nonbonded::workArrays
ComputeNonbondedWorkArrays * workArrays
Definition: ComputeNonbondedUtil.h:218
BigReal
double BigReal
Definition: common.h:123
for
for(int i=0;i< n1;++i)
Definition: ComputeFullDirectBase.h:34
TABENERGY
#define TABENERGY(X)
Definition: ComputeNonbondedBase.h:135
Generated by   doxygen 1.8.14