Merge branch 'release-4-6'
[alexxy/gromacs.git] / src / gromacs / mdlib / nbnxn_kernels / simd_4xn / nbnxn_kernel_simd_4xn_outer.h
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2012,2013, 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.
8  *
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.
13  *
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.
18  *
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.
23  *
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.
31  *
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.
34  */
35
36 {
37     const nbnxn_ci_t   *nbln;
38     const nbnxn_cj_t   *l_cj;
39     const int          *type;
40     const real         *q;
41     const real         *shiftvec;
42     const real         *x;
43     const real         *nbfp0, *nbfp1, *nbfp2 = NULL, *nbfp3 = NULL;
44     real                facel;
45     real               *nbfp_ptr;
46     int                 n, ci, ci_sh;
47     int                 ish, ish3;
48     gmx_bool            do_LJ, half_LJ, do_coul;
49     int                 sci, scix, sciy, sciz, sci2;
50     int                 cjind0, cjind1, cjind;
51     int                 ip, jp;
52
53 #ifdef ENERGY_GROUPS
54     int         Vstride_i;
55     int         egps_ishift, egps_imask;
56     int         egps_jshift, egps_jmask, egps_jstride;
57     int         egps_i;
58     real       *vvdwtp[UNROLLI];
59     real       *vctp[UNROLLI];
60 #endif
61
62     gmx_mm_pr  shX_S;
63     gmx_mm_pr  shY_S;
64     gmx_mm_pr  shZ_S;
65     gmx_mm_pr  ix_S0, iy_S0, iz_S0;
66     gmx_mm_pr  ix_S1, iy_S1, iz_S1;
67     gmx_mm_pr  ix_S2, iy_S2, iz_S2;
68     gmx_mm_pr  ix_S3, iy_S3, iz_S3;
69     gmx_mm_pr  fix_S0, fiy_S0, fiz_S0;
70     gmx_mm_pr  fix_S1, fiy_S1, fiz_S1;
71     gmx_mm_pr  fix_S2, fiy_S2, fiz_S2;
72     gmx_mm_pr  fix_S3, fiy_S3, fiz_S3;
73 #if UNROLLJ >= 4
74     /* We use an i-force SIMD register width of 4 */
75     gmx_mm_pr4 fix_S, fiy_S, fiz_S;
76 #else
77     /* We use an i-force SIMD register width of 2 */
78     gmx_mm_pr  fix0_S, fiy0_S, fiz0_S;
79     gmx_mm_pr  fix2_S, fiy2_S, fiz2_S;
80 #endif
81
82     gmx_mm_pr  diagonal_jmi_S;
83 #if UNROLLI == UNROLLJ
84     gmx_mm_pb  diagonal_mask_S0, diagonal_mask_S1, diagonal_mask_S2, diagonal_mask_S3;
85 #else
86     gmx_mm_pb  diagonal_mask0_S0, diagonal_mask0_S1, diagonal_mask0_S2, diagonal_mask0_S3;
87     gmx_mm_pb  diagonal_mask1_S0, diagonal_mask1_S1, diagonal_mask1_S2, diagonal_mask1_S3;
88 #endif
89
90     unsigned      *exclusion_filter;
91     gmx_exclfilter filter_S0, filter_S1, filter_S2, filter_S3;
92
93     gmx_mm_pr      zero_S = gmx_set1_pr(0.0);
94
95     gmx_mm_pr      one_S  = gmx_set1_pr(1.0);
96     gmx_mm_pr      iq_S0  = gmx_setzero_pr();
97     gmx_mm_pr      iq_S1  = gmx_setzero_pr();
98     gmx_mm_pr      iq_S2  = gmx_setzero_pr();
99     gmx_mm_pr      iq_S3  = gmx_setzero_pr();
100     gmx_mm_pr      mrc_3_S;
101 #ifdef CALC_ENERGIES
102     gmx_mm_pr      hrc_3_S, moh_rc_S;
103 #endif
104
105 #ifdef CALC_COUL_TAB
106     /* Coulomb table variables */
107     gmx_mm_pr   invtsp_S;
108     const real *tab_coul_F;
109 #ifndef TAB_FDV0
110     const real *tab_coul_V;
111 #endif
112     /* Thread-local working buffers for force and potential lookups */
113     int         ti0_array[2*GMX_SIMD_WIDTH_HERE], *ti0 = NULL;
114     int         ti1_array[2*GMX_SIMD_WIDTH_HERE], *ti1 = NULL;
115     int         ti2_array[2*GMX_SIMD_WIDTH_HERE], *ti2 = NULL;
116     int         ti3_array[2*GMX_SIMD_WIDTH_HERE], *ti3 = NULL;
117 #ifdef CALC_ENERGIES
118     gmx_mm_pr   mhalfsp_S;
119 #endif
120 #endif
121
122 #ifdef CALC_COUL_EWALD
123     gmx_mm_pr beta2_S, beta_S;
124 #endif
125
126 #if defined CALC_ENERGIES && (defined CALC_COUL_EWALD || defined CALC_COUL_TAB)
127     gmx_mm_pr  sh_ewald_S;
128 #endif
129
130 #ifdef LJ_COMB_LB
131     const real *ljc;
132
133     gmx_mm_pr   hsig_i_S0, seps_i_S0;
134     gmx_mm_pr   hsig_i_S1, seps_i_S1;
135     gmx_mm_pr   hsig_i_S2, seps_i_S2;
136     gmx_mm_pr   hsig_i_S3, seps_i_S3;
137 #else
138 #ifdef FIX_LJ_C
139     real        pvdw_array[2*UNROLLI*UNROLLJ+3];
140     real       *pvdw_c6, *pvdw_c12;
141     gmx_mm_pr   c6_S0, c12_S0;
142     gmx_mm_pr   c6_S1, c12_S1;
143     gmx_mm_pr   c6_S2, c12_S2;
144     gmx_mm_pr   c6_S3, c12_S3;
145 #endif
146
147 #ifdef LJ_COMB_GEOM
148     const real *ljc;
149
150     gmx_mm_pr   c6s_S0, c12s_S0;
151     gmx_mm_pr   c6s_S1, c12s_S1;
152     gmx_mm_pr   c6s_S2 = gmx_setzero_pr(), c12s_S2 = gmx_setzero_pr();
153     gmx_mm_pr   c6s_S3 = gmx_setzero_pr(), c12s_S3 = gmx_setzero_pr();
154 #endif
155 #endif /* LJ_COMB_LB */
156
157     gmx_mm_pr  vctot_S, Vvdwtot_S;
158     gmx_mm_pr  sixth_S, twelveth_S;
159
160     gmx_mm_pr  avoid_sing_S;
161     gmx_mm_pr  rc2_S;
162 #ifdef VDW_CUTOFF_CHECK
163     gmx_mm_pr  rcvdw2_S;
164 #endif
165
166 #ifdef CALC_ENERGIES
167     gmx_mm_pr  sh_invrc6_S, sh_invrc12_S;
168
169     /* cppcheck-suppress unassignedVariable */
170     real       tmpsum_array[15], *tmpsum;
171 #endif
172 #ifdef CALC_SHIFTFORCES
173     /* cppcheck-suppress unassignedVariable */
174     real       shf_array[15], *shf;
175 #endif
176
177     int ninner;
178
179 #ifdef COUNT_PAIRS
180     int npair = 0;
181 #endif
182
183 #if defined LJ_COMB_GEOM || defined LJ_COMB_LB
184     ljc = nbat->lj_comb;
185 #else
186     /* No combination rule used */
187     nbfp_ptr    = (4 == nbfp_stride) ? nbat->nbfp_s4 : nbat->nbfp;
188 #endif
189
190     /* Load j-i for the first i */
191     diagonal_jmi_S    = gmx_load_pr(nbat->simd_4xn_diagonal_j_minus_i);
192     /* Generate all the diagonal masks as comparison results */
193 #if UNROLLI == UNROLLJ
194     diagonal_mask_S0  = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
195     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
196     diagonal_mask_S1  = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
197     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
198     diagonal_mask_S2  = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
199     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
200     diagonal_mask_S3  = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
201 #else
202 #if UNROLLI == 2*UNROLLJ || 2*UNROLLI == UNROLLJ
203     diagonal_mask0_S0 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
204     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
205     diagonal_mask0_S1 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
206     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
207     diagonal_mask0_S2 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
208     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
209     diagonal_mask0_S3 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
210     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
211
212 #if UNROLLI == 2*UNROLLJ
213     /* Load j-i for the second half of the j-cluster */
214     diagonal_jmi_S    = gmx_load_pr(nbat->simd_4xn_diagonal_j_minus_i + UNROLLJ);
215 #endif
216
217     diagonal_mask1_S0 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
218     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
219     diagonal_mask1_S1 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
220     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
221     diagonal_mask1_S2 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
222     diagonal_jmi_S    = gmx_sub_pr(diagonal_jmi_S, one_S);
223     diagonal_mask1_S3 = gmx_cmplt_pr(zero_S, diagonal_jmi_S);
224 #endif
225 #endif
226
227     /* Load masks for topology exclusion masking. filter_stride is
228        static const, so the conditional will be optimized away. */
229     if (1 == filter_stride)
230     {
231         exclusion_filter = nbat->simd_exclusion_filter1;
232     }
233     else /* (2 == filter_stride) */
234     {
235         exclusion_filter = nbat->simd_exclusion_filter2;
236     }
237
238     /* Here we cast the exclusion filters from unsigned * to int * or real *.
239      * Since we only check bits, the actual value they represent does not
240      * matter, as long as both filter and mask data are treated the same way.
241      */
242     filter_S0    = gmx_load_exclusion_filter(exclusion_filter + 0*UNROLLJ*filter_stride);
243     filter_S1    = gmx_load_exclusion_filter(exclusion_filter + 1*UNROLLJ*filter_stride);
244     filter_S2    = gmx_load_exclusion_filter(exclusion_filter + 2*UNROLLJ*filter_stride);
245     filter_S3    = gmx_load_exclusion_filter(exclusion_filter + 3*UNROLLJ*filter_stride);
246
247 #ifdef CALC_COUL_TAB
248     /* Generate aligned table index pointers */
249     ti0 = prepare_table_load_buffer(ti0_array);
250     ti1 = prepare_table_load_buffer(ti1_array);
251     ti2 = prepare_table_load_buffer(ti2_array);
252     ti3 = prepare_table_load_buffer(ti3_array);
253
254     invtsp_S  = gmx_set1_pr(ic->tabq_scale);
255 #ifdef CALC_ENERGIES
256     mhalfsp_S = gmx_set1_pr(-0.5/ic->tabq_scale);
257 #endif
258
259 #ifdef TAB_FDV0
260     tab_coul_F = ic->tabq_coul_FDV0;
261 #else
262     tab_coul_F = ic->tabq_coul_F;
263     tab_coul_V = ic->tabq_coul_V;
264 #endif
265 #endif /* CALC_COUL_TAB */
266
267 #ifdef CALC_COUL_EWALD
268     beta2_S = gmx_set1_pr(ic->ewaldcoeff*ic->ewaldcoeff);
269     beta_S  = gmx_set1_pr(ic->ewaldcoeff);
270 #endif
271
272 #if (defined CALC_COUL_TAB || defined CALC_COUL_EWALD) && defined CALC_ENERGIES
273     sh_ewald_S = gmx_set1_pr(ic->sh_ewald);
274 #endif
275
276     q                   = nbat->q;
277     type                = nbat->type;
278     facel               = ic->epsfac;
279     shiftvec            = shift_vec[0];
280     x                   = nbat->x;
281
282     avoid_sing_S = gmx_set1_pr(NBNXN_AVOID_SING_R2_INC);
283
284     /* The kernel either supports rcoulomb = rvdw or rcoulomb >= rvdw */
285     rc2_S    = gmx_set1_pr(ic->rcoulomb*ic->rcoulomb);
286 #ifdef VDW_CUTOFF_CHECK
287     rcvdw2_S = gmx_set1_pr(ic->rvdw*ic->rvdw);
288 #endif
289
290 #ifdef CALC_ENERGIES
291     sixth_S      = gmx_set1_pr(1.0/6.0);
292     twelveth_S   = gmx_set1_pr(1.0/12.0);
293
294     sh_invrc6_S  = gmx_set1_pr(ic->sh_invrc6);
295     sh_invrc12_S = gmx_set1_pr(ic->sh_invrc6*ic->sh_invrc6);
296 #endif
297
298     mrc_3_S  = gmx_set1_pr(-2*ic->k_rf);
299
300 #ifdef CALC_ENERGIES
301     hrc_3_S  = gmx_set1_pr(ic->k_rf);
302
303     moh_rc_S = gmx_set1_pr(-ic->c_rf);
304 #endif
305
306 #ifdef CALC_ENERGIES
307     tmpsum   = gmx_simd_align_real(tmpsum_array);
308 #endif
309 #ifdef CALC_SHIFTFORCES
310     shf      = gmx_simd_align_real(shf_array);
311 #endif
312
313 #ifdef FIX_LJ_C
314     pvdw_c6  = gmx_simd_align_real(pvdw_array+3);
315     pvdw_c12 = pvdw_c6 + UNROLLI*UNROLLJ;
316
317     for (jp = 0; jp < UNROLLJ; jp++)
318     {
319         pvdw_c6 [0*UNROLLJ+jp] = nbat->nbfp[0*2];
320         pvdw_c6 [1*UNROLLJ+jp] = nbat->nbfp[0*2];
321         pvdw_c6 [2*UNROLLJ+jp] = nbat->nbfp[0*2];
322         pvdw_c6 [3*UNROLLJ+jp] = nbat->nbfp[0*2];
323
324         pvdw_c12[0*UNROLLJ+jp] = nbat->nbfp[0*2+1];
325         pvdw_c12[1*UNROLLJ+jp] = nbat->nbfp[0*2+1];
326         pvdw_c12[2*UNROLLJ+jp] = nbat->nbfp[0*2+1];
327         pvdw_c12[3*UNROLLJ+jp] = nbat->nbfp[0*2+1];
328     }
329     c6_S0            = gmx_load_pr(pvdw_c6 +0*UNROLLJ);
330     c6_S1            = gmx_load_pr(pvdw_c6 +1*UNROLLJ);
331     c6_S2            = gmx_load_pr(pvdw_c6 +2*UNROLLJ);
332     c6_S3            = gmx_load_pr(pvdw_c6 +3*UNROLLJ);
333
334     c12_S0           = gmx_load_pr(pvdw_c12+0*UNROLLJ);
335     c12_S1           = gmx_load_pr(pvdw_c12+1*UNROLLJ);
336     c12_S2           = gmx_load_pr(pvdw_c12+2*UNROLLJ);
337     c12_S3           = gmx_load_pr(pvdw_c12+3*UNROLLJ);
338 #endif /* FIX_LJ_C */
339
340 #ifdef ENERGY_GROUPS
341     egps_ishift  = nbat->neg_2log;
342     egps_imask   = (1<<egps_ishift) - 1;
343     egps_jshift  = 2*nbat->neg_2log;
344     egps_jmask   = (1<<egps_jshift) - 1;
345     egps_jstride = (UNROLLJ>>1)*UNROLLJ;
346     /* Major division is over i-particle energy groups, determine the stride */
347     Vstride_i    = nbat->nenergrp*(1<<nbat->neg_2log)*egps_jstride;
348 #endif
349
350     l_cj = nbl->cj;
351
352     ninner = 0;
353     for (n = 0; n < nbl->nci; n++)
354     {
355         nbln = &nbl->ci[n];
356
357         ish              = (nbln->shift & NBNXN_CI_SHIFT);
358         ish3             = ish*3;
359         cjind0           = nbln->cj_ind_start;
360         cjind1           = nbln->cj_ind_end;
361         ci               = nbln->ci;
362         ci_sh            = (ish == CENTRAL ? ci : -1);
363
364         shX_S = gmx_load1_pr(shiftvec+ish3);
365         shY_S = gmx_load1_pr(shiftvec+ish3+1);
366         shZ_S = gmx_load1_pr(shiftvec+ish3+2);
367
368 #if UNROLLJ <= 4
369         sci              = ci*STRIDE;
370         scix             = sci*DIM;
371         sci2             = sci*2;
372 #else
373         sci              = (ci>>1)*STRIDE;
374         scix             = sci*DIM + (ci & 1)*(STRIDE>>1);
375         sci2             = sci*2 + (ci & 1)*(STRIDE>>1);
376         sci             += (ci & 1)*(STRIDE>>1);
377 #endif
378
379         /* We have 5 LJ/C combinations, but use only three inner loops,
380          * as the other combinations are unlikely and/or not much faster:
381          * inner half-LJ + C for half-LJ + C / no-LJ + C
382          * inner LJ + C      for full-LJ + C
383          * inner LJ          for full-LJ + no-C / half-LJ + no-C
384          */
385         do_LJ   = (nbln->shift & NBNXN_CI_DO_LJ(0));
386         do_coul = (nbln->shift & NBNXN_CI_DO_COUL(0));
387         half_LJ = ((nbln->shift & NBNXN_CI_HALF_LJ(0)) || !do_LJ) && do_coul;
388
389 #ifdef ENERGY_GROUPS
390         egps_i = nbat->energrp[ci];
391         {
392             int ia, egp_ia;
393
394             for (ia = 0; ia < UNROLLI; ia++)
395             {
396                 egp_ia     = (egps_i >> (ia*egps_ishift)) & egps_imask;
397                 vvdwtp[ia] = Vvdw + egp_ia*Vstride_i;
398                 vctp[ia]   = Vc   + egp_ia*Vstride_i;
399             }
400         }
401 #endif
402 #if defined CALC_ENERGIES
403 #if UNROLLJ == 4
404         if (do_coul && l_cj[nbln->cj_ind_start].cj == ci_sh)
405 #endif
406 #if UNROLLJ == 2
407         if (do_coul && l_cj[nbln->cj_ind_start].cj == (ci_sh<<1))
408 #endif
409 #if UNROLLJ == 8
410         if (do_coul && l_cj[nbln->cj_ind_start].cj == (ci_sh>>1))
411 #endif
412         {
413             int  ia;
414             real Vc_sub_self;
415
416 #ifdef CALC_COUL_RF
417             Vc_sub_self = 0.5*ic->c_rf;
418 #endif
419 #ifdef CALC_COUL_TAB
420 #ifdef TAB_FDV0
421             Vc_sub_self = 0.5*tab_coul_F[2];
422 #else
423             Vc_sub_self = 0.5*tab_coul_V[0];
424 #endif
425 #endif
426 #ifdef CALC_COUL_EWALD
427             /* beta/sqrt(pi) */
428             Vc_sub_self = 0.5*ic->ewaldcoeff*M_2_SQRTPI;
429 #endif
430
431             for (ia = 0; ia < UNROLLI; ia++)
432             {
433                 real qi;
434
435                 qi = q[sci+ia];
436 #ifdef ENERGY_GROUPS
437                 vctp[ia][((egps_i>>(ia*egps_ishift)) & egps_imask)*egps_jstride]
438 #else
439                 Vc[0]
440 #endif
441                     -= facel*qi*qi*Vc_sub_self;
442             }
443         }
444 #endif
445
446         /* Load i atom data */
447         sciy             = scix + STRIDE;
448         sciz             = sciy + STRIDE;
449         ix_S0            = gmx_add_pr(gmx_load1_pr(x+scix), shX_S);
450         ix_S1            = gmx_add_pr(gmx_load1_pr(x+scix+1), shX_S);
451         ix_S2            = gmx_add_pr(gmx_load1_pr(x+scix+2), shX_S);
452         ix_S3            = gmx_add_pr(gmx_load1_pr(x+scix+3), shX_S);
453         iy_S0            = gmx_add_pr(gmx_load1_pr(x+sciy), shY_S);
454         iy_S1            = gmx_add_pr(gmx_load1_pr(x+sciy+1), shY_S);
455         iy_S2            = gmx_add_pr(gmx_load1_pr(x+sciy+2), shY_S);
456         iy_S3            = gmx_add_pr(gmx_load1_pr(x+sciy+3), shY_S);
457         iz_S0            = gmx_add_pr(gmx_load1_pr(x+sciz), shZ_S);
458         iz_S1            = gmx_add_pr(gmx_load1_pr(x+sciz+1), shZ_S);
459         iz_S2            = gmx_add_pr(gmx_load1_pr(x+sciz+2), shZ_S);
460         iz_S3            = gmx_add_pr(gmx_load1_pr(x+sciz+3), shZ_S);
461
462         if (do_coul)
463         {
464             iq_S0      = gmx_set1_pr(facel*q[sci]);
465             iq_S1      = gmx_set1_pr(facel*q[sci+1]);
466             iq_S2      = gmx_set1_pr(facel*q[sci+2]);
467             iq_S3      = gmx_set1_pr(facel*q[sci+3]);
468         }
469
470 #ifdef LJ_COMB_LB
471         hsig_i_S0      = gmx_load1_pr(ljc+sci2+0);
472         hsig_i_S1      = gmx_load1_pr(ljc+sci2+1);
473         hsig_i_S2      = gmx_load1_pr(ljc+sci2+2);
474         hsig_i_S3      = gmx_load1_pr(ljc+sci2+3);
475         seps_i_S0      = gmx_load1_pr(ljc+sci2+STRIDE+0);
476         seps_i_S1      = gmx_load1_pr(ljc+sci2+STRIDE+1);
477         seps_i_S2      = gmx_load1_pr(ljc+sci2+STRIDE+2);
478         seps_i_S3      = gmx_load1_pr(ljc+sci2+STRIDE+3);
479 #else
480 #ifdef LJ_COMB_GEOM
481         c6s_S0         = gmx_load1_pr(ljc+sci2+0);
482         c6s_S1         = gmx_load1_pr(ljc+sci2+1);
483         if (!half_LJ)
484         {
485             c6s_S2     = gmx_load1_pr(ljc+sci2+2);
486             c6s_S3     = gmx_load1_pr(ljc+sci2+3);
487         }
488         c12s_S0        = gmx_load1_pr(ljc+sci2+STRIDE+0);
489         c12s_S1        = gmx_load1_pr(ljc+sci2+STRIDE+1);
490         if (!half_LJ)
491         {
492             c12s_S2    = gmx_load1_pr(ljc+sci2+STRIDE+2);
493             c12s_S3    = gmx_load1_pr(ljc+sci2+STRIDE+3);
494         }
495 #else
496         nbfp0     = nbfp_ptr + type[sci  ]*nbat->ntype*nbfp_stride;
497         nbfp1     = nbfp_ptr + type[sci+1]*nbat->ntype*nbfp_stride;
498         if (!half_LJ)
499         {
500             nbfp2 = nbfp_ptr + type[sci+2]*nbat->ntype*nbfp_stride;
501             nbfp3 = nbfp_ptr + type[sci+3]*nbat->ntype*nbfp_stride;
502         }
503 #endif
504 #endif
505
506         /* Zero the potential energy for this list */
507         Vvdwtot_S        = gmx_setzero_pr();
508         vctot_S          = gmx_setzero_pr();
509
510         /* Clear i atom forces */
511         fix_S0           = gmx_setzero_pr();
512         fix_S1           = gmx_setzero_pr();
513         fix_S2           = gmx_setzero_pr();
514         fix_S3           = gmx_setzero_pr();
515         fiy_S0           = gmx_setzero_pr();
516         fiy_S1           = gmx_setzero_pr();
517         fiy_S2           = gmx_setzero_pr();
518         fiy_S3           = gmx_setzero_pr();
519         fiz_S0           = gmx_setzero_pr();
520         fiz_S1           = gmx_setzero_pr();
521         fiz_S2           = gmx_setzero_pr();
522         fiz_S3           = gmx_setzero_pr();
523
524         cjind = cjind0;
525
526         /* Currently all kernels use (at least half) LJ */
527 #define CALC_LJ
528         if (half_LJ)
529         {
530 #define CALC_COULOMB
531 #define HALF_LJ
532 #define CHECK_EXCLS
533             while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
534             {
535 #include "nbnxn_kernel_simd_4xn_inner.h"
536                 cjind++;
537             }
538 #undef CHECK_EXCLS
539             for (; (cjind < cjind1); cjind++)
540             {
541 #include "nbnxn_kernel_simd_4xn_inner.h"
542             }
543 #undef HALF_LJ
544 #undef CALC_COULOMB
545         }
546         else if (do_coul)
547         {
548 #define CALC_COULOMB
549 #define CHECK_EXCLS
550             while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
551             {
552 #include "nbnxn_kernel_simd_4xn_inner.h"
553                 cjind++;
554             }
555 #undef CHECK_EXCLS
556             for (; (cjind < cjind1); cjind++)
557             {
558 #include "nbnxn_kernel_simd_4xn_inner.h"
559             }
560 #undef CALC_COULOMB
561         }
562         else
563         {
564 #define CHECK_EXCLS
565             while (cjind < cjind1 && nbl->cj[cjind].excl != NBNXN_INTERACTION_MASK_ALL)
566             {
567 #include "nbnxn_kernel_simd_4xn_inner.h"
568                 cjind++;
569             }
570 #undef CHECK_EXCLS
571             for (; (cjind < cjind1); cjind++)
572             {
573 #include "nbnxn_kernel_simd_4xn_inner.h"
574             }
575         }
576 #undef CALC_LJ
577         ninner += cjind1 - cjind0;
578
579         /* Add accumulated i-forces to the force array */
580 #if UNROLLJ >= 4
581         fix_S = gmx_mm_transpose_sum4_pr(fix_S0, fix_S1, fix_S2, fix_S3);
582         gmx_store_pr4(f+scix, gmx_add_pr4(fix_S, gmx_load_pr4(f+scix)));
583
584         fiy_S = gmx_mm_transpose_sum4_pr(fiy_S0, fiy_S1, fiy_S2, fiy_S3);
585         gmx_store_pr4(f+sciy, gmx_add_pr4(fiy_S, gmx_load_pr4(f+sciy)));
586
587         fiz_S = gmx_mm_transpose_sum4_pr(fiz_S0, fiz_S1, fiz_S2, fiz_S3);
588         gmx_store_pr4(f+sciz, gmx_add_pr4(fiz_S, gmx_load_pr4(f+sciz)));
589
590 #ifdef CALC_SHIFTFORCES
591         fshift[ish3+0] += gmx_sum_simd4(fix_S, shf);
592         fshift[ish3+1] += gmx_sum_simd4(fiy_S, shf);
593         fshift[ish3+2] += gmx_sum_simd4(fiz_S, shf);
594 #endif
595 #else
596         fix0_S = gmx_mm_transpose_sum2_pr(fix_S0, fix_S1);
597         gmx_store_pr(f+scix, gmx_add_pr(fix0_S, gmx_load_pr(f+scix)));
598         fix2_S = gmx_mm_transpose_sum2_pr(fix_S2, fix_S3);
599         gmx_store_pr(f+scix+2, gmx_add_pr(fix2_S, gmx_load_pr(f+scix+2)));
600
601         fiy0_S = gmx_mm_transpose_sum2_pr(fiy_S0, fiy_S1);
602         gmx_store_pr(f+sciy, gmx_add_pr(fiy0_S, gmx_load_pr(f+sciy)));
603         fiy2_S = gmx_mm_transpose_sum2_pr(fiy_S2, fiy_S3);
604         gmx_store_pr(f+sciy+2, gmx_add_pr(fiy2_S, gmx_load_pr(f+sciy+2)));
605
606         fiz0_S = gmx_mm_transpose_sum2_pr(fiz_S0, fiz_S1);
607         gmx_store_pr(f+sciz, gmx_add_pr(fiz0_S, gmx_load_pr(f+sciz)));
608         fiz2_S = gmx_mm_transpose_sum2_pr(fiz_S2, fiz_S3);
609         gmx_store_pr(f+sciz+2, gmx_add_pr(fiz2_S, gmx_load_pr(f+sciz+2)));
610
611 #ifdef CALC_SHIFTFORCES
612         fshift[ish3+0] += gmx_sum_simd2(gmx_add_pr(fix0_S, fix2_S), shf);
613         fshift[ish3+1] += gmx_sum_simd2(gmx_add_pr(fiy0_S, fiy2_S), shf);
614         fshift[ish3+2] += gmx_sum_simd2(gmx_add_pr(fiz0_S, fiz2_S), shf);
615 #endif
616 #endif
617
618 #ifdef CALC_ENERGIES
619         if (do_coul)
620         {
621             *Vc += gmx_sum_simd(vctot_S, tmpsum);
622         }
623
624         *Vvdw += gmx_sum_simd(Vvdwtot_S, tmpsum);
625 #endif
626
627         /* Outer loop uses 6 flops/iteration */
628     }
629
630 #ifdef COUNT_PAIRS
631     printf("atom pairs %d\n", npair);
632 #endif
633 }