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