NAMD
ComputeNonbondedBase2.h
Go to the documentation of this file.
1 
7 EXCLUDED( FAST( foo bar ) )
8 EXCLUDED( MODIFIED( foo bar ) )
9 EXCLUDED( NORMAL( foo bar ) )
10 NORMAL( MODIFIED( foo bar ) )
11 ALCHPAIR( NOT_ALCHPAIR( foo bar ) )
12 
13 ALCHPAIR(
14  // get alchemical nonbonded scaling parameters (once per pairlist)
15  myLambda = ALCH1(lambdaUp) ALCH2(lambdaDown) ALCH3(lambdaUp) ALCH4(lambdaDown);
16  FEP(myLambda2 = ALCH1(lambda2Up) ALCH2(lambda2Down) ALCH3(lambda2Up) ALCH4(lambda2Down);)
17  myElecLambda = ALCH1(elecLambdaUp) ALCH2(elecLambdaDown) ALCH3(elecLambdaUp) ALCH4(elecLambdaDown);
18  FEP(myElecLambda2 = ALCH1(elecLambda2Up) ALCH2(elecLambda2Down) ALCH3(elecLambda2Up) ALCH4(elecLambda2Down);)
19  myVdwLambda = ALCH1(vdwLambdaUp) ALCH2(vdwLambdaDown) ALCH3(vdwLambdaUp) ALCH4(vdwLambdaDown);
20  FEP(myVdwLambda2 = ALCH1(vdwLambda2Up) ALCH2(vdwLambda2Down) ALCH3(vdwLambda2Up) ALCH4(vdwLambda2Down);)
21  ALCH1(myRepLambda = repLambdaUp) ALCH2(myRepLambda = repLambdaDown);
22  FEP(ALCH1(myRepLambda2 = repLambda2Up) ALCH2(myRepLambda2 = repLambda2Down);)
23  ALCH1(myVdwShift = vdwShiftUp) ALCH2(myVdwShift = vdwShiftDown);
24  FEP(ALCH1(myVdwShift2 = vdwShift2Up) ALCH2(myVdwShift2 = vdwShift2Down);)
25 )
26 
27 #ifdef ARCH_POWERPC
28  __alignx(64, table_four);
29  __alignx(32, p_1);
30 #pragma unroll(1)
31 #pragma ibm independent_loop
32 #endif
33 
34 #ifndef ARCH_POWERPC
35 #ifndef __INTEL_LLVM_COMPILER
36 #pragma ivdep
37 #endif
38 #endif
39 
40 
41 
42 #if ( FULL( EXCLUDED( SHORT( 1+ ) ) ) 0 )
43 // avoid bug in Intel 15.0 compiler
44 #pragma novector
45 #else
46 #ifdef PRAGMA_SIMD
47 #ifndef TABENERGYFLAG
48 #ifndef GOFORCES
49 #pragma omp simd SHORT(FAST(reduction(+:f_i_x,f_i_y,f_i_z)) ENERGY(FAST(reduction(+:vdwEnergy) SHORT(reduction(+:electEnergy))))) \
50  FULL(reduction(+:fullf_i_x,fullf_i_y,fullf_i_z) ENERGY(reduction(+:fullElectEnergy)))
51 #endif
52 #endif
53 #pragma loop_count avg=100
54 #else // PRAGMA_SIMD
55 #pragma loop_count avg=4
56 #endif // PRAGMA_SIMD
57 #endif
58  for (k=0; k<npairi; ++k) {
59  TABENERGY(
60  const int numtypes = simParams->tableNumTypes;
61  const float table_spacing = simParams->tableSpacing;
62  const int npertype = (int) (namdnearbyint(simParams->tableMaxDist / simParams->tableSpacing) + 1);
63  )
64 
65  int table_i = (r2iilist[2*k] >> 14) + r2_delta_expc; // table_i >= 0
66  const int j = pairlisti[k];
67  //register const CompAtom *p_j = p_1 + j;
68 #define p_j (p_1+j)
69  // const CompAtomExt *pExt_j = pExt_1 + j;
70 
71  BigReal diffa = r2list[k] - r2_table[table_i];
72  //const BigReal* const table_four_i = table_four + 16*table_i;
73 #define table_four_i (table_four + 16*table_i)
74 
75 #if ( FAST( 1 + ) TABENERGY( 1 + ) 0 ) // FAST or TABENERGY
76  //const LJTable::TableEntry * lj_pars =
77  // lj_row + 2 * p_j->vdwType MODIFIED(+ 1);
78  const int lj_index = 2 * p_j->vdwType MODIFIED(+ 1);
79 #define lj_pars (lj_row+lj_index)
80 #endif
81 
82  TABENERGY(
83  register const int tabtype = -1 - ( lj_pars->A < 0 ? lj_pars->A : 0 );
84  )
85 
86 #if ( SHORT( FAST( 1+ ) ) 0 )
87  //Force *f_j = f_1 + j;
88 #define f_j (f_1+j)
89 #endif
90 
91 #if ( FULL( 1+ ) 0 )
92  //Force *fullf_j = fullf_1 + j;
93 #define fullf_j (fullf_1+j)
94 #endif
95 
96  //Power PC aliasing and alignment constraints
97 #ifdef ARCH_POWERPC
98 #if ( FULL( 1+ ) 0 )
99 #pragma disjoint (*table_four, *fullf_1)
100 #pragma disjoint (*p_1, *fullf_1)
101 #pragma disjoint (*r2_table, *fullf_1)
102 #pragma disjoint (*r2list, *fullf_1)
103 #if ( SHORT( FAST( 1+ ) ) 0 )
104 #pragma disjoint (*f_1 , *fullf_1)
105 #pragma disjoint (*fullf_1, *f_1)
106 #endif //Short + fast
107 #endif //Full
108 
109 #if ( SHORT( FAST( 1+ ) ) 0 )
110 #pragma disjoint (*table_four, *f_1)
111 #pragma disjoint (*p_1, *f_1)
112 #pragma disjoint (*r2_table, *f_1)
113 #pragma disjoint (*r2list, *f_1)
114 #pragma disjoint (*lj_row, *f_1)
115 #endif //Short + Fast
116 
117  __alignx(64, table_four_i);
118  FAST (
119  __alignx(32, lj_pars);
120  )
121  __alignx(32, p_j);
122 #endif //ARCH_POWERPC
123 
124  /*
125  BigReal modf = 0.0;
126  int atom2 = p_j->id;
127  register char excl_flag = ( (atom2 >= excl_min && atom2 <= excl_max) ?
128  excl_flags[atom2-excl_min] : 0 );
129  if ( excl_flag ) { ++exclChecksum; }
130  SELF( if ( j < j_hgroup ) { excl_flag = EXCHCK_FULL; } )
131  if ( excl_flag ) {
132  if ( excl_flag == EXCHCK_FULL ) {
133  lj_pars = lj_null_pars;
134  modf = 1.0;
135  } else {
136  ++lj_pars;
137  modf = modf_mod;
138  }
139  }
140  */
141 
142  BigReal kqq = kq_i * p_j->charge;
143 
144 
145  LES( BigReal lambda_pair = lambda_table_i[p_j->partition]; )
146 
147  register const BigReal p_ij_x = p_i_x - p_j->position.x;
148  register const BigReal p_ij_y = p_i_y - p_j->position.y;
149  register const BigReal p_ij_z = p_i_z - p_j->position.z;
150 
151 #if ( FAST(1+) 0 )
152  const BigReal A = scaling * lj_pars->A;
153  const BigReal B = scaling * lj_pars->B;
154  BigReal vdw_d = A * table_four_i[0] - B * table_four_i[4];
155  BigReal vdw_c = A * table_four_i[1] - B * table_four_i[5];
156  BigReal vdw_b = A * table_four_i[2] - B * table_four_i[6];
157  BigReal vdw_a = A * table_four_i[3] - B * table_four_i[7];
158 
159  ALCHPAIR (
160  // Alchemical free energy calculation
161  // Pairlist 1 and 2 are for softcore atoms, while 3 and 4 are single topology atoms.
162  // Pairlists are separated so that lambda-coupled pairs are handled
163  // independently from normal nonbonded (inside ALCHPAIR macro).
164  // The separation-shifted van der Waals potential and a shifted
165  // electrostatics potential for decoupling are calculated explicitly.
166  // Would be faster with lookup tables but because only a small minority
167  // of nonbonded pairs are lambda-coupled the impact is minimal.
168  // Explicit calculation also makes things easier to modify.
169  // These are now inline functions (in ComputeNonbondedFep.C) to
170  // tidy the code
171 
172  const BigReal r2 = r2list[k] - r2_delta;
173 
174  // These are now inline functions (in ComputeNonbondedFep.C) to
175  // tidy the code
176 
177  FEP(
178  ALCH1 ( // Don't merge/recombine the ALCH 1, 2, 3 ,4. Their functions might be modified for future algorithm changes.
179  fep_vdw_forceandenergies(A, B, r2, myVdwShift, myVdwShift2,
180  switchdist2, cutoff2, switchfactor, vdwForceSwitching, myVdwLambda,
181  myVdwLambda2, alchWCAOn, myRepLambda, myRepLambda2, &alch_vdw_energy,
182  &alch_vdw_force, &alch_vdw_energy_2);)
183  ALCH2 (
184  fep_vdw_forceandenergies(A, B, r2, myVdwShift, myVdwShift2,
185  switchdist2, cutoff2, switchfactor, vdwForceSwitching, myVdwLambda,
186  myVdwLambda2, alchWCAOn, myRepLambda, myRepLambda2, &alch_vdw_energy,
187  &alch_vdw_force, &alch_vdw_energy_2);)
188  ALCH3 ( // In single topology region ALCH3 & 4, all atoms are paired so softcore potential is unnecessary.
189  ENERGY(alch_vdw_energy = -myVdwLambda * (( ( diffa * vdw_d * (1/6.)+ vdw_c * (1/4.)) * diffa + vdw_b *(1/2.)) * diffa + vdw_a);)
190  alch_vdw_energy_2 = -myVdwLambda2 * (( ( diffa * vdw_d * (1/6.)+ vdw_c * (1/4.)) * diffa + vdw_b *(1/2.)) * diffa + vdw_a);
191  alch_vdw_force = myVdwLambda * ((diffa * vdw_d + vdw_c) * diffa + vdw_b);)
192  ALCH4 (
193  ENERGY(alch_vdw_energy = -myVdwLambda * (( ( diffa * vdw_d * (1/6.)+ vdw_c * (1/4.)) * diffa + vdw_b *(1/2.)) * diffa + vdw_a);)
194  alch_vdw_energy_2 = -myVdwLambda2 * (( ( diffa * vdw_d * (1/6.)+ vdw_c * (1/4.)) * diffa + vdw_b *(1/2.)) * diffa + vdw_a);
195  alch_vdw_force = myVdwLambda * ((diffa * vdw_d + vdw_c) * diffa + vdw_b);)
196  )
197  TI(ti_vdw_force_energy_dUdl(A, B, r2, myVdwShift, switchdist2, cutoff2,
198  switchfactor, vdwForceSwitching, myVdwLambda, alchVdwShiftCoeff,
199  alchWCAOn, myRepLambda, &alch_vdw_energy, &alch_vdw_force,
200  &alch_vdw_dUdl);)
201  )
202 
203  //NOT_ALCHPAIR(
204  //TABENERGY(
205 #if (NOT_ALCHPAIR(1+) 0)
206 #if (TABENERGY(1+) 0)
207  if (tabtype >= 0) {
208  register BigReal r1;
209  r1 = sqrt(p_ij_x*p_ij_x + p_ij_y*p_ij_y + p_ij_z*p_ij_z);
210 
211  //CkPrintf("%i %i %f %f %i\n", npertype, tabtype, r1, table_spacing, (int) (namdnearbyint(r1 / table_spacing)));
212  register int eneraddress;
213  eneraddress = 2 * ((npertype * tabtype) + ((int) namdnearbyint(r1 / table_spacing)));
214  //CkPrintf("Using distance bin %i for distance %f\n", eneraddress, r1);
215  vdw_d = 0.;
216  vdw_c = 0.;
217  vdw_b = table_ener[eneraddress + 1] / r1;
218  vdw_a = (-1/2.) * diffa * vdw_b;
219  ENERGY(
220  register BigReal vdw_val = table_ener[eneraddress];
221  //CkPrintf("Found vdw energy of %f\n", vdw_val);
222  vdwEnergy += LAM(lambda_pair *) vdw_val;
223  FEP( vdwEnergy_s += d_lambda_pair * vdw_val; )
224  )
225  } else {
226  //)
227 #endif
228  ENERGY(
229  register BigReal vdw_val =
230  ( ( diffa * vdw_d * (1/6.)+ vdw_c * (1/4.)) * diffa + vdw_b *(1/2.)) * diffa + vdw_a;
231 
232  vdwEnergy -= LAM(lambda_pair *) vdw_val;
233 
234  FEP(vdwEnergy_s -= vdw_val;)
235  )
236  //TABENERGY( } ) /* endif (tabtype >= 0) */
237 #if (TABENERGY (1+) 0)
238  }
239 #endif
240  //) // NOT_ALCHPAIR
241 #endif
242 
243  ALCHPAIR(
244  ENERGY(vdwEnergy += alch_vdw_energy;)
245  FEP(vdwEnergy_s += alch_vdw_energy_2;)
246  TI(ALCH1(vdwEnergy_ti_1 += alch_vdw_dUdl;) ALCH2(vdwEnergy_ti_2 += alch_vdw_dUdl;))
247  ) // ALCHPAIR
248 
249 #endif // FAST
250 
251 #if ( FAST(1+) 0 )
252  INT(
253  register BigReal vdw_dir;
254  vdw_dir = ( diffa * vdw_d + vdw_c ) * diffa + vdw_b;
255  //BigReal force_r = LAM(lambda_pair *) vdw_dir;
256  reduction[pairVDWForceIndex_X] += force_sign * vdw_dir * p_ij_x;
257  reduction[pairVDWForceIndex_Y] += force_sign * vdw_dir * p_ij_y;
258  reduction[pairVDWForceIndex_Z] += force_sign * vdw_dir * p_ij_z;
259  )
260 
261 #if ( SHORT(1+) 0 ) // Short-range electrostatics
262 
263  NORMAL(
264  BigReal fast_d = kqq * table_four_i[8];
265  BigReal fast_c = kqq * table_four_i[9];
266  BigReal fast_b = kqq * table_four_i[10];
267  BigReal fast_a = kqq * table_four_i[11];
268  )
269  MODIFIED(
270  BigReal modfckqq = (1.0-modf_mod) * kqq;
271  BigReal fast_d = modfckqq * table_four_i[8];
272  BigReal fast_c = modfckqq * table_four_i[9];
273  BigReal fast_b = modfckqq * table_four_i[10];
274  BigReal fast_a = modfckqq * table_four_i[11];
275  )
276 
277  {
278  ENERGY(
279  register BigReal fast_val =
280  ( ( diffa * fast_d * (1/6.)+ fast_c * (1/4.)) * diffa + fast_b *(1/2.)) * diffa + fast_a;
281  NOT_ALCHPAIR (
282  electEnergy -= LAM(lambda_pair *) fast_val;
283  FEP(electEnergy_s -= fast_val;)
284  )
285  ) //ENERGY
286  ALCHPAIR(
287  ENERGY(electEnergy -= myElecLambda * fast_val;)
288  FEP(electEnergy_s -= myElecLambda2 * fast_val;)
289  TI(
290  NOENERGY(register BigReal fast_val =
291  ( ( diffa * fast_d * (1/6.)+ fast_c * (1/4.)) * diffa + fast_b *(1/2.)) * diffa + fast_a;)
292  ALCH1(electEnergy_ti_1 -= fast_val;) ALCH2(electEnergy_ti_2 -= fast_val;)
293  )
294  )
295 
296  INT(
297  register BigReal fast_dir =
298  ( diffa * fast_d + fast_c ) * diffa + fast_b;
299  // force_r -= -1.0 * LAM(lambda_pair *) fast_dir;
300  reduction[pairElectForceIndex_X] += force_sign * fast_dir * p_ij_x;
301  reduction[pairElectForceIndex_Y] += force_sign * fast_dir * p_ij_y;
302  reduction[pairElectForceIndex_Z] += force_sign * fast_dir * p_ij_z;
303  )
304  }
305 
306 
307  /***** JE - Go *****/
308  // Now Go energy should appear in VDW place -- put vdw_b back into place
309 #if ( NORMAL (1+) 0)
310 #if ( GO (1+) 0)
311 
312 
313 // JLai
314 #ifndef CODE_REDUNDANT
315 #define CODE_REDUNDANT 0
316 #endif
317 #if CODE_REDUNDANT
319  // Explicit goGroPair calculation; only calculates goGroPair if goGroPair is turned on
320  //
321  // get_gro_force has an internal checklist that sees if atom_i and atom_j are
322  // in the explicit pairlist. This is done because there is no guarantee that a
323  // processor will have atom_i and atom_j so we cannot loop over the explict atom pairs.
324  // We can only loop over all pairs.
325  //
326  // NOTE: It does not look like fast_b is not normalized by the r vector.
327  //
328  // JLai
329  BigReal groLJe = 0.0;
330  BigReal groGausse = 0.0;
331  const CompAtomExt *pExt_z = pExt_1 + j;
332  BigReal groForce = mol->get_gro_force2(p_ij_x, p_ij_y, p_ij_z,pExt_i.id,pExt_z->id,&groLJe,&groGausse);
333  NAMD_die("Failsafe. This line should never be reached\n");
334  fast_b += groForce;
335  ENERGY(
336  NOT_ALCHPAIR (
337  // JLai
338  groLJEnergy += groLJe;
339  groGaussEnergy += groGausse;
340  )
341  ) //ENERGY
342  }
343 #endif
344  BigReal goNative = 0;
345  BigReal goNonnative = 0;
346  BigReal goForce = 0;
347  register const CompAtomExt *pExt_j = pExt_1 + j;
349  goForce = mol->get_go_force2(p_ij_x, p_ij_y, p_ij_z, pExt_i.id, pExt_j->id,&goNative,&goNonnative);
350  } else {
351  // Ported by JLai -- JE - added (
352  const BigReal r2go = square(p_ij_x, p_ij_y, p_ij_z);
353  const BigReal rgo = sqrt(r2go);
354 
356  goForce = mol->get_go_force(rgo, pExt_i.id, pExt_j->id, &goNative, &goNonnative);
357  } else if (ComputeNonbondedUtil::goMethod == 3) {
358  goForce = mol->get_go_force_new(rgo, pExt_i.id, pExt_j->id, &goNative, &goNonnative);
359  } else {
360  NAMD_die("I SHOULDN'T BE HERE. DYING MELODRAMATICALLY.\n");
361  }
362  }
363 
364  fast_b += goForce;
365  {
366  ENERGY(
367  NOT_ALCHPAIR (
368  // JLai
369  goEnergyNative += goNative;
370  goEnergyNonnative += goNonnative;
371  )
372  ) //ENERGY
373  INT(
374  reduction[pairVDWForceIndex_X] += force_sign * goForce * p_ij_x;
375  reduction[pairVDWForceIndex_Y] += force_sign * goForce * p_ij_y;
376  reduction[pairVDWForceIndex_Z] += force_sign * goForce * p_ij_z;
377  )
378  }
379  // End of INT
380 
381  //DebugM(3,"rgo:" << rgo << ", pExt_i.id:" << pExt_i.id << ", pExt_j->id:" << pExt_j->id << \
382  // ", goForce:" << goForce << ", fast_b:" << fast_b << std::endl);
383 #endif // ) // End of GO macro
384  /***** JE - End Go *****/
385  // End of port JL
386 #endif //) // End of Normal MACRO
387 
388  // Combined short-range electrostatics and VdW force:
389 #if ( NOT_ALCHPAIR(1+) 0)
390  fast_d += vdw_d;
391  fast_c += vdw_c;
392  fast_b += vdw_b;
393  fast_a += vdw_a; // not used!
394 #endif
395 
396  register BigReal fast_dir =
397  (diffa * fast_d + fast_c) * diffa + fast_b;
398 
399  BigReal force_r = LAM(lambda_pair *) fast_dir;
400  ALCHPAIR(
401  force_r *= myElecLambda;
402  force_r += alch_vdw_force;
403  // special ALCH forces already multiplied by relevant lambda
404  )
405 
406  register BigReal tmp_x = force_r * p_ij_x;
407  f_i_x += tmp_x;
408  f_j->x -= tmp_x;
409 
410  register BigReal tmp_y = force_r * p_ij_y;
411  f_i_y += tmp_y;
412  f_j->y -= tmp_y;
413 
414  register BigReal tmp_z = force_r * p_ij_z;
415  f_i_z += tmp_z;
416  f_j->z -= tmp_z;
417 
418  PPROF(
419  const BigReal p_j_z = p_j->position.z;
420  int n2 = (int)floor((p_j_z-pressureProfileMin)*invThickness);
421  pp_clamp(n2, pressureProfileSlabs);
422  int p_j_partition = p_j->partition;
423 
424  pp_reduction(pressureProfileSlabs, n1, n2,
425  p_i_partition, p_j_partition, pressureProfileAtomTypes,
426  tmp_x*p_ij_x, tmp_y * p_ij_y, tmp_z*p_ij_z,
427  pressureProfileReduction);
428  )
429 
430 #endif // SHORT
431 #endif // FAST
432 
433 #if ( FULL (EXCLUDED( SHORT ( 1+ ) ) ) 0 )
434  //const BigReal* const slow_i = slow_table + 4*table_i;
435 #define slow_i (slow_table + 4*table_i)
436 
437 #ifdef ARCH_POWERPC //Alignment and aliasing constraints
438  __alignx (32, slow_table);
439 #if ( SHORT( FAST( 1+ ) ) 0 )
440 #pragma disjoint (*slow_table, *f_1)
441 #endif
442 #pragma disjoint (*slow_table, *fullf_1)
443 #endif //ARCH_POWERPC
444 
445 #endif //FULL
446 
447 
448 #if ( FULL (MODIFIED( SHORT ( 1+ ) ) ) 0 )
449  //const BigReal* const slow_i = slow_table + 4*table_i;
450 #define slow_i (slow_table + 4*table_i)
451 
452 #ifdef ARCH_POWERPC //Alignment and aliasing constraints
453  __alignx (32, slow_table);
454 #if ( SHORT( FAST( 1+ ) ) 0 )
455 #pragma disjoint (*slow_table, *f_1)
456 #endif
457 #pragma disjoint (*slow_table, *fullf_1)
458 #endif //ARCH_POWERPC
459 
460 #endif //FULL
461 
462 #if ( FULL( 1+ ) 0 )
463  BigReal slow_d = table_four_i[8 SHORT(+ 4)];
464  BigReal slow_c = table_four_i[9 SHORT(+ 4)];
465  BigReal slow_b = table_four_i[10 SHORT(+ 4)];
466  BigReal slow_a = table_four_i[11 SHORT(+ 4)];
467  EXCLUDED(
468  SHORT(
469  slow_a += slow_i[3];
470  slow_b += 2.*slow_i[2];
471  slow_c += 4.*slow_i[1];
472  slow_d += 6.*slow_i[0];
473  )
474  NOSHORT(
475  slow_d -= table_four_i[12];
476  slow_c -= table_four_i[13];
477  slow_b -= table_four_i[14];
478  slow_a -= table_four_i[15];
479  )
480  )
481  MODIFIED(
482  SHORT(
483  slow_a += modf_mod * slow_i[3];
484  slow_b += 2.*modf_mod * slow_i[2];
485  slow_c += 4.*modf_mod * slow_i[1];
486  slow_d += 6.*modf_mod * slow_i[0];
487  )
488  NOSHORT(
489  slow_d -= modf_mod * table_four_i[12];
490  slow_c -= modf_mod * table_four_i[13];
491  slow_b -= modf_mod * table_four_i[14];
492  slow_a -= modf_mod * table_four_i[15];
493  )
494  )
495  slow_d *= kqq;
496  slow_c *= kqq;
497  slow_b *= kqq;
498  slow_a *= kqq;
499 
500  ENERGY(
501  register BigReal slow_val =
502  ( ( diffa * slow_d *(1/6.)+ slow_c * (1/4.)) * diffa + slow_b *(1/2.)) * diffa + slow_a;
503 
504  NOT_ALCHPAIR (
505  fullElectEnergy -= LAM(lambda_pair *) slow_val;
506  FEP(fullElectEnergy_s -= slow_val;)
507  )
508  ) // ENERGY
509 
510  ALCHPAIR(
511  ENERGY(fullElectEnergy -= myElecLambda * slow_val;)
512  FEP(fullElectEnergy_s -= myElecLambda2 * slow_val;)
513  TI(
514  NOENERGY(register BigReal slow_val =
515  ( ( diffa * slow_d *(1/6.)+ slow_c * (1/4.)) * diffa + slow_b *(1/2.)) * diffa + slow_a;)
516  ALCH1(fullElectEnergy_ti_1 -= slow_val;) ALCH2(fullElectEnergy_ti_2 -= slow_val;)
517  )
518  )
519 
520  INT( {
521  register BigReal slow_dir =
522  ( diffa * slow_d + slow_c ) * diffa + slow_b;
523  reduction[pairElectForceIndex_X] += force_sign * slow_dir * p_ij_x;
524  reduction[pairElectForceIndex_Y] += force_sign * slow_dir * p_ij_y;
525  reduction[pairElectForceIndex_Z] += force_sign * slow_dir * p_ij_z;
526  } )
527 
528 
529 #if (NOT_ALCHPAIR (1+) 0)
530 #if (FAST(1+) 0)
531 #if (NOSHORT(1+) 0)
532  slow_d += vdw_d;
533  slow_c += vdw_c;
534  slow_b += vdw_b;
535  slow_a += vdw_a; // unused!
536 #endif
537 #endif
538 #endif
539 
540  register BigReal slow_dir = (diffa * slow_d + slow_c) * diffa + slow_b;
541  BigReal fullforce_r = slow_dir LAM(* lambda_pair);
542  ALCHPAIR (
543  fullforce_r *= myElecLambda;
544  FAST( NOSHORT(
545  fullforce_r += alch_vdw_force;
546  ))
547  )
548 
549  {
550  register BigReal ftmp_x = fullforce_r * p_ij_x;
551  fullf_i_x += ftmp_x;
552  fullf_j->x -= ftmp_x;
553  register BigReal ftmp_y = fullforce_r * p_ij_y;
554  fullf_i_y += ftmp_y;
555  fullf_j->y -= ftmp_y;
556  register BigReal ftmp_z = fullforce_r * p_ij_z;
557  fullf_i_z += ftmp_z;
558  fullf_j->z -= ftmp_z;
559 
560  PPROF(
561  const BigReal p_j_z = p_j->position.z;
562  int n2 = (int)floor((p_j_z-pressureProfileMin)*invThickness);
563  pp_clamp(n2, pressureProfileSlabs);
564  int p_j_partition = p_j->partition;
565 
566  pp_reduction(pressureProfileSlabs, n1, n2,
567  p_i_partition, p_j_partition, pressureProfileAtomTypes,
568  ftmp_x*p_ij_x, ftmp_y * p_ij_y, ftmp_z*p_ij_z,
569  pressureProfileReduction);
570 
571  )
572  }
573 #endif //FULL
574 
575  } // for pairlist
576 
577 #undef p_j
578 #undef lj_pars
579 #undef table_four_i
580 #undef slow_i
581 #undef f_j
582 #undef fullf_j
583 
#define NORMAL(X)
#define namdnearbyint(x)
Definition: common.h:85
#define f_j
#define LAM(X)
#define PPROF(X)
void pp_reduction(int nslabs, int n1, int n2, int atype1, int atype2, int numtypes, BigReal vxx, BigReal vyy, BigReal vzz, BigReal *reduction)
register BigReal electEnergy
#define NOSHORT(X)
#define lj_pars
#define TI(X)
uint32 id
Definition: NamdTypes.h:156
register const BigReal p_ij_z
#define INT(X)
void pp_clamp(int &n, int nslabs)
void fep_vdw_forceandenergies(BigReal A, BigReal B, BigReal r2, BigReal myVdwShift, BigReal myVdwShift2, BigReal switchdist2, BigReal cutoff2, BigReal switchfactor, Bool vdwForceSwitching, BigReal myVdwLambda, BigReal myVdwLambda2, Bool alchWCAOn, BigReal myRepLambda, BigReal myRepLambda2, BigReal *alch_vdw_energy, BigReal *alch_vdw_force, BigReal *alch_vdw_energy_2)
#define SHORT(X)
#define fullf_j
#define NOT_ALCHPAIR(X)
#define p_j
void NAMD_die(const char *err_msg)
Definition: common.C:147
#define table_four_i
#define FAST(X)
void ti_vdw_force_energy_dUdl(BigReal A, BigReal B, BigReal r2, BigReal myVdwShift, BigReal switchdist2, BigReal cutoff2, BigReal switchfactor, Bool vdwForceSwitching, BigReal myVdwLambda, BigReal alchVdwShiftCoeff, Bool alchWCAOn, BigReal myRepLambda, BigReal *alch_vdw_energy, BigReal *alch_vdw_force, BigReal *alch_vdw_dUdl)
#define FEP(X)
TABENERGY(register const int tabtype=-1 -(lj_pars->A< 0 ? lj_pars->A :0);) BigReal kqq
register const BigReal p_ij_x
#define simParams
Definition: Output.C:129
#define ENERGY(X)
#define FULL(X)
#define LES(X)
k< npairi;++k) { TABENERGY(const int numtypes=simParams->tableNumTypes;const float table_spacing=simParams->tableSpacing;const int npertype=(int)(namdnearbyint(simParams->tableMaxDist/simParams->tableSpacing)+1);) int table_i=(r2iilist[2 *k] >> 14)+r2_delta_expc;const int j=pairlisti[k];#define p_j BigReal diffa=r2list[k] - r2_table[table_i];#define table_four_i const int lj_index=2 *p_j-> vdwType MODIFIED(+1)
#define NOENERGY(X)
ALCHPAIR(myLambda=ALCH1(lambdaUp) ALCH2(lambdaDown) ALCH3(lambdaUp) ALCH4(lambdaDown);FEP(myLambda2=ALCH1(lambda2Up) ALCH2(lambda2Down) ALCH3(lambda2Up) ALCH4(lambda2Down);) myElecLambda=ALCH1(elecLambdaUp) ALCH2(elecLambdaDown) ALCH3(elecLambdaUp) ALCH4(elecLambdaDown);FEP(myElecLambda2=ALCH1(elecLambda2Up) ALCH2(elecLambda2Down) ALCH3(elecLambda2Up) ALCH4(elecLambda2Down);) myVdwLambda=ALCH1(vdwLambdaUp) ALCH2(vdwLambdaDown) ALCH3(vdwLambdaUp) ALCH4(vdwLambdaDown);FEP(myVdwLambda2=ALCH1(vdwLambda2Up) ALCH2(vdwLambda2Down) ALCH3(vdwLambda2Up) ALCH4(vdwLambda2Down);) ALCH1(myRepLambda=repLambdaUp) ALCH2(myRepLambda=repLambdaDown);FEP(ALCH1(myRepLambda2=repLambda2Up) ALCH2(myRepLambda2=repLambda2Down);) ALCH1(myVdwShift=vdwShiftUp) ALCH2(myVdwShift=vdwShiftDown);FEP(ALCH1(myVdwShift2=vdwShift2Up) ALCH2(myVdwShift2=vdwShift2Down);)) for(k=0
register const BigReal p_ij_y
double BigReal
Definition: common.h:123
#define EXCLUDED(X)