1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | |
10 | |
11 | |
12 | |
13 | |
14 | |
15 | |
16 | |
17 | |
18 | |
19 | |
20 | |
21 | |
22 | |
23 | |
24 | |
25 | |
26 | |
27 | |
28 | |
29 | |
30 | |
31 | |
32 | |
33 | |
34 | |
35 | |
36 | |
37 | |
38 | |
39 | #if defined CHECK_EXCLS && (defined CALC_COULOMB || defined LJ_EWALD) |
40 | #define EXCL_FORCES |
41 | #endif |
42 | |
43 | { |
44 | int cj; |
45 | #ifdef ENERGY_GROUPS |
46 | int egp_cj; |
47 | #endif |
48 | int i; |
49 | |
50 | cj = l_cj[cjind].cj; |
51 | |
52 | #ifdef ENERGY_GROUPS |
53 | egp_cj = nbat->energrp[cj]; |
54 | #endif |
55 | for (i = 0; i < UNROLLI; i++) |
56 | { |
57 | int ai; |
58 | int type_i_off; |
59 | int j; |
60 | |
61 | ai = ci*UNROLLI + i; |
62 | |
63 | type_i_off = type[ai]*ntype2; |
64 | |
65 | for (j = 0; j < UNROLLJ; j++) |
66 | { |
67 | int aj; |
68 | real dx, dy, dz; |
69 | real rsq, rinv; |
70 | real rinvsq, rinvsix; |
71 | real c6, c12; |
72 | real FrLJ6 = 0, FrLJ12 = 0, frLJ = 0, VLJ = 0; |
73 | #if defined LJ_FORCE_SWITCH || defined LJ_POT_SWITCH |
74 | real r, rsw; |
75 | #endif |
76 | |
77 | #ifdef CALC_COULOMB |
78 | real qq; |
79 | real fcoul; |
80 | #ifdef CALC_COUL_TAB |
81 | real rs, frac; |
82 | int ri; |
83 | real fexcl; |
84 | #endif |
85 | #ifdef CALC_ENERGIES |
86 | real vcoul; |
87 | #endif |
88 | #endif |
89 | real fscal; |
90 | real fx, fy, fz; |
91 | |
92 | |
93 | |
94 | |
95 | |
96 | real skipmask; |
97 | |
98 | #ifdef CHECK_EXCLS |
99 | |
100 | |
101 | |
102 | int interact; |
103 | |
104 | interact = ((l_cj[cjind].excl>>(i*UNROLLI + j)) & 1); |
105 | #ifndef EXCL_FORCES |
106 | skipmask = interact; |
107 | #else |
108 | skipmask = !(cj == ci_sh && j <= i); |
109 | #endif |
110 | #else |
111 | #define interact 1.0 |
112 | skipmask = 1.0; |
113 | #endif |
114 | |
115 | aj = cj*UNROLLJ + j; |
116 | |
117 | dx = xi[i*XI_STRIDE+XX0] - x[aj*X_STRIDE+XX0]; |
118 | dy = xi[i*XI_STRIDE+YY1] - x[aj*X_STRIDE+YY1]; |
119 | dz = xi[i*XI_STRIDE+ZZ2] - x[aj*X_STRIDE+ZZ2]; |
120 | |
121 | rsq = dx*dx + dy*dy + dz*dz; |
122 | |
123 | |
124 | skipmask = (rsq >= rcut2) ? 0 : skipmask; |
125 | |
126 | |
127 | #ifdef CHECK_EXCLS |
128 | |
129 | |
130 | |
131 | |
132 | rsq += (1 - interact)*NBNXN_AVOID_SING_R2_INC1.0e-12f; |
133 | #endif |
134 | |
135 | #ifdef COUNT_PAIRS |
136 | npair++; |
137 | #endif |
138 | |
139 | rinv = gmx_invsqrt(rsq)gmx_software_invsqrt(rsq); |
140 | |
141 | |
142 | |
143 | |
144 | |
145 | |
146 | rinv = rinv * skipmask; |
147 | |
148 | rinvsq = rinv*rinv; |
149 | |
150 | #ifdef HALF_LJ |
151 | if (i < UNROLLI/2) |
152 | #endif |
153 | { |
154 | c6 = nbfp[type_i_off+type[aj]*2 ]; |
155 | c12 = nbfp[type_i_off+type[aj]*2+1]; |
156 | |
157 | #if defined LJ_CUT || defined LJ_FORCE_SWITCH || defined LJ_POT_SWITCH |
158 | rinvsix = interact*rinvsq*rinvsq*rinvsq; |
159 | FrLJ6 = c6*rinvsix; |
160 | FrLJ12 = c12*rinvsix*rinvsix; |
161 | frLJ = FrLJ12 - FrLJ6; |
162 | |
163 | #if defined CALC_ENERGIES || defined LJ_POT_SWITCH |
164 | VLJ = (FrLJ12 + c12*ic->repulsion_shift.cpot)/12 - |
165 | (FrLJ6 + c6*ic->dispersion_shift.cpot)/6; |
166 | |
167 | #endif |
168 | #endif |
169 | |
170 | #if defined LJ_FORCE_SWITCH || defined LJ_POT_SWITCH |
171 | |
172 | r = rsq*rinv; |
173 | rsw = r - ic->rvdw_switch; |
174 | rsw = (rsw >= 0.0 ? rsw : 0.0); |
175 | #endif |
176 | #ifdef LJ_FORCE_SWITCH |
177 | frLJ += |
178 | -c6*(ic->dispersion_shift.c2 + ic->dispersion_shift.c3*rsw)*rsw*rsw*r |
179 | + c12*(ic->repulsion_shift.c2 + ic->repulsion_shift.c3*rsw)*rsw*rsw*r; |
180 | #if defined CALC_ENERGIES |
181 | VLJ += |
182 | -c6*(-ic->dispersion_shift.c2/3 - ic->dispersion_shift.c3/4*rsw)*rsw*rsw*rsw |
183 | + c12*(-ic->repulsion_shift.c2/3 - ic->repulsion_shift.c3/4*rsw)*rsw*rsw*rsw; |
184 | #endif |
185 | #endif |
186 | |
187 | #if defined CALC_ENERGIES || defined LJ_POT_SWITCH |
188 | |
189 | |
190 | |
191 | |
192 | VLJ = VLJ * interact; |
193 | #endif |
194 | |
195 | #ifdef LJ_POT_SWITCH |
196 | { |
197 | real sw, dsw; |
198 | |
199 | sw = 1.0 + (swV3 + (swV4+ swV5*rsw)*rsw)*rsw*rsw*rsw; |
200 | dsw = (swF2 + (swF3 + swF4*rsw)*rsw)*rsw*rsw; |
201 | |
202 | frLJ = frLJ*sw - r*VLJ*dsw; |
203 | VLJ *= sw; |
| Value stored to 'VLJ' is never read |
204 | } |
205 | #endif |
206 | |
207 | #ifdef LJ_EWALD |
208 | { |
209 | real c6grid, rinvsix_nm, cr2, expmcr2, poly, sh_mask; |
210 | |
211 | #ifdef LJ_EWALD_COMB_GEOM |
212 | c6grid = ljc[type[ai]*2]*ljc[type[aj]*2]; |
213 | #elif defined LJ_EWALD_COMB_LB |
214 | { |
215 | real sigma, sigma2, epsilon; |
216 | |
217 | |
218 | sigma = ljc[type[ai]*2] + ljc[type[aj]*2]; |
219 | epsilon = ljc[type[ai]*2+1]*ljc[type[aj]*2+1]; |
220 | |
221 | sigma2 = sigma*sigma; |
222 | c6grid = epsilon*sigma2*sigma2*sigma2; |
223 | } |
224 | #else |
225 | #error "No LJ Ewald combination rule defined" |
226 | #endif |
227 | |
228 | #ifdef CHECK_EXCLS |
229 | |
230 | rinvsix_nm = rinvsq*rinvsq*rinvsq; |
231 | #else |
232 | rinvsix_nm = rinvsix; |
233 | #endif |
234 | cr2 = lje_coeff2*rsq; |
235 | #ifdef GMX_DOUBLE |
236 | expmcr2 = exp(-cr2); |
237 | #else |
238 | expmcr2 = expf(-cr2); |
239 | #endif |
240 | poly = 1 + cr2 + 0.5*cr2*cr2; |
241 | |
242 | |
243 | frLJ += c6grid*(rinvsix_nm - expmcr2*(rinvsix_nm*poly + lje_coeff6_6)); |
244 | #ifdef CALC_ENERGIES |
245 | |
246 | sh_mask = lje_vc*interact; |
247 | |
248 | VLJ += c6grid/6*(rinvsix_nm*(1 - expmcr2*poly) + sh_mask); |
249 | #endif |
250 | } |
251 | #endif /* LJ_EWALD */ |
252 | |
253 | #ifdef VDW_CUTOFF_CHECK |
254 | |
255 | { |
256 | real skipmask_rvdw; |
257 | |
258 | skipmask_rvdw = (rsq < rvdw2); |
259 | frLJ *= skipmask_rvdw; |
260 | #ifdef CALC_ENERGIES |
261 | VLJ *= skipmask_rvdw; |
262 | #endif |
263 | } |
264 | #else |
265 | #if defined CALC_ENERGIES |
266 | |
267 | VLJ = VLJ * skipmask; |
268 | |
269 | #endif |
270 | #endif /* VDW_CUTOFF_CHECK */ |
271 | |
272 | |
273 | #ifdef CALC_ENERGIES |
274 | #ifdef ENERGY_GROUPS |
275 | Vvdw[egp_sh_i[i]+((egp_cj>>(nbat->neg_2log*j)) & egp_mask)] += VLJ; |
276 | #else |
277 | Vvdw_ci += VLJ; |
278 | |
279 | #endif |
280 | #endif |
281 | } |
282 | |
283 | #ifdef CALC_COULOMB |
284 | |
285 | |
286 | |
287 | |
288 | |
289 | |
290 | |
291 | |
292 | |
293 | qq = skipmask * qi[i] * q[aj]; |
294 | |
295 | #ifdef CALC_COUL_RF |
296 | fcoul = qq*(interact*rinv*rinvsq - k_rf2); |
297 | |
298 | #ifdef CALC_ENERGIES |
299 | vcoul = qq*(interact*rinv + k_rf*rsq - c_rf); |
300 | |
301 | #endif |
302 | #endif |
303 | |
304 | #ifdef CALC_COUL_TAB |
305 | rs = rsq*rinv*ic->tabq_scale; |
306 | ri = (int)rs; |
307 | frac = rs - ri; |
308 | #ifndef GMX_DOUBLE |
309 | |
310 | fexcl = tab_coul_FDV0[ri*4] + frac*tab_coul_FDV0[ri*4+1]; |
311 | #else |
312 | |
313 | fexcl = (1 - frac)*tab_coul_F[ri] + frac*tab_coul_F[ri+1]; |
314 | #endif |
315 | fcoul = interact*rinvsq - fexcl; |
316 | |
317 | #ifdef CALC_ENERGIES |
318 | #ifndef GMX_DOUBLE |
319 | vcoul = qq*(interact*(rinv - ic->sh_ewald) |
320 | -(tab_coul_FDV0[ri*4+2] |
321 | -halfsp*frac*(tab_coul_FDV0[ri*4] + fexcl))); |
322 | |
323 | #else |
324 | vcoul = qq*(interact*(rinv - ic->sh_ewald) |
325 | -(tab_coul_V[ri] |
326 | -halfsp*frac*(tab_coul_F[ri] + fexcl))); |
327 | #endif |
328 | #endif |
329 | fcoul *= qq*rinv; |
330 | #endif |
331 | |
332 | #ifdef CALC_ENERGIES |
333 | #ifdef ENERGY_GROUPS |
334 | Vc[egp_sh_i[i]+((egp_cj>>(nbat->neg_2log*j)) & egp_mask)] += vcoul; |
335 | #else |
336 | Vc_ci += vcoul; |
337 | |
338 | #endif |
339 | #endif |
340 | #endif |
341 | |
342 | #ifdef CALC_COULOMB |
343 | #ifdef HALF_LJ |
344 | if (i < UNROLLI/2) |
345 | #endif |
346 | { |
347 | fscal = frLJ*rinvsq + fcoul; |
348 | |
349 | } |
350 | #ifdef HALF_LJ |
351 | else |
352 | { |
353 | fscal = fcoul; |
354 | } |
355 | #endif |
356 | #else |
357 | fscal = frLJ*rinvsq; |
358 | #endif |
359 | fx = fscal*dx; |
360 | fy = fscal*dy; |
361 | fz = fscal*dz; |
362 | |
363 | |
364 | fi[i*FI_STRIDE+XX0] += fx; |
365 | fi[i*FI_STRIDE+YY1] += fy; |
366 | fi[i*FI_STRIDE+ZZ2] += fz; |
367 | |
368 | f[aj*F_STRIDE+XX0] -= fx; |
369 | f[aj*F_STRIDE+YY1] -= fy; |
370 | f[aj*F_STRIDE+ZZ2] -= fz; |
371 | |
372 | } |
373 | } |
374 | } |
375 | |
376 | #undef interact |
377 | #undef EXCL_FORCES |