NAMD
ComputeNonbondedBase2KNL.h
Go to the documentation of this file.
1 
7 #ifndef KNL_MAKE_DEPENDS_INCLUDE
8 
9 #if 0
10 #if (VDW_SWITCH_MODE == VDW_SWITCH_MODE_FORCE) || (VDW_SWITCH_MODE == VDW_SWITCH_MODE_MARTINI)
11 // REMOVE WHEN r2list NO LONGER NEEDED!
12 #pragma ivdep
13 for (k=0; k<npairi; ++k) {
14  r2list[k] = r2list_f[k] + r2_delta;
15 }
16 #endif
17 #endif
18 
19 EXCLUDED( FAST( foo bar ) )
20 EXCLUDED( MODIFIED( foo bar ) )
21 EXCLUDED( NORMAL( foo bar ) )
22 NORMAL( MODIFIED( foo bar ) )
23 ALCHPAIR( NOT_ALCHPAIR( foo bar ) )
24 
25 EXCLUDED( foo bar )
26 MODIFIED( foo bar )
27 ALCHPAIR( foo bar )
28 TABENERGY( foo bar )
29 NOFAST( foo bar )
30 
31 #pragma ivdep
32 
33 #if ( FULL( EXCLUDED( SHORT( 1+ ) ) ) 0 )
34 // avoid bug in Intel 15.0 compiler
35 #pragma novector
36 #else
37 #ifdef PRAGMA_SIMD
38 #ifndef TABENERGYFLAG
39 #if __INTEL_COMPILER_BUILD_DATE == 20160721
40 #warning disabled simd pragma on innner loop due to compiler segfault
41 #else
42 #pragma omp simd SHORT(FAST(reduction(+:f_i_x,f_i_y,f_i_z)) ENERGY(FAST(reduction(+:vdwEnergy) SHORT(reduction(+:electEnergy))))) \
43  FULL(reduction(+:fullf_i_x,fullf_i_y,fullf_i_z) ENERGY(reduction(+:fullElectEnergy)))
44 #endif
45 #endif
46 #pragma loop_count avg=100
47 #else // PRAGMA_SIMD
48 #pragma loop_count avg=4
49 #endif // PRAGMA_SIMD
50 #endif
51  for (k=0; k<npairi; ++k) {
52 
53  const float r2 = r2list_f[k];
54  const float r_1 = 1.f / sqrtf(r2);
55  const float r_2 = r_1 * r_1;
56  const float knl_table_r_1 =
57  r_1 > KNL_TABLE_MAX_R_1 ? KNL_TABLE_MAX_R_1 : r_1;
58  const float knl_table_f = (KNL_TABLE_FACTOR-2) * knl_table_r_1;
59  const int knl_table_i = knl_table_f;
60  const float knl_diff = knl_table_f - knl_table_i;
61 
62  const int j = pairlisti[k];
63  //register const CompAtom *p_j = p_1 + j;
64 #define p_j (p_1+j)
65 #define pFlt_j (pFlt_1+j)
66 
67 #if (VDW_SWITCH_MODE == VDW_SWITCH_MODE_FORCE) || (VDW_SWITCH_MODE == VDW_SWITCH_MODE_MARTINI)
68 #if 0
69  int table_i = (r2iilist[2*k] >> 14) + r2_delta_expc; // table_i >= 0
70 
71  float 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 #endif
75 #endif
76 
77  //const LJTable::TableEntry * lj_pars =
78  // lj_row + 2 * p_j->vdwType MODIFIED(+ 1);
79  const int lj_index = 2 * pFlt_j->vdwType MODIFIED(+ 1);
80 #define lj_pars (lj_row+lj_index)
81 
82 #if ( SHORT( 1+ ) 0 )
83  //Force *f_j = f_1 + j;
84 #define f_j (f_1+j)
85 #endif
86 
87 #if ( FULL( 1+ ) 0 )
88  //Force *fullf_j = fullf_1 + j;
89 #define fullf_j (fullf_1+j)
90 #endif
91 
92  float kqq = kq_i_f * p_j->charge;
93 
94  LES( float lambda_pair = lambda_table_i[p_j->partition]; )
95 
96  register const float p_ij_x = xlist[k];
97  register const float p_ij_y = ylist[k];
98  register const float p_ij_z = zlist[k];
99 
100  const float A = scaling_f * lj_pars->A;
101  const float B = scaling_f * lj_pars->B;
102 
103 #if VDW_SWITCH_MODE == VDW_SWITCH_MODE_FORCE
104  { int vdw_switch_mode_force; } // for preprocessor debugging only
105  float vdw_b = 0.f;
106  {
107  const float r_6 = r_2 * r_2 * r_2;
108  float vdwa_energy, vdwb_energy, vdwa_gradient, vdwb_gradient;
109  // from Steinbach & Brooks, JCC 15, pgs 667-683, 1994, eqns 10-13
110  if ( r2 > switchOn2_f ) {
111  const float tmpa = r_6 - cutoff_6_f;
112  vdwa_energy = k_vdwa_f * tmpa * tmpa;
113  const float tmpb = r_1 * r_2 - cutoff_3_f;
114  vdwb_energy = k_vdwb_f * tmpb * tmpb;
115  vdwa_gradient = -6.f * k_vdwa_f * tmpa * r_2 * r_6;
116  vdwb_gradient = -3.f * k_vdwb_f * tmpb * r_2 * r_2 * r_1;
117  } else {
118  const float r_12 = r_6 * r_6;
119  vdwa_energy = r_12 + v_vdwa_f;
120  vdwb_energy = r_6 + v_vdwb_f;
121  vdwa_gradient = -6.f * r_2 * r_12;
122  vdwb_gradient = -3.f * r_2 * r_6;
123  }
124  vdw_b = -2.f * ( A * vdwa_gradient - B * vdwb_gradient );
125  ENERGY(
126  vdwEnergy += A * vdwa_energy - B * vdwb_energy;
127  )
128  }
129 #elif VDW_SWITCH_MODE == VDW_SWITCH_MODE_MARTINI
130  { int vdw_switch_mode_martini; } // for preprocessor debugging only
131  float vdw_b = 0.f;
132  {
133  const float r = r2 * r_1;
134  const float r12 = (r-switchOn_f)*(r-switchOn_f);
135  const float r13 = (r-switchOn_f)*(r-switchOn_f)*(r-switchOn_f);
136 
137  ENERGY(
138  const float LJshifttempA = -(1.f/3.f)*A12_f*r13 - (1.f/4.f)*B12_f*r12*r12 - C12_f;
139  const float LJshifttempB = -(1.f/3.f)*A6_f*r13 - (1.f/4.f)*B6_f*r12*r12 - C6_f;
140  const float shiftValA = ( r > switchOn_f ? LJshifttempA : -C12_f);
141  const float shiftValB = ( r > switchOn_f ? LJshifttempB : -C6_f);
142  )
143 
144  const float LJdshifttempA = -A12_f*r12 - B12_f*r13;
145  const float LJdshifttempB = -A6_f*r12 - B6_f*r13;
146  const float dshiftValA = ( r > switchOn_f ? LJdshifttempA*0.5f*r_1 : 0 );
147  const float dshiftValB = ( r > switchOn_f ? LJdshifttempB*0.5f*r_1 : 0 );
148 
149  const float r_6 = r_2 * r_2 * r_2;
150  const float r_12 = r_6 * r_6;
151 
152  ENERGY(
153  const float vdwa_energy = r_12 + shiftValA;
154  const float vdwb_energy = r_6 + shiftValB;
155  )
156 
157  const float vdwa_gradient = -6.f * r_2 * r_12 + dshiftValA ;
158  const float vdwb_gradient = -3.f * r_2 * r_6 + dshiftValB;
159 
160  vdw_b = -2.f * ( A * vdwa_gradient - B * vdwb_gradient );
161  ENERGY(
162  vdwEnergy += A * vdwa_energy - B * vdwb_energy;
163  )
164  }
165 #elif VDW_SWITCH_MODE == VDW_SWITCH_MODE_ENERGY
166  { int vdw_switch_mode_energy; } // for preprocessor debugging only
167  float vdw_b = 0.f;
168  {
169  const float r_6 = r_2 * r_2 * r_2;
170  const float r_12 = r_6 * r_6;
171  const float c2 = cutoff2_f-r2;
172  const float c4 = c2*(c3_f-2.f*c2);
173  const float switchVal = // used for Lennard-Jones
174  ( r2 > switchOn2_f ? c2*c4*c1_f : 1.f );
175  const float dSwitchVal = // d switchVal / d r2
176  ( r2 > switchOn2_f ? 2.f*c1_f*(c2*c2-c4) : 0.f );
177  const float vdwa_gradient = ( dSwitchVal - 6.f * switchVal * r_2 ) * r_12;
178  const float vdwb_gradient = ( dSwitchVal - 3.f * switchVal * r_2 ) * r_6;
179  vdw_b = -2.f * ( A * vdwa_gradient - B * vdwb_gradient );
180  ENERGY(
181  vdwEnergy += switchVal * ( A * r_12 - B * r_6 );
182  )
183  }
184 #else
185 #error VDW_SWITCH_MODE not recognized
186 #endif // VDW_SWITCH_MODE
187 
188 #if ( SHORT(1+) 0 ) // Short-range electrostatics
189 
190  NORMAL(
191  float fast_b = kqq * ( knl_fast_grad_table[knl_table_i] * (1.f-knl_diff) +
192  knl_fast_grad_table[knl_table_i+1] * knl_diff );
193  )
194 
195  {
196  ENERGY(
197  float fast_val = kqq * ( knl_fast_ener_table[knl_table_i] * (1.f-knl_diff) +
198  knl_fast_ener_table[knl_table_i+1] * knl_diff );
199  electEnergy -= LAM(lambda_pair *) fast_val;
200  ) //ENERGY
201  }
202 
203  // Combined short-range electrostatics and VdW force:
204  fast_b += vdw_b;
205 
206  float fast_dir = fast_b;
207 
208  float force_r = LAM(lambda_pair *) fast_dir;
209 
210  BigReal tmp_x = force_r * p_ij_x;
211  f_i_x += tmp_x;
212  f_j->x -= tmp_x;
213 
214  BigReal tmp_y = force_r * p_ij_y;
215  f_i_y += tmp_y;
216  f_j->y -= tmp_y;
217 
218  BigReal tmp_z = force_r * p_ij_z;
219  f_i_z += tmp_z;
220  f_j->z -= tmp_z;
221 
222 #endif // SHORT
223 
224 #if ( FULL( 1+ ) 0 )
225  #if ( SHORT( 1+ ) 0 )
226  float slow_b = kqq * ( knl_scor_grad_table[knl_table_i] * (1.f-knl_diff) +
227  knl_scor_grad_table[knl_table_i+1] * knl_diff );
228  ENERGY(
229  float slow_val = kqq * ( knl_scor_ener_table[knl_table_i] * (1.f-knl_diff) +
230  knl_scor_ener_table[knl_table_i+1] * knl_diff );
231  )
232  #else
233  float slow_b = kqq * ( knl_corr_grad_table[knl_table_i] * (1.f-knl_diff) +
234  knl_corr_grad_table[knl_table_i+1] * knl_diff );
235  ENERGY(
236  float slow_val = kqq * ( knl_corr_ener_table[knl_table_i] * (1.f-knl_diff) +
237  knl_corr_ener_table[knl_table_i+1] * knl_diff );
238  )
239  #endif
240 
241  ENERGY(
242  fullElectEnergy -= LAM(lambda_pair *) slow_val;
243  ) // ENERGY
244 
245 #if (NOSHORT(1+) 0)
246  slow_b += vdw_b;
247 #endif
248 
249  register float slow_dir = slow_b;
250  float fullforce_r = slow_dir LAM(* lambda_pair);
251 
252  {
253  BigReal ftmp_x = fullforce_r * p_ij_x;
254  fullf_i_x += ftmp_x;
255  fullf_j->x -= ftmp_x;
256  BigReal ftmp_y = fullforce_r * p_ij_y;
257  fullf_i_y += ftmp_y;
258  fullf_j->y -= ftmp_y;
259  BigReal ftmp_z = fullforce_r * p_ij_z;
260  fullf_i_z += ftmp_z;
261  fullf_j->z -= ftmp_z;
262  }
263 #endif //FULL
264 
265  } // for pairlist
266 
267 #undef p_j
268 #undef lj_pars
269 #undef table_four_i
270 #undef slow_i
271 #undef f_j
272 #undef fullf_j
273 
274 #endif // KNL_MAKE_DEPENDS_INCLUDE
275 
#define NORMAL(X)
#define lj_pars
register BigReal fast_dir
BigReal force_r
#define LAM(X)
const BigReal A
#define pFlt_j
register BigReal electEnergy
#define VDW_SWITCH_MODE_MARTINI
#define NOSHORT(X)
BigReal vdw_b
#define ALCHPAIR(X)
#define f_j
register const BigReal p_ij_z
#define p_j
#define NOFAST(X)
#define SHORT(X)
#define NOT_ALCHPAIR(X)
#define FAST(X)
#define MODIFIED(X)
register const BigReal p_ij_x
#define ENERGY(X)
const BigReal B
#define VDW_SWITCH_MODE_ENERGY
#define LES(X)
register const BigReal p_ij_y
double BigReal
Definition: common.h:114
#define EXCLUDED(X)
#define TABENERGY(X)