2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
39 const nbnxn_cj_t *l_cj;
42 gmx_bool do_LJ, half_LJ, do_coul;
43 int cjind0, cjind1, cjind;
47 int egps_ishift, egps_imask;
48 int egps_jshift, egps_jmask, egps_jstride;
50 real *vvdwtp[UNROLLI];
57 SimdReal ix_S0, iy_S0, iz_S0;
58 SimdReal ix_S1, iy_S1, iz_S1;
59 SimdReal ix_S2, iy_S2, iz_S2;
60 SimdReal ix_S3, iy_S3, iz_S3;
61 SimdReal fix_S0, fiy_S0, fiz_S0;
62 SimdReal fix_S1, fiy_S1, fiz_S1;
63 SimdReal fix_S2, fiy_S2, fiz_S2;
64 SimdReal fix_S3, fiy_S3, fiz_S3;
66 SimdReal diagonal_jmi_S;
67 #if UNROLLI == UNROLLJ
68 SimdBool diagonal_mask_S0, diagonal_mask_S1, diagonal_mask_S2, diagonal_mask_S3;
70 SimdBool diagonal_mask0_S0, diagonal_mask0_S1, diagonal_mask0_S2, diagonal_mask0_S3;
71 SimdBool diagonal_mask1_S0, diagonal_mask1_S1, diagonal_mask1_S2, diagonal_mask1_S3;
74 SimdBitMask filter_S0, filter_S1, filter_S2, filter_S3;
79 SimdReal iq_S0 = setZero();
80 SimdReal iq_S1 = setZero();
81 SimdReal iq_S2 = setZero();
82 SimdReal iq_S3 = setZero();
87 SimdReal hrc_3_S, moh_rc_S;
92 /* Coulomb table variables */
94 const real * tab_coul_F;
96 const real gmx_unused *tab_coul_V;
103 #ifdef CALC_COUL_EWALD
104 SimdReal beta2_S, beta_S;
107 #if defined CALC_ENERGIES && (defined CALC_COUL_EWALD || defined CALC_COUL_TAB)
111 #if defined LJ_CUT && defined CALC_ENERGIES
112 SimdReal p6_cpot_S, p12_cpot_S;
116 SimdReal swV3_S, swV4_S, swV5_S;
117 SimdReal swF2_S, swF3_S, swF4_S;
119 #ifdef LJ_FORCE_SWITCH
121 SimdReal p6_fc2_S, p6_fc3_S;
122 SimdReal p12_fc2_S, p12_fc3_S;
124 SimdReal p6_vc3_S, p6_vc4_S;
125 SimdReal p12_vc3_S, p12_vc4_S;
126 SimdReal p6_6cpot_S, p12_12cpot_S;
130 real lj_ewaldcoeff2, lj_ewaldcoeff6_6;
131 SimdReal mone_S, half_S, lje_c2_S, lje_c6_6_S;
135 SimdReal hsig_i_S0, seps_i_S0;
136 SimdReal hsig_i_S1, seps_i_S1;
137 SimdReal hsig_i_S2, seps_i_S2;
138 SimdReal hsig_i_S3, seps_i_S3;
139 #endif /* LJ_COMB_LB */
143 #ifdef VDW_CUTOFF_CHECK
153 const nbnxn_atomdata_t::Params &nbatParams = nbat->params();
155 #if defined LJ_COMB_GEOM || defined LJ_COMB_LB || defined LJ_EWALD_GEOM
156 const real * gmx_restrict ljc = nbatParams.lj_comb.data();
158 #if !(defined LJ_COMB_GEOM || defined LJ_COMB_LB || defined FIX_LJ_C)
159 /* No combination rule used */
160 const real * gmx_restrict nbfp_ptr = nbatParams.nbfp_aligned.data();
161 const int * gmx_restrict type = nbatParams.type.data();
164 /* Load j-i for the first i */
165 diagonal_jmi_S = load<SimdReal>(nbat->simdMasks.diagonal_4xn_j_minus_i.data());
166 /* Generate all the diagonal masks as comparison results */
167 #if UNROLLI == UNROLLJ
168 diagonal_mask_S0 = (zero_S < diagonal_jmi_S);
169 diagonal_jmi_S = diagonal_jmi_S - one_S;
170 diagonal_mask_S1 = (zero_S < diagonal_jmi_S);
171 diagonal_jmi_S = diagonal_jmi_S - one_S;
172 diagonal_mask_S2 = (zero_S < diagonal_jmi_S);
173 diagonal_jmi_S = diagonal_jmi_S - one_S;
174 diagonal_mask_S3 = (zero_S < diagonal_jmi_S);
176 #if UNROLLI == 2*UNROLLJ || 2*UNROLLI == UNROLLJ
177 diagonal_mask0_S0 = (zero_S < diagonal_jmi_S);
178 diagonal_jmi_S = diagonal_jmi_S - one_S;
179 diagonal_mask0_S1 = (zero_S < diagonal_jmi_S);
180 diagonal_jmi_S = diagonal_jmi_S - one_S;
181 diagonal_mask0_S2 = (zero_S < diagonal_jmi_S);
182 diagonal_jmi_S = diagonal_jmi_S - one_S;
183 diagonal_mask0_S3 = (zero_S < diagonal_jmi_S);
184 diagonal_jmi_S = diagonal_jmi_S - one_S;
186 #if UNROLLI == 2*UNROLLJ
187 /* Load j-i for the second half of the j-cluster */
188 diagonal_jmi_S = load<SimdReal>(nbat->simdMasks.diagonal_4xn_j_minus_i.data() + UNROLLJ);
191 diagonal_mask1_S0 = (zero_S < diagonal_jmi_S);
192 diagonal_jmi_S = diagonal_jmi_S - one_S;
193 diagonal_mask1_S1 = (zero_S < diagonal_jmi_S);
194 diagonal_jmi_S = diagonal_jmi_S - one_S;
195 diagonal_mask1_S2 = (zero_S < diagonal_jmi_S);
196 diagonal_jmi_S = diagonal_jmi_S - one_S;
197 diagonal_mask1_S3 = (zero_S < diagonal_jmi_S);
201 #if GMX_DOUBLE && !GMX_SIMD_HAVE_INT32_LOGICAL
202 const std::uint64_t * gmx_restrict exclusion_filter = nbat->simdMasks.exclusion_filter64.data();
204 const std::uint32_t * gmx_restrict exclusion_filter = nbat->simdMasks.exclusion_filter.data();
207 /* Here we cast the exclusion filters from unsigned * to int * or real *.
208 * Since we only check bits, the actual value they represent does not
209 * matter, as long as both filter and mask data are treated the same way.
211 #if GMX_SIMD_HAVE_INT32_LOGICAL
212 filter_S0 = load<SimdBitMask>(reinterpret_cast<const int *>(exclusion_filter + 0*UNROLLJ));
213 filter_S1 = load<SimdBitMask>(reinterpret_cast<const int *>(exclusion_filter + 1*UNROLLJ));
214 filter_S2 = load<SimdBitMask>(reinterpret_cast<const int *>(exclusion_filter + 2*UNROLLJ));
215 filter_S3 = load<SimdBitMask>(reinterpret_cast<const int *>(exclusion_filter + 3*UNROLLJ));
217 filter_S0 = load<SimdBitMask>(reinterpret_cast<const real *>(exclusion_filter + 0*UNROLLJ));
218 filter_S1 = load<SimdBitMask>(reinterpret_cast<const real *>(exclusion_filter + 1*UNROLLJ));
219 filter_S2 = load<SimdBitMask>(reinterpret_cast<const real *>(exclusion_filter + 2*UNROLLJ));
220 filter_S3 = load<SimdBitMask>(reinterpret_cast<const real *>(exclusion_filter + 3*UNROLLJ));
224 /* Reaction-field constants */
225 mrc_3_S = SimdReal(-2*ic->k_rf);
227 hrc_3_S = SimdReal(ic->k_rf);
228 moh_rc_S = SimdReal(-ic->c_rf);
234 invtsp_S = SimdReal(ic->tabq_scale);
236 mhalfsp_S = SimdReal(-0.5/ic->tabq_scale);
240 tab_coul_F = ic->tabq_coul_FDV0;
242 tab_coul_F = ic->tabq_coul_F;
243 tab_coul_V = ic->tabq_coul_V;
245 #endif /* CALC_COUL_TAB */
247 #ifdef CALC_COUL_EWALD
248 beta2_S = SimdReal(ic->ewaldcoeff_q*ic->ewaldcoeff_q);
249 beta_S = SimdReal(ic->ewaldcoeff_q);
252 #if (defined CALC_COUL_TAB || defined CALC_COUL_EWALD) && defined CALC_ENERGIES
253 sh_ewald_S = SimdReal(ic->sh_ewald);
256 /* LJ function constants */
257 #if defined CALC_ENERGIES || defined LJ_POT_SWITCH
258 SimdReal sixth_S(1.0/6.0);
259 SimdReal twelveth_S(1.0/12.0);
262 #if defined LJ_CUT && defined CALC_ENERGIES
263 /* We shift the potential by cpot, which can be zero */
264 p6_cpot_S = SimdReal(ic->dispersion_shift.cpot);
265 p12_cpot_S = SimdReal(ic->repulsion_shift.cpot);
268 rswitch_S = SimdReal(ic->rvdw_switch);
269 swV3_S = SimdReal(ic->vdw_switch.c3);
270 swV4_S = SimdReal(ic->vdw_switch.c4);
271 swV5_S = SimdReal(ic->vdw_switch.c5);
272 swF2_S = SimdReal(3*ic->vdw_switch.c3);
273 swF3_S = SimdReal(4*ic->vdw_switch.c4);
274 swF4_S = SimdReal(5*ic->vdw_switch.c5);
276 #ifdef LJ_FORCE_SWITCH
277 rswitch_S = SimdReal(ic->rvdw_switch);
278 p6_fc2_S = SimdReal(ic->dispersion_shift.c2);
279 p6_fc3_S = SimdReal(ic->dispersion_shift.c3);
280 p12_fc2_S = SimdReal(ic->repulsion_shift.c2);
281 p12_fc3_S = SimdReal(ic->repulsion_shift.c3);
284 SimdReal mthird_S(-1.0/3.0);
285 SimdReal mfourth_S(-1.0/4.0);
287 p6_vc3_S = mthird_S * p6_fc2_S;
288 p6_vc4_S = mfourth_S * p6_fc3_S;
289 p6_6cpot_S = SimdReal(ic->dispersion_shift.cpot/6);
290 p12_vc3_S = mthird_S * p12_fc2_S;
291 p12_vc4_S = mfourth_S * p12_fc3_S;
292 p12_12cpot_S = SimdReal(ic->repulsion_shift.cpot/12);
297 mone_S = SimdReal(-1.0);
298 half_S = SimdReal(0.5);
299 lj_ewaldcoeff2 = ic->ewaldcoeff_lj*ic->ewaldcoeff_lj;
300 lj_ewaldcoeff6_6 = lj_ewaldcoeff2*lj_ewaldcoeff2*lj_ewaldcoeff2/6;
301 lje_c2_S = SimdReal(lj_ewaldcoeff2);
302 lje_c6_6_S = SimdReal(lj_ewaldcoeff6_6);
304 /* Determine the grid potential at the cut-off */
305 SimdReal lje_vc_S(ic->sh_lj_ewald);
309 /* The kernel either supports rcoulomb = rvdw or rcoulomb >= rvdw */
310 rc2_S = SimdReal(ic->rcoulomb*ic->rcoulomb);
311 #ifdef VDW_CUTOFF_CHECK
312 rcvdw2_S = SimdReal(ic->rvdw*ic->rvdw);
315 minRsq_S = SimdReal(NBNXN_MIN_RSQ);
317 const real * gmx_restrict q = nbatParams.q.data();
318 const real facel = ic->epsfac;
319 const real * gmx_restrict shiftvec = shift_vec[0];
320 const real * gmx_restrict x = nbat->x().data();
323 alignas(GMX_SIMD_ALIGNMENT) real pvdw_c6[2*UNROLLI*UNROLLJ];
324 real *pvdw_c12 = pvdw_c6 + UNROLLI*UNROLLJ;
326 for (int jp = 0; jp < UNROLLJ; jp++)
328 pvdw_c6 [0*UNROLLJ+jp] = nbat->nbfp[0*2];
329 pvdw_c6 [1*UNROLLJ+jp] = nbat->nbfp[0*2];
330 pvdw_c6 [2*UNROLLJ+jp] = nbat->nbfp[0*2];
331 pvdw_c6 [3*UNROLLJ+jp] = nbat->nbfp[0*2];
333 pvdw_c12[0*UNROLLJ+jp] = nbat->nbfp[0*2+1];
334 pvdw_c12[1*UNROLLJ+jp] = nbat->nbfp[0*2+1];
335 pvdw_c12[2*UNROLLJ+jp] = nbat->nbfp[0*2+1];
336 pvdw_c12[3*UNROLLJ+jp] = nbat->nbfp[0*2+1];
338 SimdReal c6_S0 = load<SimdReal>(pvdw_c6 +0*UNROLLJ);
339 SimdReal c6_S1 = load<SimdReal>(pvdw_c6 +1*UNROLLJ);
340 SimdReal c6_S2 = load<SimdReal>(pvdw_c6 +2*UNROLLJ);
341 SimdReal c6_S3 = load<SimdReal>(pvdw_c6 +3*UNROLLJ);
343 SimdReal c12_S0 = load<SimdReal>(pvdw_c12+0*UNROLLJ);
344 SimdReal c12_S1 = load<SimdReal>(pvdw_c12+1*UNROLLJ);
345 SimdReal c12_S2 = load<SimdReal>(pvdw_c12+2*UNROLLJ);
346 SimdReal c12_S3 = load<SimdReal>(pvdw_c12+3*UNROLLJ);
347 #endif /* FIX_LJ_C */
350 egps_ishift = nbatParams.neg_2log;
351 egps_imask = (1<<egps_ishift) - 1;
352 egps_jshift = 2*nbatParams.neg_2log;
353 egps_jmask = (1<<egps_jshift) - 1;
354 egps_jstride = (UNROLLJ>>1)*UNROLLJ;
355 /* Major division is over i-particle energy groups, determine the stride */
356 Vstride_i = nbatParams.nenergrp*(1 << nbatParams.neg_2log)*egps_jstride;
359 l_cj = nbl->cj.data();
363 for (const nbnxn_ci_t &ciEntry : nbl->ci)
365 ish = (ciEntry.shift & NBNXN_CI_SHIFT);
367 cjind0 = ciEntry.cj_ind_start;
368 cjind1 = ciEntry.cj_ind_end;
370 ci_sh = (ish == CENTRAL ? ci : -1);
372 shX_S = SimdReal(shiftvec[ish3]);
373 shY_S = SimdReal(shiftvec[ish3+1]);
374 shZ_S = SimdReal(shiftvec[ish3+2]);
379 #if defined LJ_COMB_LB || defined LJ_COMB_GEOM || defined LJ_EWALD_GEOM
383 int sci = (ci>>1)*STRIDE;
384 int scix = sci*DIM + (ci & 1)*(STRIDE>>1);
385 #if defined LJ_COMB_LB || defined LJ_COMB_GEOM || defined LJ_EWALD_GEOM
386 int sci2 = sci*2 + (ci & 1)*(STRIDE>>1);
388 sci += (ci & 1)*(STRIDE>>1);
391 /* We have 5 LJ/C combinations, but use only three inner loops,
392 * as the other combinations are unlikely and/or not much faster:
393 * inner half-LJ + C for half-LJ + C / no-LJ + C
394 * inner LJ + C for full-LJ + C
395 * inner LJ for full-LJ + no-C / half-LJ + no-C
397 do_LJ = ((ciEntry.shift & NBNXN_CI_DO_LJ(0)) != 0);
398 do_coul = ((ciEntry.shift & NBNXN_CI_DO_COUL(0)) != 0);
399 half_LJ = (((ciEntry.shift & NBNXN_CI_HALF_LJ(0)) != 0) || !do_LJ) && do_coul;
402 egps_i = nbatParams.energrp[ci];
406 for (ia = 0; ia < UNROLLI; ia++)
408 egp_ia = (egps_i >> (ia*egps_ishift)) & egps_imask;
409 vvdwtp[ia] = Vvdw + egp_ia*Vstride_i;
410 vctp[ia] = Vc + egp_ia*Vstride_i;
417 gmx_bool do_self = TRUE;
419 gmx_bool do_self = do_coul;
422 if (do_self && l_cj[ciEntry.cj_ind_start].cj == ci_sh)
425 if (do_self && l_cj[ciEntry.cj_ind_start].cj == (ci_sh<<1))
428 if (do_self && l_cj[ciEntry.cj_ind_start].cj == (ci_sh>>1))
437 Vc_sub_self = 0.5*ic->c_rf;
441 Vc_sub_self = 0.5*tab_coul_F[2];
443 Vc_sub_self = 0.5*tab_coul_V[0];
446 #ifdef CALC_COUL_EWALD
448 Vc_sub_self = 0.5*ic->ewaldcoeff_q*M_2_SQRTPI;
451 for (ia = 0; ia < UNROLLI; ia++)
457 vctp[ia][((egps_i>>(ia*egps_ishift)) & egps_imask)*egps_jstride]
461 -= facel*qi*qi*Vc_sub_self;
469 for (ia = 0; ia < UNROLLI; ia++)
473 c6_i = nbatParams.nbfp[nbatParams.type[sci+ia]*(nbatParams.numTypes + 1)*2]/6;
475 vvdwtp[ia][((egps_i>>(ia*egps_ishift)) & egps_imask)*egps_jstride]
479 += 0.5*c6_i*lj_ewaldcoeff6_6;
482 #endif /* LJ_EWALD_GEOM */
486 /* Load i atom data */
487 int sciy = scix + STRIDE;
488 int sciz = sciy + STRIDE;
489 ix_S0 = SimdReal(x[scix]) + shX_S;
490 ix_S1 = SimdReal(x[scix+1]) + shX_S;
491 ix_S2 = SimdReal(x[scix+2]) + shX_S;
492 ix_S3 = SimdReal(x[scix+3]) + shX_S;
493 iy_S0 = SimdReal(x[sciy]) + shY_S;
494 iy_S1 = SimdReal(x[sciy+1]) + shY_S;
495 iy_S2 = SimdReal(x[sciy+2]) + shY_S;
496 iy_S3 = SimdReal(x[sciy+3]) + shY_S;
497 iz_S0 = SimdReal(x[sciz]) + shZ_S;
498 iz_S1 = SimdReal(x[sciz+1]) + shZ_S;
499 iz_S2 = SimdReal(x[sciz+2]) + shZ_S;
500 iz_S3 = SimdReal(x[sciz+3]) + shZ_S;
504 iq_S0 = SimdReal(facel*q[sci]);
505 iq_S1 = SimdReal(facel*q[sci+1]);
506 iq_S2 = SimdReal(facel*q[sci+2]);
507 iq_S3 = SimdReal(facel*q[sci+3]);
511 hsig_i_S0 = SimdReal(ljc[sci2+0]);
512 hsig_i_S1 = SimdReal(ljc[sci2+1]);
513 hsig_i_S2 = SimdReal(ljc[sci2+2]);
514 hsig_i_S3 = SimdReal(ljc[sci2+3]);
515 seps_i_S0 = SimdReal(ljc[sci2+STRIDE+0]);
516 seps_i_S1 = SimdReal(ljc[sci2+STRIDE+1]);
517 seps_i_S2 = SimdReal(ljc[sci2+STRIDE+2]);
518 seps_i_S3 = SimdReal(ljc[sci2+STRIDE+3]);
521 SimdReal c6s_S0 = SimdReal(ljc[sci2+0]);
522 SimdReal c6s_S1 = SimdReal(ljc[sci2+1]);
523 SimdReal c6s_S2, c6s_S3;
526 c6s_S2 = SimdReal(ljc[sci2+2]);
527 c6s_S3 = SimdReal(ljc[sci2+3]);
534 SimdReal c12s_S0 = SimdReal(ljc[sci2+STRIDE+0]);
535 SimdReal c12s_S1 = SimdReal(ljc[sci2+STRIDE+1]);
536 SimdReal c12s_S2, c12s_S3;
539 c12s_S2 = SimdReal(ljc[sci2+STRIDE+2]);
540 c12s_S3 = SimdReal(ljc[sci2+STRIDE+3]);
548 const int numTypes = nbatParams.numTypes;
549 const real *nbfp0 = nbfp_ptr + type[sci ]*numTypes*c_simdBestPairAlignment;
550 const real *nbfp1 = nbfp_ptr + type[sci+1]*numTypes*c_simdBestPairAlignment;
551 const real *nbfp2 = nullptr, *nbfp3 = nullptr;
554 nbfp2 = nbfp_ptr + type[sci+2]*numTypes*c_simdBestPairAlignment;
555 nbfp3 = nbfp_ptr + type[sci+3]*numTypes*c_simdBestPairAlignment;
560 /* We need the geometrically combined C6 for the PME grid correction */
561 SimdReal c6s_S0 = SimdReal(ljc[sci2+0]);
562 SimdReal c6s_S1 = SimdReal(ljc[sci2+1]);
563 SimdReal c6s_S2, c6s_S3;
566 c6s_S2 = SimdReal(ljc[sci2+2]);
567 c6s_S3 = SimdReal(ljc[sci2+3]);
571 /* Zero the potential energy for this list */
573 SimdReal Vvdwtot_S = setZero();
574 SimdReal vctot_S = setZero();
577 /* Clear i atom forces */
593 /* Currently all kernels use (at least half) LJ */
597 /* Coulomb: all i-atoms, LJ: first half i-atoms */
601 while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
603 #include "gromacs/nbnxm/kernels_simd_4xm/kernel_inner.h"
607 for (; (cjind < cjind1); cjind++)
609 #include "gromacs/nbnxm/kernels_simd_4xm/kernel_inner.h"
616 /* Coulomb: all i-atoms, LJ: all i-atoms */
619 while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
621 #include "gromacs/nbnxm/kernels_simd_4xm/kernel_inner.h"
625 for (; (cjind < cjind1); cjind++)
627 #include "gromacs/nbnxm/kernels_simd_4xm/kernel_inner.h"
633 /* Coulomb: none, LJ: all i-atoms */
635 while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
637 #include "gromacs/nbnxm/kernels_simd_4xm/kernel_inner.h"
641 for (; (cjind < cjind1); cjind++)
643 #include "gromacs/nbnxm/kernels_simd_4xm/kernel_inner.h"
647 ninner += cjind1 - cjind0;
649 /* Add accumulated i-forces to the force array */
650 real fShiftX = reduceIncr4ReturnSum(f+scix, fix_S0, fix_S1, fix_S2, fix_S3);
651 real fShiftY = reduceIncr4ReturnSum(f+sciy, fiy_S0, fiy_S1, fiy_S2, fiy_S3);
652 real fShiftZ = reduceIncr4ReturnSum(f+sciz, fiz_S0, fiz_S1, fiz_S2, fiz_S3);
655 #ifdef CALC_SHIFTFORCES
656 fshift[ish3+0] += fShiftX;
657 fshift[ish3+1] += fShiftY;
658 fshift[ish3+2] += fShiftZ;
664 *Vc += reduce(vctot_S);
667 *Vvdw += reduce(Vvdwtot_S);
670 /* Outer loop uses 6 flops/iteration */
674 printf("atom pairs %d\n", npair);