2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
41 /* Unpack pointers for output */
42 real
* f
= out
->f
.data();
43 real
* fshift
= out
->fshift
.data();
46 real
* Vvdw
= out
->VSvdw
.data();
47 real
* Vc
= out
->VSc
.data();
49 real
* Vvdw
= out
->Vvdw
.data();
50 real
* Vc
= out
->Vc
.data();
54 const nbnxn_cj_t
* l_cj
;
57 gmx_bool do_LJ
, half_LJ
, do_coul
;
58 int cjind0
, cjind1
, cjind
;
62 int egps_ishift
, egps_imask
;
63 int egps_jshift
, egps_jmask
, egps_jstride
;
65 real
* vvdwtp
[UNROLLI
];
72 SimdReal ix_S0
, iy_S0
, iz_S0
;
73 SimdReal ix_S2
, iy_S2
, iz_S2
;
74 SimdReal fix_S0
, fiy_S0
, fiz_S0
;
75 SimdReal fix_S2
, fiy_S2
, fiz_S2
;
77 SimdReal diagonal_jmi_S
;
78 #if UNROLLI == UNROLLJ
79 SimdBool diagonal_mask_S0
, diagonal_mask_S2
;
81 SimdBool diagonal_mask0_S0
, diagonal_mask0_S2
;
82 SimdBool diagonal_mask1_S0
, diagonal_mask1_S2
;
85 SimdBitMask filter_S0
, filter_S2
;
90 SimdReal iq_S0
= setZero();
91 SimdReal iq_S2
= setZero();
96 SimdReal hrc_3_S
, moh_rc_S
;
101 /* Coulomb table variables */
103 const real
* tab_coul_F
;
104 # if defined CALC_ENERGIES && !defined TAB_FDV0
105 const real
* tab_coul_V
;
108 # ifdef CALC_ENERGIES
113 #ifdef CALC_COUL_EWALD
114 SimdReal beta2_S
, beta_S
;
117 #if defined CALC_ENERGIES && (defined CALC_COUL_EWALD || defined CALC_COUL_TAB)
121 #if defined LJ_CUT && defined CALC_ENERGIES
122 SimdReal p6_cpot_S
, p12_cpot_S
;
126 SimdReal swV3_S
, swV4_S
, swV5_S
;
127 SimdReal swF2_S
, swF3_S
, swF4_S
;
129 #ifdef LJ_FORCE_SWITCH
131 SimdReal p6_fc2_S
, p6_fc3_S
;
132 SimdReal p12_fc2_S
, p12_fc3_S
;
133 # ifdef CALC_ENERGIES
134 SimdReal p6_vc3_S
, p6_vc4_S
;
135 SimdReal p12_vc3_S
, p12_vc4_S
;
136 SimdReal p6_6cpot_S
, p12_12cpot_S
;
140 real lj_ewaldcoeff2
, lj_ewaldcoeff6_6
;
141 SimdReal mone_S
, half_S
, lje_c2_S
, lje_c6_6_S
;
145 SimdReal hsig_i_S0
, seps_i_S0
;
146 SimdReal hsig_i_S2
, seps_i_S2
;
149 alignas(GMX_SIMD_ALIGNMENT
) real pvdw_c6
[2 * UNROLLI
* UNROLLJ
];
150 real
* pvdw_c12
= pvdw_c6
+ UNROLLI
* UNROLLJ
;
152 #endif /* LJ_COMB_LB */
156 #ifdef VDW_CUTOFF_CHECK
166 const nbnxn_atomdata_t::Params
& nbatParams
= nbat
->params();
168 #if defined LJ_COMB_GEOM || defined LJ_COMB_LB || defined LJ_EWALD_GEOM
169 const real
* gmx_restrict ljc
= nbatParams
.lj_comb
.data();
171 #if !(defined LJ_COMB_GEOM || defined LJ_COMB_LB || defined FIX_LJ_C)
172 /* No combination rule used */
173 const real
* gmx_restrict nbfp_ptr
= nbatParams
.nbfp_aligned
.data();
174 const int* gmx_restrict type
= nbatParams
.type
.data();
177 /* Load j-i for the first i */
178 diagonal_jmi_S
= load
<SimdReal
>(nbat
->simdMasks
.diagonal_2xnn_j_minus_i
.data());
179 /* Generate all the diagonal masks as comparison results */
180 #if UNROLLI == UNROLLJ
181 diagonal_mask_S0
= (zero_S
< diagonal_jmi_S
);
182 diagonal_jmi_S
= diagonal_jmi_S
- one_S
;
183 diagonal_jmi_S
= diagonal_jmi_S
- one_S
;
184 diagonal_mask_S2
= (zero_S
< diagonal_jmi_S
);
186 # if 2 * UNROLLI == UNROLLJ
187 diagonal_mask0_S0
= (zero_S
< diagonal_jmi_S
);
188 diagonal_jmi_S
= diagonal_jmi_S
- one_S
;
189 diagonal_jmi_S
= diagonal_jmi_S
- one_S
;
190 diagonal_mask0_S2
= (zero_S
< diagonal_jmi_S
);
191 diagonal_jmi_S
= diagonal_jmi_S
- one_S
;
192 diagonal_jmi_S
= diagonal_jmi_S
- one_S
;
193 diagonal_mask1_S0
= (zero_S
< diagonal_jmi_S
);
194 diagonal_jmi_S
= diagonal_jmi_S
- one_S
;
195 diagonal_jmi_S
= diagonal_jmi_S
- one_S
;
196 diagonal_mask1_S2
= (zero_S
< diagonal_jmi_S
);
200 /* Load masks for topology exclusion masking. filter_stride is
201 static const, so the conditional will be optimized away. */
202 #if GMX_DOUBLE && !GMX_SIMD_HAVE_INT32_LOGICAL
203 const std::uint64_t* gmx_restrict exclusion_filter
= nbat
->simdMasks
.exclusion_filter64
.data();
205 const std::uint32_t* gmx_restrict exclusion_filter
= nbat
->simdMasks
.exclusion_filter
.data();
208 /* Here we cast the exclusion filters from unsigned * to int * or real *.
209 * Since we only check bits, the actual value they represent does not
210 * matter, as long as both filter and mask data are treated the same way.
212 #if GMX_SIMD_HAVE_INT32_LOGICAL
213 filter_S0
= load
<SimdBitMask
>(reinterpret_cast<const int*>(exclusion_filter
+ 0 * UNROLLJ
));
214 filter_S2
= load
<SimdBitMask
>(reinterpret_cast<const int*>(exclusion_filter
+ 2 * UNROLLJ
));
216 filter_S0
= load
<SimdBitMask
>(reinterpret_cast<const real
*>(exclusion_filter
+ 0 * UNROLLJ
));
217 filter_S2
= load
<SimdBitMask
>(reinterpret_cast<const real
*>(exclusion_filter
+ 2 * UNROLLJ
));
221 /* Reaction-field constants */
222 mrc_3_S
= SimdReal(-2 * ic
->k_rf
);
223 # ifdef CALC_ENERGIES
224 hrc_3_S
= SimdReal(ic
->k_rf
);
225 moh_rc_S
= SimdReal(-ic
->c_rf
);
231 invtsp_S
= SimdReal(ic
->coulombEwaldTables
->scale
);
232 # ifdef CALC_ENERGIES
233 mhalfsp_S
= SimdReal(-0.5_real
/ ic
->coulombEwaldTables
->scale
);
237 tab_coul_F
= ic
->coulombEwaldTables
->tableFDV0
.data();
239 tab_coul_F
= ic
->coulombEwaldTables
->tableF
.data();
240 # ifdef CALC_ENERGIES
241 tab_coul_V
= ic
->coulombEwaldTables
->tableV
.data();
244 #endif /* CALC_COUL_TAB */
246 #ifdef CALC_COUL_EWALD
247 beta2_S
= SimdReal(ic
->ewaldcoeff_q
* ic
->ewaldcoeff_q
);
248 beta_S
= SimdReal(ic
->ewaldcoeff_q
);
251 #if (defined CALC_COUL_TAB || defined CALC_COUL_EWALD) && defined CALC_ENERGIES
252 sh_ewald_S
= SimdReal(ic
->sh_ewald
);
255 /* LJ function constants */
256 #if defined CALC_ENERGIES || defined LJ_POT_SWITCH
257 SimdReal sixth_S
= SimdReal(1.0 / 6.0);
258 SimdReal twelveth_S
= SimdReal(1.0 / 12.0);
261 #if defined LJ_CUT && defined CALC_ENERGIES
262 /* We shift the potential by cpot, which can be zero */
263 p6_cpot_S
= SimdReal(ic
->dispersion_shift
.cpot
);
264 p12_cpot_S
= SimdReal(ic
->repulsion_shift
.cpot
);
267 rswitch_S
= SimdReal(ic
->rvdw_switch
);
268 swV3_S
= SimdReal(ic
->vdw_switch
.c3
);
269 swV4_S
= SimdReal(ic
->vdw_switch
.c4
);
270 swV5_S
= SimdReal(ic
->vdw_switch
.c5
);
271 swF2_S
= SimdReal(3 * ic
->vdw_switch
.c3
);
272 swF3_S
= SimdReal(4 * ic
->vdw_switch
.c4
);
273 swF4_S
= SimdReal(5 * ic
->vdw_switch
.c5
);
275 #ifdef LJ_FORCE_SWITCH
276 rswitch_S
= SimdReal(ic
->rvdw_switch
);
277 p6_fc2_S
= SimdReal(ic
->dispersion_shift
.c2
);
278 p6_fc3_S
= SimdReal(ic
->dispersion_shift
.c3
);
279 p12_fc2_S
= SimdReal(ic
->repulsion_shift
.c2
);
280 p12_fc3_S
= SimdReal(ic
->repulsion_shift
.c3
);
281 # ifdef CALC_ENERGIES
283 SimdReal mthird_S
= SimdReal(-1.0 / 3.0);
284 SimdReal mfourth_S
= SimdReal(-1.0 / 4.0);
286 p6_vc3_S
= mthird_S
* p6_fc2_S
;
287 p6_vc4_S
= mfourth_S
* p6_fc3_S
;
288 p6_6cpot_S
= SimdReal(ic
->dispersion_shift
.cpot
/ 6);
289 p12_vc3_S
= mthird_S
* p12_fc2_S
;
290 p12_vc4_S
= mfourth_S
* p12_fc3_S
;
291 p12_12cpot_S
= SimdReal(ic
->repulsion_shift
.cpot
/ 12);
296 mone_S
= SimdReal(-1.0);
297 half_S
= SimdReal(0.5);
298 lj_ewaldcoeff2
= ic
->ewaldcoeff_lj
* ic
->ewaldcoeff_lj
;
299 lj_ewaldcoeff6_6
= lj_ewaldcoeff2
* lj_ewaldcoeff2
* lj_ewaldcoeff2
/ 6;
300 lje_c2_S
= SimdReal(lj_ewaldcoeff2
);
301 lje_c6_6_S
= SimdReal(lj_ewaldcoeff6_6
);
302 # ifdef CALC_ENERGIES
303 /* Determine the grid potential at the cut-off */
304 SimdReal lje_vc_S
= SimdReal(ic
->sh_lj_ewald
);
308 /* The kernel either supports rcoulomb = rvdw or rcoulomb >= rvdw */
309 rc2_S
= SimdReal(ic
->rcoulomb
* ic
->rcoulomb
);
310 #ifdef VDW_CUTOFF_CHECK
311 rcvdw2_S
= SimdReal(ic
->rvdw
* ic
->rvdw
);
314 minRsq_S
= SimdReal(NBNXN_MIN_RSQ
);
316 const real
* gmx_restrict q
= nbatParams
.q
.data();
317 const real facel
= ic
->epsfac
;
318 const real
* gmx_restrict shiftvec
= shift_vec
[0];
319 const real
* gmx_restrict x
= nbat
->x().data();
323 for (jp
= 0; jp
< UNROLLJ
; jp
++)
325 pvdw_c6
[0 * UNROLLJ
+ jp
] = nbat
->nbfp
[0 * 2];
326 pvdw_c6
[1 * UNROLLJ
+ jp
] = nbat
->nbfp
[0 * 2];
327 pvdw_c6
[2 * UNROLLJ
+ jp
] = nbat
->nbfp
[0 * 2];
328 pvdw_c6
[3 * UNROLLJ
+ jp
] = nbat
->nbfp
[0 * 2];
330 pvdw_c12
[0 * UNROLLJ
+ jp
] = nbat
->nbfp
[0 * 2 + 1];
331 pvdw_c12
[1 * UNROLLJ
+ jp
] = nbat
->nbfp
[0 * 2 + 1];
332 pvdw_c12
[2 * UNROLLJ
+ jp
] = nbat
->nbfp
[0 * 2 + 1];
333 pvdw_c12
[3 * UNROLLJ
+ jp
] = nbat
->nbfp
[0 * 2 + 1];
335 SimdReal c6_S0
= load
<SimdReal
>(pvdw_c6
+ 0 * UNROLLJ
);
336 SimdReal c6_S1
= load
<SimdReal
>(pvdw_c6
+ 1 * UNROLLJ
);
337 SimdReal c6_S2
= load
<SimdReal
>(pvdw_c6
+ 2 * UNROLLJ
);
338 SimdReal c6_S3
= load
<SimdReal
>(pvdw_c6
+ 3 * UNROLLJ
);
340 SimdReal c12_S0
= load
<SimdReal
>(pvdw_c12
+ 0 * UNROLLJ
);
341 SimdReal c12_S1
= load
<SimdReal
>(pvdw_c12
+ 1 * UNROLLJ
);
342 SimdReal c12_S2
= load
<SimdReal
>(pvdw_c12
+ 2 * UNROLLJ
);
343 SimdReal c12_S3
= load
<SimdReal
>(pvdw_c12
+ 3 * UNROLLJ
);
344 #endif /* FIX_LJ_C */
347 egps_ishift
= nbatParams
.neg_2log
;
348 egps_imask
= (1 << egps_ishift
) - 1;
349 egps_jshift
= 2 * nbatParams
.neg_2log
;
350 egps_jmask
= (1 << egps_jshift
) - 1;
351 egps_jstride
= (UNROLLJ
>> 1) * UNROLLJ
;
352 /* Major division is over i-particle energy groups, determine the stride */
353 Vstride_i
= nbatParams
.nenergrp
* (1 << nbatParams
.neg_2log
) * egps_jstride
;
356 l_cj
= nbl
->cj
.data();
359 for (const nbnxn_ci_t
& ciEntry
: nbl
->ci
)
361 ish
= (ciEntry
.shift
& NBNXN_CI_SHIFT
);
363 cjind0
= ciEntry
.cj_ind_start
;
364 cjind1
= ciEntry
.cj_ind_end
;
366 ci_sh
= (ish
== CENTRAL
? ci
: -1);
368 shX_S
= SimdReal(shiftvec
[ish3
]);
369 shY_S
= SimdReal(shiftvec
[ish3
+ 1]);
370 shZ_S
= SimdReal(shiftvec
[ish3
+ 2]);
373 int sci
= ci
* STRIDE
;
374 int scix
= sci
* DIM
;
375 # if defined LJ_COMB_LB || defined LJ_COMB_GEOM || defined LJ_EWALD_GEOM
379 int sci
= (ci
>> 1) * STRIDE
;
380 int scix
= sci
* DIM
+ (ci
& 1) * (STRIDE
>> 1);
381 # if defined LJ_COMB_LB || defined LJ_COMB_GEOM || defined LJ_EWALD_GEOM
382 int sci2
= sci
* 2 + (ci
& 1) * (STRIDE
>> 1);
384 sci
+= (ci
& 1) * (STRIDE
>> 1);
387 /* We have 5 LJ/C combinations, but use only three inner loops,
388 * as the other combinations are unlikely and/or not much faster:
389 * inner half-LJ + C for half-LJ + C / no-LJ + C
390 * inner LJ + C for full-LJ + C
391 * inner LJ for full-LJ + no-C / half-LJ + no-C
393 do_LJ
= ((ciEntry
.shift
& NBNXN_CI_DO_LJ(0)) != 0);
394 do_coul
= ((ciEntry
.shift
& NBNXN_CI_DO_COUL(0)) != 0);
395 half_LJ
= (((ciEntry
.shift
& NBNXN_CI_HALF_LJ(0)) != 0) || !do_LJ
) && do_coul
;
398 egps_i
= nbatParams
.energrp
[ci
];
402 for (ia
= 0; ia
< UNROLLI
; ia
++)
404 egp_ia
= (egps_i
>> (ia
* egps_ishift
)) & egps_imask
;
405 vvdwtp
[ia
] = Vvdw
+ egp_ia
* Vstride_i
;
406 vctp
[ia
] = Vc
+ egp_ia
* Vstride_i
;
412 # ifdef LJ_EWALD_GEOM
413 gmx_bool do_self
= TRUE
;
415 gmx_bool do_self
= do_coul
;
418 if (do_self
&& l_cj
[ciEntry
.cj_ind_start
].cj
== ci_sh
)
421 if (do_self
&& l_cj
[ciEntry
.cj_ind_start
].cj
== (ci_sh
>> 1))
430 Vc_sub_self
= 0.5 * ic
->c_rf
;
432 # ifdef CALC_COUL_TAB
434 Vc_sub_self
= 0.5 * tab_coul_F
[2];
436 Vc_sub_self
= 0.5 * tab_coul_V
[0];
439 # ifdef CALC_COUL_EWALD
441 Vc_sub_self
= 0.5 * ic
->ewaldcoeff_q
* M_2_SQRTPI
;
444 for (ia
= 0; ia
< UNROLLI
; ia
++)
449 # ifdef ENERGY_GROUPS
450 vctp
[ia
][((egps_i
>> (ia
* egps_ishift
)) & egps_imask
) * egps_jstride
]
454 -= facel
* qi
* qi
* Vc_sub_self
;
458 # ifdef LJ_EWALD_GEOM
462 for (ia
= 0; ia
< UNROLLI
; ia
++)
466 c6_i
= nbatParams
.nbfp
[nbatParams
.type
[sci
+ ia
] * (nbatParams
.numTypes
+ 1) * 2]
468 # ifdef ENERGY_GROUPS
469 vvdwtp
[ia
][((egps_i
>> (ia
* egps_ishift
)) & egps_imask
) * egps_jstride
]
473 += 0.5 * c6_i
* lj_ewaldcoeff6_6
;
476 # endif /* LJ_EWALD */
480 /* Load i atom data */
481 int sciy
= scix
+ STRIDE
;
482 int sciz
= sciy
+ STRIDE
;
483 ix_S0
= loadU1DualHsimd(x
+ scix
);
484 ix_S2
= loadU1DualHsimd(x
+ scix
+ 2);
485 iy_S0
= loadU1DualHsimd(x
+ sciy
);
486 iy_S2
= loadU1DualHsimd(x
+ sciy
+ 2);
487 iz_S0
= loadU1DualHsimd(x
+ sciz
);
488 iz_S2
= loadU1DualHsimd(x
+ sciz
+ 2);
489 ix_S0
= ix_S0
+ shX_S
;
490 ix_S2
= ix_S2
+ shX_S
;
491 iy_S0
= iy_S0
+ shY_S
;
492 iy_S2
= iy_S2
+ shY_S
;
493 iz_S0
= iz_S0
+ shZ_S
;
494 iz_S2
= iz_S2
+ shZ_S
;
500 facel_S
= SimdReal(facel
);
502 iq_S0
= loadU1DualHsimd(q
+ sci
);
503 iq_S2
= loadU1DualHsimd(q
+ sci
+ 2);
504 iq_S0
= facel_S
* iq_S0
;
505 iq_S2
= facel_S
* iq_S2
;
509 hsig_i_S0
= loadU1DualHsimd(ljc
+ sci2
);
510 hsig_i_S2
= loadU1DualHsimd(ljc
+ sci2
+ 2);
511 seps_i_S0
= loadU1DualHsimd(ljc
+ sci2
+ STRIDE
);
512 seps_i_S2
= loadU1DualHsimd(ljc
+ sci2
+ STRIDE
+ 2);
515 SimdReal c6s_S0
, c12s_S0
;
516 SimdReal c6s_S2
, c12s_S2
;
518 c6s_S0
= loadU1DualHsimd(ljc
+ sci2
);
522 c6s_S2
= loadU1DualHsimd(ljc
+ sci2
+ 2);
524 c12s_S0
= loadU1DualHsimd(ljc
+ sci2
+ STRIDE
);
527 c12s_S2
= loadU1DualHsimd(ljc
+ sci2
+ STRIDE
+ 2);
529 # elif !defined LJ_COMB_LB && !defined FIX_LJ_C
530 const int numTypes
= nbatParams
.numTypes
;
531 const real
* nbfp0
= nbfp_ptr
+ type
[sci
] * numTypes
* c_simdBestPairAlignment
;
532 const real
* nbfp1
= nbfp_ptr
+ type
[sci
+ 1] * numTypes
* c_simdBestPairAlignment
;
533 const real
*nbfp2
= nullptr, *nbfp3
= nullptr;
536 nbfp2
= nbfp_ptr
+ type
[sci
+ 2] * numTypes
* c_simdBestPairAlignment
;
537 nbfp3
= nbfp_ptr
+ type
[sci
+ 3] * numTypes
* c_simdBestPairAlignment
;
542 /* We need the geometrically combined C6 for the PME grid correction */
543 SimdReal c6s_S0
, c6s_S2
;
544 c6s_S0
= loadU1DualHsimd(ljc
+ sci2
);
547 c6s_S2
= loadU1DualHsimd(ljc
+ sci2
+ 2);
551 /* Zero the potential energy for this list */
553 SimdReal Vvdwtot_S
= setZero();
554 SimdReal vctot_S
= setZero();
557 /* Clear i atom forces */
567 /* Currently all kernels use (at least half) LJ */
571 /* Coulomb: all i-atoms, LJ: first half i-atoms */
575 while (cjind
< cjind1
&& nbl
->cj
[cjind
].excl
!= NBNXN_INTERACTION_MASK_ALL
)
577 #include "kernel_inner.h"
581 for (; (cjind
< cjind1
); cjind
++)
583 #include "kernel_inner.h"
590 /* Coulomb: all i-atoms, LJ: all i-atoms */
593 while (cjind
< cjind1
&& nbl
->cj
[cjind
].excl
!= NBNXN_INTERACTION_MASK_ALL
)
595 #include "kernel_inner.h"
599 for (; (cjind
< cjind1
); cjind
++)
601 #include "kernel_inner.h"
607 /* Coulomb: none, LJ: all i-atoms */
609 while (cjind
< cjind1
&& nbl
->cj
[cjind
].excl
!= NBNXN_INTERACTION_MASK_ALL
)
611 #include "kernel_inner.h"
615 for (; (cjind
< cjind1
); cjind
++)
617 #include "kernel_inner.h"
621 ninner
+= cjind1
- cjind0
;
623 /* Add accumulated i-forces to the force array */
624 real fShiftX
= reduceIncr4ReturnSumHsimd(f
+ scix
, fix_S0
, fix_S2
);
625 real fShiftY
= reduceIncr4ReturnSumHsimd(f
+ sciy
, fiy_S0
, fiy_S2
);
626 real fShiftZ
= reduceIncr4ReturnSumHsimd(f
+ sciz
, fiz_S0
, fiz_S2
);
628 #ifdef CALC_SHIFTFORCES
629 fshift
[ish3
+ 0] += fShiftX
;
630 fshift
[ish3
+ 1] += fShiftY
;
631 fshift
[ish3
+ 2] += fShiftZ
;
637 *Vc
+= reduce(vctot_S
);
639 *Vvdw
+= reduce(Vvdwtot_S
);
642 /* Outer loop uses 6 flops/iteration */
646 printf("atom pairs %d\n", npair
);