?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Update instructions in containers.rst
[gromacs.git] / src / gromacs / nbnxm / kernels_simd_4xm / kernel_outer.h
blob19ea6304ba15f17bacf99d28213a7b51261b6983
1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 2012,2013,2014,2015,2016 by the GROMACS development team.
5 * Copyright (c) 2017,2018,2019,2020, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
9 *
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
14 *
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
19 *
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 *
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
32 *
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
35 */
36
37
38 {
39 using namespace gmx;
40
41 /* Unpack pointers for output */
42 real* f = out->f.data();
43 real* fshift = out->fshift.data();
44 #ifdef CALC_ENERGIES
45 # ifdef ENERGY_GROUPS
46 real* Vvdw = out->VSvdw.data();
47 real* Vc = out->VSc.data();
48 # else
49 real* Vvdw = out->Vvdw.data();
50 real* Vc = out->Vc.data();
51 # endif
52 #endif
53
54 const nbnxn_cj_t* l_cj;
55 int ci, ci_sh;
56 int ish, ish3;
57 gmx_bool do_LJ, half_LJ, do_coul;
58 int cjind0, cjind1, cjind;
59
60 #ifdef ENERGY_GROUPS
61 int Vstride_i;
62 int egps_ishift, egps_imask;
63 int egps_jshift, egps_jmask, egps_jstride;
64 int egps_i;
65 real* vvdwtp[UNROLLI];
66 real* vctp[UNROLLI];
67 #endif
68
69 SimdReal shX_S;
70 SimdReal shY_S;
71 SimdReal shZ_S;
72 SimdReal ix_S0, iy_S0, iz_S0;
73 SimdReal ix_S1, iy_S1, iz_S1;
74 SimdReal ix_S2, iy_S2, iz_S2;
75 SimdReal ix_S3, iy_S3, iz_S3;
76 SimdReal fix_S0, fiy_S0, fiz_S0;
77 SimdReal fix_S1, fiy_S1, fiz_S1;
78 SimdReal fix_S2, fiy_S2, fiz_S2;
79 SimdReal fix_S3, fiy_S3, fiz_S3;
80
81 SimdReal diagonal_jmi_S;
82 #if UNROLLI == UNROLLJ
83 SimdBool diagonal_mask_S0, diagonal_mask_S1, diagonal_mask_S2, diagonal_mask_S3;
84 #else
85 SimdBool diagonal_mask0_S0, diagonal_mask0_S1, diagonal_mask0_S2, diagonal_mask0_S3;
86 SimdBool diagonal_mask1_S0, diagonal_mask1_S1, diagonal_mask1_S2, diagonal_mask1_S3;
87 #endif
88
89 SimdBitMask filter_S0, filter_S1, filter_S2, filter_S3;
90
91 SimdReal zero_S(0.0);
92
93 SimdReal one_S(1.0);
94 SimdReal iq_S0 = setZero();
95 SimdReal iq_S1 = setZero();
96 SimdReal iq_S2 = setZero();
97 SimdReal iq_S3 = setZero();
98
99 #ifdef CALC_COUL_RF
100 SimdReal mrc_3_S;
101 # ifdef CALC_ENERGIES
102 SimdReal hrc_3_S, moh_rc_S;
103 # endif
104 #endif
105
106 #ifdef CALC_COUL_TAB
107 /* Coulomb table variables */
108 SimdReal invtsp_S;
109 const real* tab_coul_F;
110 # if defined CALC_ENERGIES && !defined TAB_FDV0
111 const real gmx_unused* tab_coul_V;
112 # endif
113 # ifdef CALC_ENERGIES
114 SimdReal mhalfsp_S;
115 # endif
116 #endif
117
118 #ifdef CALC_COUL_EWALD
119 SimdReal beta2_S, beta_S;
120 #endif
121
122 #if defined CALC_ENERGIES && (defined CALC_COUL_EWALD || defined CALC_COUL_TAB)
123 SimdReal sh_ewald_S;
124 #endif
125
126 #if defined LJ_CUT && defined CALC_ENERGIES
127 SimdReal p6_cpot_S, p12_cpot_S;
128 #endif
129 #ifdef LJ_POT_SWITCH
130 SimdReal rswitch_S;
131 SimdReal swV3_S, swV4_S, swV5_S;
132 SimdReal swF2_S, swF3_S, swF4_S;
133 #endif
134 #ifdef LJ_FORCE_SWITCH
135 SimdReal rswitch_S;
136 SimdReal p6_fc2_S, p6_fc3_S;
137 SimdReal p12_fc2_S, p12_fc3_S;
138 # ifdef CALC_ENERGIES
139 SimdReal p6_vc3_S, p6_vc4_S;
140 SimdReal p12_vc3_S, p12_vc4_S;
141 SimdReal p6_6cpot_S, p12_12cpot_S;
142 # endif
143 #endif
144 #ifdef LJ_EWALD_GEOM
145 real lj_ewaldcoeff2, lj_ewaldcoeff6_6;
146 SimdReal half_S, lje_c2_S, lje_c6_6_S;
147 #endif
148
149 #ifdef LJ_COMB_LB
150 SimdReal hsig_i_S0, seps_i_S0;
151 SimdReal hsig_i_S1, seps_i_S1;
152 SimdReal hsig_i_S2, seps_i_S2;
153 SimdReal hsig_i_S3, seps_i_S3;
154 #endif /* LJ_COMB_LB */
155
156 SimdReal minRsq_S;
157 SimdReal rc2_S;
158 #ifdef VDW_CUTOFF_CHECK
159 SimdReal rcvdw2_S;
160 #endif
161
162 int ninner;
163
164 #ifdef COUNT_PAIRS
165 int npair = 0;
166 #endif
167
168 const nbnxn_atomdata_t::Params& nbatParams = nbat->params();
169
170 #if defined LJ_COMB_GEOM || defined LJ_COMB_LB || defined LJ_EWALD_GEOM
171 const real* gmx_restrict ljc = nbatParams.lj_comb.data();
172 #endif
173 #if !(defined LJ_COMB_GEOM || defined LJ_COMB_LB || defined FIX_LJ_C)
174 /* No combination rule used */
175 const real* gmx_restrict nbfp_ptr = nbatParams.nbfp_aligned.data();
176 const int* gmx_restrict type = nbatParams.type.data();
177 #endif
178
179 /* Load j-i for the first i */
180 diagonal_jmi_S = load<SimdReal>(nbat->simdMasks.diagonal_4xn_j_minus_i.data());
181 /* Generate all the diagonal masks as comparison results */
182 #if UNROLLI == UNROLLJ
183 diagonal_mask_S0 = (zero_S < diagonal_jmi_S);
184 diagonal_jmi_S = diagonal_jmi_S - one_S;
185 diagonal_mask_S1 = (zero_S < diagonal_jmi_S);
186 diagonal_jmi_S = diagonal_jmi_S - one_S;
187 diagonal_mask_S2 = (zero_S < diagonal_jmi_S);
188 diagonal_jmi_S = diagonal_jmi_S - one_S;
189 diagonal_mask_S3 = (zero_S < diagonal_jmi_S);
190 #else
191 # if UNROLLI == 2 * UNROLLJ || 2 * UNROLLI == UNROLLJ
192 diagonal_mask0_S0 = (zero_S < diagonal_jmi_S);
193 diagonal_jmi_S = diagonal_jmi_S - one_S;
194 diagonal_mask0_S1 = (zero_S < diagonal_jmi_S);
195 diagonal_jmi_S = diagonal_jmi_S - one_S;
196 diagonal_mask0_S2 = (zero_S < diagonal_jmi_S);
197 diagonal_jmi_S = diagonal_jmi_S - one_S;
198 diagonal_mask0_S3 = (zero_S < diagonal_jmi_S);
199 diagonal_jmi_S = diagonal_jmi_S - one_S;
200
201 # if UNROLLI == 2 * UNROLLJ
202 /* Load j-i for the second half of the j-cluster */
203 diagonal_jmi_S = load<SimdReal>(nbat->simdMasks.diagonal_4xn_j_minus_i.data() + UNROLLJ);
204 # endif
205
206 diagonal_mask1_S0 = (zero_S < diagonal_jmi_S);
207 diagonal_jmi_S = diagonal_jmi_S - one_S;
208 diagonal_mask1_S1 = (zero_S < diagonal_jmi_S);
209 diagonal_jmi_S = diagonal_jmi_S - one_S;
210 diagonal_mask1_S2 = (zero_S < diagonal_jmi_S);
211 diagonal_jmi_S = diagonal_jmi_S - one_S;
212 diagonal_mask1_S3 = (zero_S < diagonal_jmi_S);
213 # endif
214 #endif
215
216 #if GMX_DOUBLE && !GMX_SIMD_HAVE_INT32_LOGICAL
217 const std::uint64_t* gmx_restrict exclusion_filter = nbat->simdMasks.exclusion_filter64.data();
218 #else
219 const std::uint32_t* gmx_restrict exclusion_filter = nbat->simdMasks.exclusion_filter.data();
220 #endif
221
222 /* Here we cast the exclusion filters from unsigned * to int * or real *.
223 * Since we only check bits, the actual value they represent does not
224 * matter, as long as both filter and mask data are treated the same way.
225 */
226 #if GMX_SIMD_HAVE_INT32_LOGICAL
227 filter_S0 = load<SimdBitMask>(reinterpret_cast<const int*>(exclusion_filter + 0 * UNROLLJ));
228 filter_S1 = load<SimdBitMask>(reinterpret_cast<const int*>(exclusion_filter + 1 * UNROLLJ));
229 filter_S2 = load<SimdBitMask>(reinterpret_cast<const int*>(exclusion_filter + 2 * UNROLLJ));
230 filter_S3 = load<SimdBitMask>(reinterpret_cast<const int*>(exclusion_filter + 3 * UNROLLJ));
231 #else
232 filter_S0 = load<SimdBitMask>(reinterpret_cast<const real*>(exclusion_filter + 0 * UNROLLJ));
233 filter_S1 = load<SimdBitMask>(reinterpret_cast<const real*>(exclusion_filter + 1 * UNROLLJ));
234 filter_S2 = load<SimdBitMask>(reinterpret_cast<const real*>(exclusion_filter + 2 * UNROLLJ));
235 filter_S3 = load<SimdBitMask>(reinterpret_cast<const real*>(exclusion_filter + 3 * UNROLLJ));
236 #endif
237
238 #ifdef CALC_COUL_RF
239 /* Reaction-field constants */
240 mrc_3_S = SimdReal(-2 * ic->k_rf);
241 # ifdef CALC_ENERGIES
242 hrc_3_S = SimdReal(ic->k_rf);
243 moh_rc_S = SimdReal(-ic->c_rf);
244 # endif
245 #endif
246
247 #ifdef CALC_COUL_TAB
248
249 invtsp_S = SimdReal(ic->coulombEwaldTables->scale);
250 # ifdef CALC_ENERGIES
251 mhalfsp_S = SimdReal(-0.5_real / ic->coulombEwaldTables->scale);
252 # endif
253
254 # ifdef TAB_FDV0
255 tab_coul_F = ic->coulombEwaldTables->tableFDV0.data();
256 # else
257 tab_coul_F = ic->coulombEwaldTables->tableF.data();
258 # ifdef CALC_ENERGIES
259 tab_coul_V = ic->coulombEwaldTables->tableV.data();
260 # endif
261 # endif
262 #endif /* CALC_COUL_TAB */
263
264 #ifdef CALC_COUL_EWALD
265 beta2_S = SimdReal(ic->ewaldcoeff_q * ic->ewaldcoeff_q);
266 beta_S = SimdReal(ic->ewaldcoeff_q);
267 #endif
268
269 #if (defined CALC_COUL_TAB || defined CALC_COUL_EWALD) && defined CALC_ENERGIES
270 sh_ewald_S = SimdReal(ic->sh_ewald);
271 #endif
272
273 /* LJ function constants */
274 #if defined CALC_ENERGIES || defined LJ_POT_SWITCH
275 SimdReal sixth_S(1.0 / 6.0);
276 SimdReal twelveth_S(1.0 / 12.0);
277 #endif
278
279 #if defined LJ_CUT && defined CALC_ENERGIES
280 /* We shift the potential by cpot, which can be zero */
281 p6_cpot_S = SimdReal(ic->dispersion_shift.cpot);
282 p12_cpot_S = SimdReal(ic->repulsion_shift.cpot);
283 #endif
284 #ifdef LJ_POT_SWITCH
285 rswitch_S = SimdReal(ic->rvdw_switch);
286 swV3_S = SimdReal(ic->vdw_switch.c3);
287 swV4_S = SimdReal(ic->vdw_switch.c4);
288 swV5_S = SimdReal(ic->vdw_switch.c5);
289 swF2_S = SimdReal(3 * ic->vdw_switch.c3);
290 swF3_S = SimdReal(4 * ic->vdw_switch.c4);
291 swF4_S = SimdReal(5 * ic->vdw_switch.c5);
292 #endif
293 #ifdef LJ_FORCE_SWITCH
294 rswitch_S = SimdReal(ic->rvdw_switch);
295 p6_fc2_S = SimdReal(ic->dispersion_shift.c2);
296 p6_fc3_S = SimdReal(ic->dispersion_shift.c3);
297 p12_fc2_S = SimdReal(ic->repulsion_shift.c2);
298 p12_fc3_S = SimdReal(ic->repulsion_shift.c3);
299 # ifdef CALC_ENERGIES
300 {
301 SimdReal mthird_S(-1.0 / 3.0);
302 SimdReal mfourth_S(-1.0 / 4.0);
303
304 p6_vc3_S = mthird_S * p6_fc2_S;
305 p6_vc4_S = mfourth_S * p6_fc3_S;
306 p6_6cpot_S = SimdReal(ic->dispersion_shift.cpot / 6);
307 p12_vc3_S = mthird_S * p12_fc2_S;
308 p12_vc4_S = mfourth_S * p12_fc3_S;
309 p12_12cpot_S = SimdReal(ic->repulsion_shift.cpot / 12);
310 }
311 # endif
312 #endif
313 #ifdef LJ_EWALD_GEOM
314 half_S = SimdReal(0.5);
315 lj_ewaldcoeff2 = ic->ewaldcoeff_lj * ic->ewaldcoeff_lj;
316 lj_ewaldcoeff6_6 = lj_ewaldcoeff2 * lj_ewaldcoeff2 * lj_ewaldcoeff2 / 6;
317 lje_c2_S = SimdReal(lj_ewaldcoeff2);
318 lje_c6_6_S = SimdReal(lj_ewaldcoeff6_6);
319 # ifdef CALC_ENERGIES
320 /* Determine the grid potential at the cut-off */
321 SimdReal lje_vc_S(ic->sh_lj_ewald);
322 # endif
323 #endif
324
325 /* The kernel either supports rcoulomb = rvdw or rcoulomb >= rvdw */
326 rc2_S = SimdReal(ic->rcoulomb * ic->rcoulomb);
327 #ifdef VDW_CUTOFF_CHECK
328 rcvdw2_S = SimdReal(ic->rvdw * ic->rvdw);
329 #endif
330
331 minRsq_S = SimdReal(c_nbnxnMinDistanceSquared);
332
333 const real* gmx_restrict q = nbatParams.q.data();
334 const real facel = ic->epsfac;
335 const real* gmx_restrict shiftvec = shift_vec[0];
336 const real* gmx_restrict x = nbat->x().data();
337
338 #ifdef FIX_LJ_C
339 alignas(GMX_SIMD_ALIGNMENT) real pvdw_c6[2 * UNROLLI * UNROLLJ];
340 real* pvdw_c12 = pvdw_c6 + UNROLLI * UNROLLJ;
341
342 for (int jp = 0; jp < UNROLLJ; jp++)
343 {
344 pvdw_c6[0 * UNROLLJ + jp] = nbat->nbfp[0 * 2];
345 pvdw_c6[1 * UNROLLJ + jp] = nbat->nbfp[0 * 2];
346 pvdw_c6[2 * UNROLLJ + jp] = nbat->nbfp[0 * 2];
347 pvdw_c6[3 * UNROLLJ + jp] = nbat->nbfp[0 * 2];
348
349 pvdw_c12[0 * UNROLLJ + jp] = nbat->nbfp[0 * 2 + 1];
350 pvdw_c12[1 * UNROLLJ + jp] = nbat->nbfp[0 * 2 + 1];
351 pvdw_c12[2 * UNROLLJ + jp] = nbat->nbfp[0 * 2 + 1];
352 pvdw_c12[3 * UNROLLJ + jp] = nbat->nbfp[0 * 2 + 1];
353 }
354 SimdReal c6_S0 = load<SimdReal>(pvdw_c6 + 0 * UNROLLJ);
355 SimdReal c6_S1 = load<SimdReal>(pvdw_c6 + 1 * UNROLLJ);
356 SimdReal c6_S2 = load<SimdReal>(pvdw_c6 + 2 * UNROLLJ);
357 SimdReal c6_S3 = load<SimdReal>(pvdw_c6 + 3 * UNROLLJ);
358
359 SimdReal c12_S0 = load<SimdReal>(pvdw_c12 + 0 * UNROLLJ);
360 SimdReal c12_S1 = load<SimdReal>(pvdw_c12 + 1 * UNROLLJ);
361 SimdReal c12_S2 = load<SimdReal>(pvdw_c12 + 2 * UNROLLJ);
362 SimdReal c12_S3 = load<SimdReal>(pvdw_c12 + 3 * UNROLLJ);
363 #endif /* FIX_LJ_C */
364
365 #ifdef ENERGY_GROUPS
366 egps_ishift = nbatParams.neg_2log;
367 egps_imask = (1 << egps_ishift) - 1;
368 egps_jshift = 2 * nbatParams.neg_2log;
369 egps_jmask = (1 << egps_jshift) - 1;
370 egps_jstride = (UNROLLJ >> 1) * UNROLLJ;
371 /* Major division is over i-particle energy groups, determine the stride */
372 Vstride_i = nbatParams.nenergrp * (1 << nbatParams.neg_2log) * egps_jstride;
373 #endif
374
375 l_cj = nbl->cj.data();
376
377 ninner = 0;
378
379 for (const nbnxn_ci_t& ciEntry : nbl->ci)
380 {
381 ish = (ciEntry.shift & NBNXN_CI_SHIFT);
382 ish3 = ish * 3;
383 cjind0 = ciEntry.cj_ind_start;
384 cjind1 = ciEntry.cj_ind_end;
385 ci = ciEntry.ci;
386 ci_sh = (ish == CENTRAL ? ci : -1);
387
388 shX_S = SimdReal(shiftvec[ish3]);
389 shY_S = SimdReal(shiftvec[ish3 + 1]);
390 shZ_S = SimdReal(shiftvec[ish3 + 2]);
391
392 #if UNROLLJ <= 4
393 int sci = ci * STRIDE;
394 int scix = sci * DIM;
395 # if defined LJ_COMB_LB || defined LJ_COMB_GEOM || defined LJ_EWALD_GEOM
396 int sci2 = sci * 2;
397 # endif
398 #else
399 int sci = (ci >> 1) * STRIDE;
400 int scix = sci * DIM + (ci & 1) * (STRIDE >> 1);
401 # if defined LJ_COMB_LB || defined LJ_COMB_GEOM || defined LJ_EWALD_GEOM
402 int sci2 = sci * 2 + (ci & 1) * (STRIDE >> 1);
403 # endif
404 sci += (ci & 1) * (STRIDE >> 1);
405 #endif
406
407 /* We have 5 LJ/C combinations, but use only three inner loops,
408 * as the other combinations are unlikely and/or not much faster:
409 * inner half-LJ + C for half-LJ + C / no-LJ + C
410 * inner LJ + C for full-LJ + C
411 * inner LJ for full-LJ + no-C / half-LJ + no-C
412 */
413 do_LJ = ((ciEntry.shift & NBNXN_CI_DO_LJ(0)) != 0);
414 do_coul = ((ciEntry.shift & NBNXN_CI_DO_COUL(0)) != 0);
415 half_LJ = (((ciEntry.shift & NBNXN_CI_HALF_LJ(0)) != 0) || !do_LJ) && do_coul;
416
417 #ifdef ENERGY_GROUPS
418 egps_i = nbatParams.energrp[ci];
419 {
420 int ia, egp_ia;
421
422 for (ia = 0; ia < UNROLLI; ia++)
423 {
424 egp_ia = (egps_i >> (ia * egps_ishift)) & egps_imask;
425 vvdwtp[ia] = Vvdw + egp_ia * Vstride_i;
426 vctp[ia] = Vc + egp_ia * Vstride_i;
427 }
428 }
429 #endif
430
431 #ifdef CALC_ENERGIES
432 # ifdef LJ_EWALD_GEOM
433 gmx_bool do_self = TRUE;
434 # else
435 gmx_bool do_self = do_coul;
436 # endif
437 # if UNROLLJ == 4
438 if (do_self && l_cj[ciEntry.cj_ind_start].cj == ci_sh)
439 # endif
440 # if UNROLLJ == 2
441 if (do_self && l_cj[ciEntry.cj_ind_start].cj == (ci_sh << 1))
442 # endif
443 # if UNROLLJ == 8
444 if (do_self && l_cj[ciEntry.cj_ind_start].cj == (ci_sh >> 1))
445 # endif
446 {
447 if (do_coul)
448 {
449 real Vc_sub_self;
450 int ia;
451
452 # ifdef CALC_COUL_RF
453 Vc_sub_self = 0.5 * ic->c_rf;
454 # endif
455 # ifdef CALC_COUL_TAB
456 # ifdef TAB_FDV0
457 Vc_sub_self = 0.5 * tab_coul_F[2];
458 # else
459 Vc_sub_self = 0.5 * tab_coul_V[0];
460 # endif
461 # endif
462 # ifdef CALC_COUL_EWALD
463 /* beta/sqrt(pi) */
464 Vc_sub_self = 0.5 * ic->ewaldcoeff_q * M_2_SQRTPI;
465 # endif
466
467 for (ia = 0; ia < UNROLLI; ia++)
468 {
469 real qi;
470
471 qi = q[sci + ia];
472 # ifdef ENERGY_GROUPS
473 vctp[ia][((egps_i >> (ia * egps_ishift)) & egps_imask) * egps_jstride]
474 # else
475 Vc[0]
476 # endif
477 -= facel * qi * qi * Vc_sub_self;
478 }
479 }
480
481 # ifdef LJ_EWALD_GEOM
482 {
483 int ia;
484
485 for (ia = 0; ia < UNROLLI; ia++)
486 {
487 real c6_i;
488
489 c6_i = nbatParams.nbfp[nbatParams.type[sci + ia] * (nbatParams.numTypes + 1) * 2]
490 / 6;
491 # ifdef ENERGY_GROUPS
492 vvdwtp[ia][((egps_i >> (ia * egps_ishift)) & egps_imask) * egps_jstride]
493 # else
494 Vvdw[0]
495 # endif
496 += 0.5 * c6_i * lj_ewaldcoeff6_6;
497 }
498 }
499 # endif /* LJ_EWALD_GEOM */
500 }
501 #endif
502
503 /* Load i atom data */
504 int sciy = scix + STRIDE;
505 int sciz = sciy + STRIDE;
506 ix_S0 = SimdReal(x[scix]) + shX_S;
507 ix_S1 = SimdReal(x[scix + 1]) + shX_S;
508 ix_S2 = SimdReal(x[scix + 2]) + shX_S;
509 ix_S3 = SimdReal(x[scix + 3]) + shX_S;
510 iy_S0 = SimdReal(x[sciy]) + shY_S;
511 iy_S1 = SimdReal(x[sciy + 1]) + shY_S;
512 iy_S2 = SimdReal(x[sciy + 2]) + shY_S;
513 iy_S3 = SimdReal(x[sciy + 3]) + shY_S;
514 iz_S0 = SimdReal(x[sciz]) + shZ_S;
515 iz_S1 = SimdReal(x[sciz + 1]) + shZ_S;
516 iz_S2 = SimdReal(x[sciz + 2]) + shZ_S;
517 iz_S3 = SimdReal(x[sciz + 3]) + shZ_S;
518
519 if (do_coul)
520 {
521 iq_S0 = SimdReal(facel * q[sci]);
522 iq_S1 = SimdReal(facel * q[sci + 1]);
523 iq_S2 = SimdReal(facel * q[sci + 2]);
524 iq_S3 = SimdReal(facel * q[sci + 3]);
525 }
526
527 #ifdef LJ_COMB_LB
528 hsig_i_S0 = SimdReal(ljc[sci2 + 0]);
529 hsig_i_S1 = SimdReal(ljc[sci2 + 1]);
530 hsig_i_S2 = SimdReal(ljc[sci2 + 2]);
531 hsig_i_S3 = SimdReal(ljc[sci2 + 3]);
532 seps_i_S0 = SimdReal(ljc[sci2 + STRIDE + 0]);
533 seps_i_S1 = SimdReal(ljc[sci2 + STRIDE + 1]);
534 seps_i_S2 = SimdReal(ljc[sci2 + STRIDE + 2]);
535 seps_i_S3 = SimdReal(ljc[sci2 + STRIDE + 3]);
536 #else
537 # ifdef LJ_COMB_GEOM
538 SimdReal c6s_S0 = SimdReal(ljc[sci2 + 0]);
539 SimdReal c6s_S1 = SimdReal(ljc[sci2 + 1]);
540 SimdReal c6s_S2, c6s_S3;
541 if (!half_LJ)
542 {
543 c6s_S2 = SimdReal(ljc[sci2 + 2]);
544 c6s_S3 = SimdReal(ljc[sci2 + 3]);
545 }
546 else
547 {
548 c6s_S2 = setZero();
549 c6s_S3 = setZero();
550 }
551 SimdReal c12s_S0 = SimdReal(ljc[sci2 + STRIDE + 0]);
552 SimdReal c12s_S1 = SimdReal(ljc[sci2 + STRIDE + 1]);
553 SimdReal c12s_S2, c12s_S3;
554 if (!half_LJ)
555 {
556 c12s_S2 = SimdReal(ljc[sci2 + STRIDE + 2]);
557 c12s_S3 = SimdReal(ljc[sci2 + STRIDE + 3]);
558 }
559 else
560 {
561 c12s_S2 = setZero();
562 c12s_S3 = setZero();
563 }
564 # else
565 const int numTypes = nbatParams.numTypes;
566 const real* nbfp0 = nbfp_ptr + type[sci] * numTypes * c_simdBestPairAlignment;
567 const real* nbfp1 = nbfp_ptr + type[sci + 1] * numTypes * c_simdBestPairAlignment;
568 const real *nbfp2 = nullptr, *nbfp3 = nullptr;
569 if (!half_LJ)
570 {
571 nbfp2 = nbfp_ptr + type[sci + 2] * numTypes * c_simdBestPairAlignment;
572 nbfp3 = nbfp_ptr + type[sci + 3] * numTypes * c_simdBestPairAlignment;
573 }
574 # endif
575 #endif
576 #ifdef LJ_EWALD_GEOM
577 /* We need the geometrically combined C6 for the PME grid correction */
578 SimdReal c6s_S0 = SimdReal(ljc[sci2 + 0]);
579 SimdReal c6s_S1 = SimdReal(ljc[sci2 + 1]);
580 SimdReal c6s_S2, c6s_S3;
581 if (!half_LJ)
582 {
583 c6s_S2 = SimdReal(ljc[sci2 + 2]);
584 c6s_S3 = SimdReal(ljc[sci2 + 3]);
585 }
586 #endif
587
588 /* Zero the potential energy for this list */
589 #ifdef CALC_ENERGIES
590 SimdReal Vvdwtot_S = setZero();
591 SimdReal vctot_S = setZero();
592 #endif
593
594 /* Clear i atom forces */
595 fix_S0 = setZero();
596 fix_S1 = setZero();
597 fix_S2 = setZero();
598 fix_S3 = setZero();
599 fiy_S0 = setZero();
600 fiy_S1 = setZero();
601 fiy_S2 = setZero();
602 fiy_S3 = setZero();
603 fiz_S0 = setZero();
604 fiz_S1 = setZero();
605 fiz_S2 = setZero();
606 fiz_S3 = setZero();
607
608 cjind = cjind0;
609
610 /* Currently all kernels use (at least half) LJ */
611 #define CALC_LJ
612 if (half_LJ)
613 {
614 /* Coulomb: all i-atoms, LJ: first half i-atoms */
615 #define CALC_COULOMB
616 #define HALF_LJ
617 #define CHECK_EXCLS
618 while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
619 {
620 #include "kernel_inner.h"
621 cjind++;
622 }
623 #undef CHECK_EXCLS
624 for (; (cjind < cjind1); cjind++)
625 {
626 #include "kernel_inner.h"
627 }
628 #undef HALF_LJ
629 #undef CALC_COULOMB
630 }
631 else if (do_coul)
632 {
633 /* Coulomb: all i-atoms, LJ: all i-atoms */
634 #define CALC_COULOMB
635 #define CHECK_EXCLS
636 while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
637 {
638 #include "kernel_inner.h"
639 cjind++;
640 }
641 #undef CHECK_EXCLS
642 for (; (cjind < cjind1); cjind++)
643 {
644 #include "kernel_inner.h"
645 }
646 #undef CALC_COULOMB
647 }
648 else
649 {
650 /* Coulomb: none, LJ: all i-atoms */
651 #define CHECK_EXCLS
652 while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
653 {
654 #include "kernel_inner.h"
655 cjind++;
656 }
657 #undef CHECK_EXCLS
658 for (; (cjind < cjind1); cjind++)
659 {
660 #include "kernel_inner.h"
661 }
662 }
663 #undef CALC_LJ
664 ninner += cjind1 - cjind0;
665
666 /* Add accumulated i-forces to the force array */
667 real fShiftX = reduceIncr4ReturnSum(f + scix, fix_S0, fix_S1, fix_S2, fix_S3);
668 real fShiftY = reduceIncr4ReturnSum(f + sciy, fiy_S0, fiy_S1, fiy_S2, fiy_S3);
669 real fShiftZ = reduceIncr4ReturnSum(f + sciz, fiz_S0, fiz_S1, fiz_S2, fiz_S3);
670
671
672 #ifdef CALC_SHIFTFORCES
673 fshift[ish3 + 0] += fShiftX;
674 fshift[ish3 + 1] += fShiftY;
675 fshift[ish3 + 2] += fShiftZ;
676 #endif
677
678 #ifdef CALC_ENERGIES
679 if (do_coul)
680 {
681 *Vc += reduce(vctot_S);
682 }
683
684 *Vvdw += reduce(Vvdwtot_S);
685 #endif
686
687 /* Outer loop uses 6 flops/iteration */
688 }
689
690 #ifdef COUNT_PAIRS
691 printf("atom pairs %d\n", npair);
692 #endif
693 }
The ultimate molecular dynamics simulation package