Merge release-4-6 into master
[alexxy/gromacs.git] / src / gromacs / gmxlib / nonbonded / nb_kernel_avx_256_single / nb_kernel_ElecRFCut_VdwNone_GeomW4P1_avx_256_single.c
1 /*
2  * Note: this file was generated by the Gromacs avx_256_single kernel generator.
3  *
4  *                This source code is part of
5  *
6  *                 G   R   O   M   A   C   S
7  *
8  * Copyright (c) 2001-2012, The GROMACS Development Team
9  *
10  * Gromacs is a library for molecular simulation and trajectory analysis,
11  * written by Erik Lindahl, David van der Spoel, Berk Hess, and others - for
12  * a full list of developers and information, check out http://www.gromacs.org
13  *
14  * This program is free software; you can redistribute it and/or modify it under
15  * the terms of the GNU Lesser General Public License as published by the Free
16  * Software Foundation; either version 2 of the License, or (at your option) any
17  * later version.
18  *
19  * To help fund GROMACS development, we humbly ask that you cite
20  * the papers people have written on it - you can find them on the website.
21  */
22 #ifdef HAVE_CONFIG_H
23 #include <config.h>
24 #endif
25
26 #include <math.h>
27
28 #include "../nb_kernel.h"
29 #include "types/simple.h"
30 #include "vec.h"
31 #include "nrnb.h"
32
33 #include "gmx_math_x86_avx_256_single.h"
34 #include "kernelutil_x86_avx_256_single.h"
35
36 /*
37  * Gromacs nonbonded kernel:   nb_kernel_ElecRFCut_VdwNone_GeomW4P1_VF_avx_256_single
38  * Electrostatics interaction: ReactionField
39  * VdW interaction:            None
40  * Geometry:                   Water4-Particle
41  * Calculate force/pot:        PotentialAndForce
42  */
43 void
44 nb_kernel_ElecRFCut_VdwNone_GeomW4P1_VF_avx_256_single
45                     (t_nblist * gmx_restrict                nlist,
46                      rvec * gmx_restrict                    xx,
47                      rvec * gmx_restrict                    ff,
48                      t_forcerec * gmx_restrict              fr,
49                      t_mdatoms * gmx_restrict               mdatoms,
50                      nb_kernel_data_t * gmx_restrict        kernel_data,
51                      t_nrnb * gmx_restrict                  nrnb)
52 {
53     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or 
54      * just 0 for non-waters.
55      * Suffixes A,B,C,D,E,F,G,H refer to j loop unrolling done with AVX, e.g. for the eight different
56      * jnr indices corresponding to data put in the four positions in the SIMD register.
57      */
58     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
59     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
60     int              jnrA,jnrB,jnrC,jnrD;
61     int              jnrE,jnrF,jnrG,jnrH;
62     int              jnrlistA,jnrlistB,jnrlistC,jnrlistD;
63     int              jnrlistE,jnrlistF,jnrlistG,jnrlistH;
64     int              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
65     int              j_coord_offsetE,j_coord_offsetF,j_coord_offsetG,j_coord_offsetH;
66     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
67     real             rcutoff_scalar;
68     real             *shiftvec,*fshift,*x,*f;
69     real             *fjptrA,*fjptrB,*fjptrC,*fjptrD,*fjptrE,*fjptrF,*fjptrG,*fjptrH;
70     real             scratch[4*DIM];
71     __m256           tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
72     real *           vdwioffsetptr1;
73     __m256           ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
74     real *           vdwioffsetptr2;
75     __m256           ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
76     real *           vdwioffsetptr3;
77     __m256           ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
78     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D,vdwjidx0E,vdwjidx0F,vdwjidx0G,vdwjidx0H;
79     __m256           jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
80     __m256           dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
81     __m256           dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
82     __m256           dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
83     __m256           velec,felec,velecsum,facel,crf,krf,krf2;
84     real             *charge;
85     __m256           dummy_mask,cutoff_mask;
86     __m256           signbit = _mm256_castsi256_ps( _mm256_set1_epi32(0x80000000) );
87     __m256           one     = _mm256_set1_ps(1.0);
88     __m256           two     = _mm256_set1_ps(2.0);
89     x                = xx[0];
90     f                = ff[0];
91
92     nri              = nlist->nri;
93     iinr             = nlist->iinr;
94     jindex           = nlist->jindex;
95     jjnr             = nlist->jjnr;
96     shiftidx         = nlist->shift;
97     gid              = nlist->gid;
98     shiftvec         = fr->shift_vec[0];
99     fshift           = fr->fshift[0];
100     facel            = _mm256_set1_ps(fr->epsfac);
101     charge           = mdatoms->chargeA;
102     krf              = _mm256_set1_ps(fr->ic->k_rf);
103     krf2             = _mm256_set1_ps(fr->ic->k_rf*2.0);
104     crf              = _mm256_set1_ps(fr->ic->c_rf);
105
106     /* Setup water-specific parameters */
107     inr              = nlist->iinr[0];
108     iq1              = _mm256_mul_ps(facel,_mm256_set1_ps(charge[inr+1]));
109     iq2              = _mm256_mul_ps(facel,_mm256_set1_ps(charge[inr+2]));
110     iq3              = _mm256_mul_ps(facel,_mm256_set1_ps(charge[inr+3]));
111
112     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
113     rcutoff_scalar   = fr->rcoulomb;
114     rcutoff          = _mm256_set1_ps(rcutoff_scalar);
115     rcutoff2         = _mm256_mul_ps(rcutoff,rcutoff);
116
117     /* Avoid stupid compiler warnings */
118     jnrA = jnrB = jnrC = jnrD = jnrE = jnrF = jnrG = jnrH = 0;
119     j_coord_offsetA = 0;
120     j_coord_offsetB = 0;
121     j_coord_offsetC = 0;
122     j_coord_offsetD = 0;
123     j_coord_offsetE = 0;
124     j_coord_offsetF = 0;
125     j_coord_offsetG = 0;
126     j_coord_offsetH = 0;
127
128     outeriter        = 0;
129     inneriter        = 0;
130
131     for(iidx=0;iidx<4*DIM;iidx++)
132     {
133         scratch[iidx] = 0.0;
134     }
135
136     /* Start outer loop over neighborlists */
137     for(iidx=0; iidx<nri; iidx++)
138     {
139         /* Load shift vector for this list */
140         i_shift_offset   = DIM*shiftidx[iidx];
141
142         /* Load limits for loop over neighbors */
143         j_index_start    = jindex[iidx];
144         j_index_end      = jindex[iidx+1];
145
146         /* Get outer coordinate index */
147         inr              = iinr[iidx];
148         i_coord_offset   = DIM*inr;
149
150         /* Load i particle coords and add shift vector */
151         gmx_mm256_load_shift_and_3rvec_broadcast_ps(shiftvec+i_shift_offset,x+i_coord_offset+DIM,
152                                                     &ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
153
154         fix1             = _mm256_setzero_ps();
155         fiy1             = _mm256_setzero_ps();
156         fiz1             = _mm256_setzero_ps();
157         fix2             = _mm256_setzero_ps();
158         fiy2             = _mm256_setzero_ps();
159         fiz2             = _mm256_setzero_ps();
160         fix3             = _mm256_setzero_ps();
161         fiy3             = _mm256_setzero_ps();
162         fiz3             = _mm256_setzero_ps();
163
164         /* Reset potential sums */
165         velecsum         = _mm256_setzero_ps();
166
167         /* Start inner kernel loop */
168         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+7]>=0; jidx+=8)
169         {
170
171             /* Get j neighbor index, and coordinate index */
172             jnrA             = jjnr[jidx];
173             jnrB             = jjnr[jidx+1];
174             jnrC             = jjnr[jidx+2];
175             jnrD             = jjnr[jidx+3];
176             jnrE             = jjnr[jidx+4];
177             jnrF             = jjnr[jidx+5];
178             jnrG             = jjnr[jidx+6];
179             jnrH             = jjnr[jidx+7];
180             j_coord_offsetA  = DIM*jnrA;
181             j_coord_offsetB  = DIM*jnrB;
182             j_coord_offsetC  = DIM*jnrC;
183             j_coord_offsetD  = DIM*jnrD;
184             j_coord_offsetE  = DIM*jnrE;
185             j_coord_offsetF  = DIM*jnrF;
186             j_coord_offsetG  = DIM*jnrG;
187             j_coord_offsetH  = DIM*jnrH;
188
189             /* load j atom coordinates */
190             gmx_mm256_load_1rvec_8ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
191                                                  x+j_coord_offsetC,x+j_coord_offsetD,
192                                                  x+j_coord_offsetE,x+j_coord_offsetF,
193                                                  x+j_coord_offsetG,x+j_coord_offsetH,
194                                                  &jx0,&jy0,&jz0);
195
196             /* Calculate displacement vector */
197             dx10             = _mm256_sub_ps(ix1,jx0);
198             dy10             = _mm256_sub_ps(iy1,jy0);
199             dz10             = _mm256_sub_ps(iz1,jz0);
200             dx20             = _mm256_sub_ps(ix2,jx0);
201             dy20             = _mm256_sub_ps(iy2,jy0);
202             dz20             = _mm256_sub_ps(iz2,jz0);
203             dx30             = _mm256_sub_ps(ix3,jx0);
204             dy30             = _mm256_sub_ps(iy3,jy0);
205             dz30             = _mm256_sub_ps(iz3,jz0);
206
207             /* Calculate squared distance and things based on it */
208             rsq10            = gmx_mm256_calc_rsq_ps(dx10,dy10,dz10);
209             rsq20            = gmx_mm256_calc_rsq_ps(dx20,dy20,dz20);
210             rsq30            = gmx_mm256_calc_rsq_ps(dx30,dy30,dz30);
211
212             rinv10           = gmx_mm256_invsqrt_ps(rsq10);
213             rinv20           = gmx_mm256_invsqrt_ps(rsq20);
214             rinv30           = gmx_mm256_invsqrt_ps(rsq30);
215
216             rinvsq10         = _mm256_mul_ps(rinv10,rinv10);
217             rinvsq20         = _mm256_mul_ps(rinv20,rinv20);
218             rinvsq30         = _mm256_mul_ps(rinv30,rinv30);
219
220             /* Load parameters for j particles */
221             jq0              = gmx_mm256_load_8real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
222                                                                  charge+jnrC+0,charge+jnrD+0,
223                                                                  charge+jnrE+0,charge+jnrF+0,
224                                                                  charge+jnrG+0,charge+jnrH+0);
225
226             fjx0             = _mm256_setzero_ps();
227             fjy0             = _mm256_setzero_ps();
228             fjz0             = _mm256_setzero_ps();
229
230             /**************************
231              * CALCULATE INTERACTIONS *
232              **************************/
233
234             if (gmx_mm256_any_lt(rsq10,rcutoff2))
235             {
236
237             /* Compute parameters for interactions between i and j atoms */
238             qq10             = _mm256_mul_ps(iq1,jq0);
239
240             /* REACTION-FIELD ELECTROSTATICS */
241             velec            = _mm256_mul_ps(qq10,_mm256_sub_ps(_mm256_add_ps(rinv10,_mm256_mul_ps(krf,rsq10)),crf));
242             felec            = _mm256_mul_ps(qq10,_mm256_sub_ps(_mm256_mul_ps(rinv10,rinvsq10),krf2));
243
244             cutoff_mask      = _mm256_cmp_ps(rsq10,rcutoff2,_CMP_LT_OQ);
245
246             /* Update potential sum for this i atom from the interaction with this j atom. */
247             velec            = _mm256_and_ps(velec,cutoff_mask);
248             velecsum         = _mm256_add_ps(velecsum,velec);
249
250             fscal            = felec;
251
252             fscal            = _mm256_and_ps(fscal,cutoff_mask);
253
254             /* Calculate temporary vectorial force */
255             tx               = _mm256_mul_ps(fscal,dx10);
256             ty               = _mm256_mul_ps(fscal,dy10);
257             tz               = _mm256_mul_ps(fscal,dz10);
258
259             /* Update vectorial force */
260             fix1             = _mm256_add_ps(fix1,tx);
261             fiy1             = _mm256_add_ps(fiy1,ty);
262             fiz1             = _mm256_add_ps(fiz1,tz);
263
264             fjx0             = _mm256_add_ps(fjx0,tx);
265             fjy0             = _mm256_add_ps(fjy0,ty);
266             fjz0             = _mm256_add_ps(fjz0,tz);
267
268             }
269
270             /**************************
271              * CALCULATE INTERACTIONS *
272              **************************/
273
274             if (gmx_mm256_any_lt(rsq20,rcutoff2))
275             {
276
277             /* Compute parameters for interactions between i and j atoms */
278             qq20             = _mm256_mul_ps(iq2,jq0);
279
280             /* REACTION-FIELD ELECTROSTATICS */
281             velec            = _mm256_mul_ps(qq20,_mm256_sub_ps(_mm256_add_ps(rinv20,_mm256_mul_ps(krf,rsq20)),crf));
282             felec            = _mm256_mul_ps(qq20,_mm256_sub_ps(_mm256_mul_ps(rinv20,rinvsq20),krf2));
283
284             cutoff_mask      = _mm256_cmp_ps(rsq20,rcutoff2,_CMP_LT_OQ);
285
286             /* Update potential sum for this i atom from the interaction with this j atom. */
287             velec            = _mm256_and_ps(velec,cutoff_mask);
288             velecsum         = _mm256_add_ps(velecsum,velec);
289
290             fscal            = felec;
291
292             fscal            = _mm256_and_ps(fscal,cutoff_mask);
293
294             /* Calculate temporary vectorial force */
295             tx               = _mm256_mul_ps(fscal,dx20);
296             ty               = _mm256_mul_ps(fscal,dy20);
297             tz               = _mm256_mul_ps(fscal,dz20);
298
299             /* Update vectorial force */
300             fix2             = _mm256_add_ps(fix2,tx);
301             fiy2             = _mm256_add_ps(fiy2,ty);
302             fiz2             = _mm256_add_ps(fiz2,tz);
303
304             fjx0             = _mm256_add_ps(fjx0,tx);
305             fjy0             = _mm256_add_ps(fjy0,ty);
306             fjz0             = _mm256_add_ps(fjz0,tz);
307
308             }
309
310             /**************************
311              * CALCULATE INTERACTIONS *
312              **************************/
313
314             if (gmx_mm256_any_lt(rsq30,rcutoff2))
315             {
316
317             /* Compute parameters for interactions between i and j atoms */
318             qq30             = _mm256_mul_ps(iq3,jq0);
319
320             /* REACTION-FIELD ELECTROSTATICS */
321             velec            = _mm256_mul_ps(qq30,_mm256_sub_ps(_mm256_add_ps(rinv30,_mm256_mul_ps(krf,rsq30)),crf));
322             felec            = _mm256_mul_ps(qq30,_mm256_sub_ps(_mm256_mul_ps(rinv30,rinvsq30),krf2));
323
324             cutoff_mask      = _mm256_cmp_ps(rsq30,rcutoff2,_CMP_LT_OQ);
325
326             /* Update potential sum for this i atom from the interaction with this j atom. */
327             velec            = _mm256_and_ps(velec,cutoff_mask);
328             velecsum         = _mm256_add_ps(velecsum,velec);
329
330             fscal            = felec;
331
332             fscal            = _mm256_and_ps(fscal,cutoff_mask);
333
334             /* Calculate temporary vectorial force */
335             tx               = _mm256_mul_ps(fscal,dx30);
336             ty               = _mm256_mul_ps(fscal,dy30);
337             tz               = _mm256_mul_ps(fscal,dz30);
338
339             /* Update vectorial force */
340             fix3             = _mm256_add_ps(fix3,tx);
341             fiy3             = _mm256_add_ps(fiy3,ty);
342             fiz3             = _mm256_add_ps(fiz3,tz);
343
344             fjx0             = _mm256_add_ps(fjx0,tx);
345             fjy0             = _mm256_add_ps(fjy0,ty);
346             fjz0             = _mm256_add_ps(fjz0,tz);
347
348             }
349
350             fjptrA             = f+j_coord_offsetA;
351             fjptrB             = f+j_coord_offsetB;
352             fjptrC             = f+j_coord_offsetC;
353             fjptrD             = f+j_coord_offsetD;
354             fjptrE             = f+j_coord_offsetE;
355             fjptrF             = f+j_coord_offsetF;
356             fjptrG             = f+j_coord_offsetG;
357             fjptrH             = f+j_coord_offsetH;
358
359             gmx_mm256_decrement_1rvec_8ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,fjptrE,fjptrF,fjptrG,fjptrH,fjx0,fjy0,fjz0);
360
361             /* Inner loop uses 111 flops */
362         }
363
364         if(jidx<j_index_end)
365         {
366
367             /* Get j neighbor index, and coordinate index */
368             jnrlistA         = jjnr[jidx];
369             jnrlistB         = jjnr[jidx+1];
370             jnrlistC         = jjnr[jidx+2];
371             jnrlistD         = jjnr[jidx+3];
372             jnrlistE         = jjnr[jidx+4];
373             jnrlistF         = jjnr[jidx+5];
374             jnrlistG         = jjnr[jidx+6];
375             jnrlistH         = jjnr[jidx+7];
376             /* Sign of each element will be negative for non-real atoms.
377              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
378              * so use it as val = _mm_andnot_ps(mask,val) to clear dummy entries.
379              */
380             dummy_mask = gmx_mm256_set_m128(gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx+4)),_mm_setzero_si128())),
381                                             gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128())));
382                                             
383             jnrA       = (jnrlistA>=0) ? jnrlistA : 0;
384             jnrB       = (jnrlistB>=0) ? jnrlistB : 0;
385             jnrC       = (jnrlistC>=0) ? jnrlistC : 0;
386             jnrD       = (jnrlistD>=0) ? jnrlistD : 0;
387             jnrE       = (jnrlistE>=0) ? jnrlistE : 0;
388             jnrF       = (jnrlistF>=0) ? jnrlistF : 0;
389             jnrG       = (jnrlistG>=0) ? jnrlistG : 0;
390             jnrH       = (jnrlistH>=0) ? jnrlistH : 0;
391             j_coord_offsetA  = DIM*jnrA;
392             j_coord_offsetB  = DIM*jnrB;
393             j_coord_offsetC  = DIM*jnrC;
394             j_coord_offsetD  = DIM*jnrD;
395             j_coord_offsetE  = DIM*jnrE;
396             j_coord_offsetF  = DIM*jnrF;
397             j_coord_offsetG  = DIM*jnrG;
398             j_coord_offsetH  = DIM*jnrH;
399
400             /* load j atom coordinates */
401             gmx_mm256_load_1rvec_8ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
402                                                  x+j_coord_offsetC,x+j_coord_offsetD,
403                                                  x+j_coord_offsetE,x+j_coord_offsetF,
404                                                  x+j_coord_offsetG,x+j_coord_offsetH,
405                                                  &jx0,&jy0,&jz0);
406
407             /* Calculate displacement vector */
408             dx10             = _mm256_sub_ps(ix1,jx0);
409             dy10             = _mm256_sub_ps(iy1,jy0);
410             dz10             = _mm256_sub_ps(iz1,jz0);
411             dx20             = _mm256_sub_ps(ix2,jx0);
412             dy20             = _mm256_sub_ps(iy2,jy0);
413             dz20             = _mm256_sub_ps(iz2,jz0);
414             dx30             = _mm256_sub_ps(ix3,jx0);
415             dy30             = _mm256_sub_ps(iy3,jy0);
416             dz30             = _mm256_sub_ps(iz3,jz0);
417
418             /* Calculate squared distance and things based on it */
419             rsq10            = gmx_mm256_calc_rsq_ps(dx10,dy10,dz10);
420             rsq20            = gmx_mm256_calc_rsq_ps(dx20,dy20,dz20);
421             rsq30            = gmx_mm256_calc_rsq_ps(dx30,dy30,dz30);
422
423             rinv10           = gmx_mm256_invsqrt_ps(rsq10);
424             rinv20           = gmx_mm256_invsqrt_ps(rsq20);
425             rinv30           = gmx_mm256_invsqrt_ps(rsq30);
426
427             rinvsq10         = _mm256_mul_ps(rinv10,rinv10);
428             rinvsq20         = _mm256_mul_ps(rinv20,rinv20);
429             rinvsq30         = _mm256_mul_ps(rinv30,rinv30);
430
431             /* Load parameters for j particles */
432             jq0              = gmx_mm256_load_8real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
433                                                                  charge+jnrC+0,charge+jnrD+0,
434                                                                  charge+jnrE+0,charge+jnrF+0,
435                                                                  charge+jnrG+0,charge+jnrH+0);
436
437             fjx0             = _mm256_setzero_ps();
438             fjy0             = _mm256_setzero_ps();
439             fjz0             = _mm256_setzero_ps();
440
441             /**************************
442              * CALCULATE INTERACTIONS *
443              **************************/
444
445             if (gmx_mm256_any_lt(rsq10,rcutoff2))
446             {
447
448             /* Compute parameters for interactions between i and j atoms */
449             qq10             = _mm256_mul_ps(iq1,jq0);
450
451             /* REACTION-FIELD ELECTROSTATICS */
452             velec            = _mm256_mul_ps(qq10,_mm256_sub_ps(_mm256_add_ps(rinv10,_mm256_mul_ps(krf,rsq10)),crf));
453             felec            = _mm256_mul_ps(qq10,_mm256_sub_ps(_mm256_mul_ps(rinv10,rinvsq10),krf2));
454
455             cutoff_mask      = _mm256_cmp_ps(rsq10,rcutoff2,_CMP_LT_OQ);
456
457             /* Update potential sum for this i atom from the interaction with this j atom. */
458             velec            = _mm256_and_ps(velec,cutoff_mask);
459             velec            = _mm256_andnot_ps(dummy_mask,velec);
460             velecsum         = _mm256_add_ps(velecsum,velec);
461
462             fscal            = felec;
463
464             fscal            = _mm256_and_ps(fscal,cutoff_mask);
465
466             fscal            = _mm256_andnot_ps(dummy_mask,fscal);
467
468             /* Calculate temporary vectorial force */
469             tx               = _mm256_mul_ps(fscal,dx10);
470             ty               = _mm256_mul_ps(fscal,dy10);
471             tz               = _mm256_mul_ps(fscal,dz10);
472
473             /* Update vectorial force */
474             fix1             = _mm256_add_ps(fix1,tx);
475             fiy1             = _mm256_add_ps(fiy1,ty);
476             fiz1             = _mm256_add_ps(fiz1,tz);
477
478             fjx0             = _mm256_add_ps(fjx0,tx);
479             fjy0             = _mm256_add_ps(fjy0,ty);
480             fjz0             = _mm256_add_ps(fjz0,tz);
481
482             }
483
484             /**************************
485              * CALCULATE INTERACTIONS *
486              **************************/
487
488             if (gmx_mm256_any_lt(rsq20,rcutoff2))
489             {
490
491             /* Compute parameters for interactions between i and j atoms */
492             qq20             = _mm256_mul_ps(iq2,jq0);
493
494             /* REACTION-FIELD ELECTROSTATICS */
495             velec            = _mm256_mul_ps(qq20,_mm256_sub_ps(_mm256_add_ps(rinv20,_mm256_mul_ps(krf,rsq20)),crf));
496             felec            = _mm256_mul_ps(qq20,_mm256_sub_ps(_mm256_mul_ps(rinv20,rinvsq20),krf2));
497
498             cutoff_mask      = _mm256_cmp_ps(rsq20,rcutoff2,_CMP_LT_OQ);
499
500             /* Update potential sum for this i atom from the interaction with this j atom. */
501             velec            = _mm256_and_ps(velec,cutoff_mask);
502             velec            = _mm256_andnot_ps(dummy_mask,velec);
503             velecsum         = _mm256_add_ps(velecsum,velec);
504
505             fscal            = felec;
506
507             fscal            = _mm256_and_ps(fscal,cutoff_mask);
508
509             fscal            = _mm256_andnot_ps(dummy_mask,fscal);
510
511             /* Calculate temporary vectorial force */
512             tx               = _mm256_mul_ps(fscal,dx20);
513             ty               = _mm256_mul_ps(fscal,dy20);
514             tz               = _mm256_mul_ps(fscal,dz20);
515
516             /* Update vectorial force */
517             fix2             = _mm256_add_ps(fix2,tx);
518             fiy2             = _mm256_add_ps(fiy2,ty);
519             fiz2             = _mm256_add_ps(fiz2,tz);
520
521             fjx0             = _mm256_add_ps(fjx0,tx);
522             fjy0             = _mm256_add_ps(fjy0,ty);
523             fjz0             = _mm256_add_ps(fjz0,tz);
524
525             }
526
527             /**************************
528              * CALCULATE INTERACTIONS *
529              **************************/
530
531             if (gmx_mm256_any_lt(rsq30,rcutoff2))
532             {
533
534             /* Compute parameters for interactions between i and j atoms */
535             qq30             = _mm256_mul_ps(iq3,jq0);
536
537             /* REACTION-FIELD ELECTROSTATICS */
538             velec            = _mm256_mul_ps(qq30,_mm256_sub_ps(_mm256_add_ps(rinv30,_mm256_mul_ps(krf,rsq30)),crf));
539             felec            = _mm256_mul_ps(qq30,_mm256_sub_ps(_mm256_mul_ps(rinv30,rinvsq30),krf2));
540
541             cutoff_mask      = _mm256_cmp_ps(rsq30,rcutoff2,_CMP_LT_OQ);
542
543             /* Update potential sum for this i atom from the interaction with this j atom. */
544             velec            = _mm256_and_ps(velec,cutoff_mask);
545             velec            = _mm256_andnot_ps(dummy_mask,velec);
546             velecsum         = _mm256_add_ps(velecsum,velec);
547
548             fscal            = felec;
549
550             fscal            = _mm256_and_ps(fscal,cutoff_mask);
551
552             fscal            = _mm256_andnot_ps(dummy_mask,fscal);
553
554             /* Calculate temporary vectorial force */
555             tx               = _mm256_mul_ps(fscal,dx30);
556             ty               = _mm256_mul_ps(fscal,dy30);
557             tz               = _mm256_mul_ps(fscal,dz30);
558
559             /* Update vectorial force */
560             fix3             = _mm256_add_ps(fix3,tx);
561             fiy3             = _mm256_add_ps(fiy3,ty);
562             fiz3             = _mm256_add_ps(fiz3,tz);
563
564             fjx0             = _mm256_add_ps(fjx0,tx);
565             fjy0             = _mm256_add_ps(fjy0,ty);
566             fjz0             = _mm256_add_ps(fjz0,tz);
567
568             }
569
570             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
571             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
572             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
573             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
574             fjptrE             = (jnrlistE>=0) ? f+j_coord_offsetE : scratch;
575             fjptrF             = (jnrlistF>=0) ? f+j_coord_offsetF : scratch;
576             fjptrG             = (jnrlistG>=0) ? f+j_coord_offsetG : scratch;
577             fjptrH             = (jnrlistH>=0) ? f+j_coord_offsetH : scratch;
578
579             gmx_mm256_decrement_1rvec_8ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,fjptrE,fjptrF,fjptrG,fjptrH,fjx0,fjy0,fjz0);
580
581             /* Inner loop uses 111 flops */
582         }
583
584         /* End of innermost loop */
585
586         gmx_mm256_update_iforce_3atom_swizzle_ps(fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
587                                                  f+i_coord_offset+DIM,fshift+i_shift_offset);
588
589         ggid                        = gid[iidx];
590         /* Update potential energies */
591         gmx_mm256_update_1pot_ps(velecsum,kernel_data->energygrp_elec+ggid);
592
593         /* Increment number of inner iterations */
594         inneriter                  += j_index_end - j_index_start;
595
596         /* Outer loop uses 19 flops */
597     }
598
599     /* Increment number of outer iterations */
600     outeriter        += nri;
601
602     /* Update outer/inner flops */
603
604     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_W4_VF,outeriter*19 + inneriter*111);
605 }
606 /*
607  * Gromacs nonbonded kernel:   nb_kernel_ElecRFCut_VdwNone_GeomW4P1_F_avx_256_single
608  * Electrostatics interaction: ReactionField
609  * VdW interaction:            None
610  * Geometry:                   Water4-Particle
611  * Calculate force/pot:        Force
612  */
613 void
614 nb_kernel_ElecRFCut_VdwNone_GeomW4P1_F_avx_256_single
615                     (t_nblist * gmx_restrict                nlist,
616                      rvec * gmx_restrict                    xx,
617                      rvec * gmx_restrict                    ff,
618                      t_forcerec * gmx_restrict              fr,
619                      t_mdatoms * gmx_restrict               mdatoms,
620                      nb_kernel_data_t * gmx_restrict        kernel_data,
621                      t_nrnb * gmx_restrict                  nrnb)
622 {
623     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or 
624      * just 0 for non-waters.
625      * Suffixes A,B,C,D,E,F,G,H refer to j loop unrolling done with AVX, e.g. for the eight different
626      * jnr indices corresponding to data put in the four positions in the SIMD register.
627      */
628     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
629     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
630     int              jnrA,jnrB,jnrC,jnrD;
631     int              jnrE,jnrF,jnrG,jnrH;
632     int              jnrlistA,jnrlistB,jnrlistC,jnrlistD;
633     int              jnrlistE,jnrlistF,jnrlistG,jnrlistH;
634     int              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
635     int              j_coord_offsetE,j_coord_offsetF,j_coord_offsetG,j_coord_offsetH;
636     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
637     real             rcutoff_scalar;
638     real             *shiftvec,*fshift,*x,*f;
639     real             *fjptrA,*fjptrB,*fjptrC,*fjptrD,*fjptrE,*fjptrF,*fjptrG,*fjptrH;
640     real             scratch[4*DIM];
641     __m256           tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
642     real *           vdwioffsetptr1;
643     __m256           ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
644     real *           vdwioffsetptr2;
645     __m256           ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
646     real *           vdwioffsetptr3;
647     __m256           ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
648     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D,vdwjidx0E,vdwjidx0F,vdwjidx0G,vdwjidx0H;
649     __m256           jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
650     __m256           dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
651     __m256           dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
652     __m256           dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
653     __m256           velec,felec,velecsum,facel,crf,krf,krf2;
654     real             *charge;
655     __m256           dummy_mask,cutoff_mask;
656     __m256           signbit = _mm256_castsi256_ps( _mm256_set1_epi32(0x80000000) );
657     __m256           one     = _mm256_set1_ps(1.0);
658     __m256           two     = _mm256_set1_ps(2.0);
659     x                = xx[0];
660     f                = ff[0];
661
662     nri              = nlist->nri;
663     iinr             = nlist->iinr;
664     jindex           = nlist->jindex;
665     jjnr             = nlist->jjnr;
666     shiftidx         = nlist->shift;
667     gid              = nlist->gid;
668     shiftvec         = fr->shift_vec[0];
669     fshift           = fr->fshift[0];
670     facel            = _mm256_set1_ps(fr->epsfac);
671     charge           = mdatoms->chargeA;
672     krf              = _mm256_set1_ps(fr->ic->k_rf);
673     krf2             = _mm256_set1_ps(fr->ic->k_rf*2.0);
674     crf              = _mm256_set1_ps(fr->ic->c_rf);
675
676     /* Setup water-specific parameters */
677     inr              = nlist->iinr[0];
678     iq1              = _mm256_mul_ps(facel,_mm256_set1_ps(charge[inr+1]));
679     iq2              = _mm256_mul_ps(facel,_mm256_set1_ps(charge[inr+2]));
680     iq3              = _mm256_mul_ps(facel,_mm256_set1_ps(charge[inr+3]));
681
682     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
683     rcutoff_scalar   = fr->rcoulomb;
684     rcutoff          = _mm256_set1_ps(rcutoff_scalar);
685     rcutoff2         = _mm256_mul_ps(rcutoff,rcutoff);
686
687     /* Avoid stupid compiler warnings */
688     jnrA = jnrB = jnrC = jnrD = jnrE = jnrF = jnrG = jnrH = 0;
689     j_coord_offsetA = 0;
690     j_coord_offsetB = 0;
691     j_coord_offsetC = 0;
692     j_coord_offsetD = 0;
693     j_coord_offsetE = 0;
694     j_coord_offsetF = 0;
695     j_coord_offsetG = 0;
696     j_coord_offsetH = 0;
697
698     outeriter        = 0;
699     inneriter        = 0;
700
701     for(iidx=0;iidx<4*DIM;iidx++)
702     {
703         scratch[iidx] = 0.0;
704     }
705
706     /* Start outer loop over neighborlists */
707     for(iidx=0; iidx<nri; iidx++)
708     {
709         /* Load shift vector for this list */
710         i_shift_offset   = DIM*shiftidx[iidx];
711
712         /* Load limits for loop over neighbors */
713         j_index_start    = jindex[iidx];
714         j_index_end      = jindex[iidx+1];
715
716         /* Get outer coordinate index */
717         inr              = iinr[iidx];
718         i_coord_offset   = DIM*inr;
719
720         /* Load i particle coords and add shift vector */
721         gmx_mm256_load_shift_and_3rvec_broadcast_ps(shiftvec+i_shift_offset,x+i_coord_offset+DIM,
722                                                     &ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
723
724         fix1             = _mm256_setzero_ps();
725         fiy1             = _mm256_setzero_ps();
726         fiz1             = _mm256_setzero_ps();
727         fix2             = _mm256_setzero_ps();
728         fiy2             = _mm256_setzero_ps();
729         fiz2             = _mm256_setzero_ps();
730         fix3             = _mm256_setzero_ps();
731         fiy3             = _mm256_setzero_ps();
732         fiz3             = _mm256_setzero_ps();
733
734         /* Start inner kernel loop */
735         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+7]>=0; jidx+=8)
736         {
737
738             /* Get j neighbor index, and coordinate index */
739             jnrA             = jjnr[jidx];
740             jnrB             = jjnr[jidx+1];
741             jnrC             = jjnr[jidx+2];
742             jnrD             = jjnr[jidx+3];
743             jnrE             = jjnr[jidx+4];
744             jnrF             = jjnr[jidx+5];
745             jnrG             = jjnr[jidx+6];
746             jnrH             = jjnr[jidx+7];
747             j_coord_offsetA  = DIM*jnrA;
748             j_coord_offsetB  = DIM*jnrB;
749             j_coord_offsetC  = DIM*jnrC;
750             j_coord_offsetD  = DIM*jnrD;
751             j_coord_offsetE  = DIM*jnrE;
752             j_coord_offsetF  = DIM*jnrF;
753             j_coord_offsetG  = DIM*jnrG;
754             j_coord_offsetH  = DIM*jnrH;
755
756             /* load j atom coordinates */
757             gmx_mm256_load_1rvec_8ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
758                                                  x+j_coord_offsetC,x+j_coord_offsetD,
759                                                  x+j_coord_offsetE,x+j_coord_offsetF,
760                                                  x+j_coord_offsetG,x+j_coord_offsetH,
761                                                  &jx0,&jy0,&jz0);
762
763             /* Calculate displacement vector */
764             dx10             = _mm256_sub_ps(ix1,jx0);
765             dy10             = _mm256_sub_ps(iy1,jy0);
766             dz10             = _mm256_sub_ps(iz1,jz0);
767             dx20             = _mm256_sub_ps(ix2,jx0);
768             dy20             = _mm256_sub_ps(iy2,jy0);
769             dz20             = _mm256_sub_ps(iz2,jz0);
770             dx30             = _mm256_sub_ps(ix3,jx0);
771             dy30             = _mm256_sub_ps(iy3,jy0);
772             dz30             = _mm256_sub_ps(iz3,jz0);
773
774             /* Calculate squared distance and things based on it */
775             rsq10            = gmx_mm256_calc_rsq_ps(dx10,dy10,dz10);
776             rsq20            = gmx_mm256_calc_rsq_ps(dx20,dy20,dz20);
777             rsq30            = gmx_mm256_calc_rsq_ps(dx30,dy30,dz30);
778
779             rinv10           = gmx_mm256_invsqrt_ps(rsq10);
780             rinv20           = gmx_mm256_invsqrt_ps(rsq20);
781             rinv30           = gmx_mm256_invsqrt_ps(rsq30);
782
783             rinvsq10         = _mm256_mul_ps(rinv10,rinv10);
784             rinvsq20         = _mm256_mul_ps(rinv20,rinv20);
785             rinvsq30         = _mm256_mul_ps(rinv30,rinv30);
786
787             /* Load parameters for j particles */
788             jq0              = gmx_mm256_load_8real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
789                                                                  charge+jnrC+0,charge+jnrD+0,
790                                                                  charge+jnrE+0,charge+jnrF+0,
791                                                                  charge+jnrG+0,charge+jnrH+0);
792
793             fjx0             = _mm256_setzero_ps();
794             fjy0             = _mm256_setzero_ps();
795             fjz0             = _mm256_setzero_ps();
796
797             /**************************
798              * CALCULATE INTERACTIONS *
799              **************************/
800
801             if (gmx_mm256_any_lt(rsq10,rcutoff2))
802             {
803
804             /* Compute parameters for interactions between i and j atoms */
805             qq10             = _mm256_mul_ps(iq1,jq0);
806
807             /* REACTION-FIELD ELECTROSTATICS */
808             felec            = _mm256_mul_ps(qq10,_mm256_sub_ps(_mm256_mul_ps(rinv10,rinvsq10),krf2));
809
810             cutoff_mask      = _mm256_cmp_ps(rsq10,rcutoff2,_CMP_LT_OQ);
811
812             fscal            = felec;
813
814             fscal            = _mm256_and_ps(fscal,cutoff_mask);
815
816             /* Calculate temporary vectorial force */
817             tx               = _mm256_mul_ps(fscal,dx10);
818             ty               = _mm256_mul_ps(fscal,dy10);
819             tz               = _mm256_mul_ps(fscal,dz10);
820
821             /* Update vectorial force */
822             fix1             = _mm256_add_ps(fix1,tx);
823             fiy1             = _mm256_add_ps(fiy1,ty);
824             fiz1             = _mm256_add_ps(fiz1,tz);
825
826             fjx0             = _mm256_add_ps(fjx0,tx);
827             fjy0             = _mm256_add_ps(fjy0,ty);
828             fjz0             = _mm256_add_ps(fjz0,tz);
829
830             }
831
832             /**************************
833              * CALCULATE INTERACTIONS *
834              **************************/
835
836             if (gmx_mm256_any_lt(rsq20,rcutoff2))
837             {
838
839             /* Compute parameters for interactions between i and j atoms */
840             qq20             = _mm256_mul_ps(iq2,jq0);
841
842             /* REACTION-FIELD ELECTROSTATICS */
843             felec            = _mm256_mul_ps(qq20,_mm256_sub_ps(_mm256_mul_ps(rinv20,rinvsq20),krf2));
844
845             cutoff_mask      = _mm256_cmp_ps(rsq20,rcutoff2,_CMP_LT_OQ);
846
847             fscal            = felec;
848
849             fscal            = _mm256_and_ps(fscal,cutoff_mask);
850
851             /* Calculate temporary vectorial force */
852             tx               = _mm256_mul_ps(fscal,dx20);
853             ty               = _mm256_mul_ps(fscal,dy20);
854             tz               = _mm256_mul_ps(fscal,dz20);
855
856             /* Update vectorial force */
857             fix2             = _mm256_add_ps(fix2,tx);
858             fiy2             = _mm256_add_ps(fiy2,ty);
859             fiz2             = _mm256_add_ps(fiz2,tz);
860
861             fjx0             = _mm256_add_ps(fjx0,tx);
862             fjy0             = _mm256_add_ps(fjy0,ty);
863             fjz0             = _mm256_add_ps(fjz0,tz);
864
865             }
866
867             /**************************
868              * CALCULATE INTERACTIONS *
869              **************************/
870
871             if (gmx_mm256_any_lt(rsq30,rcutoff2))
872             {
873
874             /* Compute parameters for interactions between i and j atoms */
875             qq30             = _mm256_mul_ps(iq3,jq0);
876
877             /* REACTION-FIELD ELECTROSTATICS */
878             felec            = _mm256_mul_ps(qq30,_mm256_sub_ps(_mm256_mul_ps(rinv30,rinvsq30),krf2));
879
880             cutoff_mask      = _mm256_cmp_ps(rsq30,rcutoff2,_CMP_LT_OQ);
881
882             fscal            = felec;
883
884             fscal            = _mm256_and_ps(fscal,cutoff_mask);
885
886             /* Calculate temporary vectorial force */
887             tx               = _mm256_mul_ps(fscal,dx30);
888             ty               = _mm256_mul_ps(fscal,dy30);
889             tz               = _mm256_mul_ps(fscal,dz30);
890
891             /* Update vectorial force */
892             fix3             = _mm256_add_ps(fix3,tx);
893             fiy3             = _mm256_add_ps(fiy3,ty);
894             fiz3             = _mm256_add_ps(fiz3,tz);
895
896             fjx0             = _mm256_add_ps(fjx0,tx);
897             fjy0             = _mm256_add_ps(fjy0,ty);
898             fjz0             = _mm256_add_ps(fjz0,tz);
899
900             }
901
902             fjptrA             = f+j_coord_offsetA;
903             fjptrB             = f+j_coord_offsetB;
904             fjptrC             = f+j_coord_offsetC;
905             fjptrD             = f+j_coord_offsetD;
906             fjptrE             = f+j_coord_offsetE;
907             fjptrF             = f+j_coord_offsetF;
908             fjptrG             = f+j_coord_offsetG;
909             fjptrH             = f+j_coord_offsetH;
910
911             gmx_mm256_decrement_1rvec_8ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,fjptrE,fjptrF,fjptrG,fjptrH,fjx0,fjy0,fjz0);
912
913             /* Inner loop uses 93 flops */
914         }
915
916         if(jidx<j_index_end)
917         {
918
919             /* Get j neighbor index, and coordinate index */
920             jnrlistA         = jjnr[jidx];
921             jnrlistB         = jjnr[jidx+1];
922             jnrlistC         = jjnr[jidx+2];
923             jnrlistD         = jjnr[jidx+3];
924             jnrlistE         = jjnr[jidx+4];
925             jnrlistF         = jjnr[jidx+5];
926             jnrlistG         = jjnr[jidx+6];
927             jnrlistH         = jjnr[jidx+7];
928             /* Sign of each element will be negative for non-real atoms.
929              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
930              * so use it as val = _mm_andnot_ps(mask,val) to clear dummy entries.
931              */
932             dummy_mask = gmx_mm256_set_m128(gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx+4)),_mm_setzero_si128())),
933                                             gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128())));
934                                             
935             jnrA       = (jnrlistA>=0) ? jnrlistA : 0;
936             jnrB       = (jnrlistB>=0) ? jnrlistB : 0;
937             jnrC       = (jnrlistC>=0) ? jnrlistC : 0;
938             jnrD       = (jnrlistD>=0) ? jnrlistD : 0;
939             jnrE       = (jnrlistE>=0) ? jnrlistE : 0;
940             jnrF       = (jnrlistF>=0) ? jnrlistF : 0;
941             jnrG       = (jnrlistG>=0) ? jnrlistG : 0;
942             jnrH       = (jnrlistH>=0) ? jnrlistH : 0;
943             j_coord_offsetA  = DIM*jnrA;
944             j_coord_offsetB  = DIM*jnrB;
945             j_coord_offsetC  = DIM*jnrC;
946             j_coord_offsetD  = DIM*jnrD;
947             j_coord_offsetE  = DIM*jnrE;
948             j_coord_offsetF  = DIM*jnrF;
949             j_coord_offsetG  = DIM*jnrG;
950             j_coord_offsetH  = DIM*jnrH;
951
952             /* load j atom coordinates */
953             gmx_mm256_load_1rvec_8ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
954                                                  x+j_coord_offsetC,x+j_coord_offsetD,
955                                                  x+j_coord_offsetE,x+j_coord_offsetF,
956                                                  x+j_coord_offsetG,x+j_coord_offsetH,
957                                                  &jx0,&jy0,&jz0);
958
959             /* Calculate displacement vector */
960             dx10             = _mm256_sub_ps(ix1,jx0);
961             dy10             = _mm256_sub_ps(iy1,jy0);
962             dz10             = _mm256_sub_ps(iz1,jz0);
963             dx20             = _mm256_sub_ps(ix2,jx0);
964             dy20             = _mm256_sub_ps(iy2,jy0);
965             dz20             = _mm256_sub_ps(iz2,jz0);
966             dx30             = _mm256_sub_ps(ix3,jx0);
967             dy30             = _mm256_sub_ps(iy3,jy0);
968             dz30             = _mm256_sub_ps(iz3,jz0);
969
970             /* Calculate squared distance and things based on it */
971             rsq10            = gmx_mm256_calc_rsq_ps(dx10,dy10,dz10);
972             rsq20            = gmx_mm256_calc_rsq_ps(dx20,dy20,dz20);
973             rsq30            = gmx_mm256_calc_rsq_ps(dx30,dy30,dz30);
974
975             rinv10           = gmx_mm256_invsqrt_ps(rsq10);
976             rinv20           = gmx_mm256_invsqrt_ps(rsq20);
977             rinv30           = gmx_mm256_invsqrt_ps(rsq30);
978
979             rinvsq10         = _mm256_mul_ps(rinv10,rinv10);
980             rinvsq20         = _mm256_mul_ps(rinv20,rinv20);
981             rinvsq30         = _mm256_mul_ps(rinv30,rinv30);
982
983             /* Load parameters for j particles */
984             jq0              = gmx_mm256_load_8real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
985                                                                  charge+jnrC+0,charge+jnrD+0,
986                                                                  charge+jnrE+0,charge+jnrF+0,
987                                                                  charge+jnrG+0,charge+jnrH+0);
988
989             fjx0             = _mm256_setzero_ps();
990             fjy0             = _mm256_setzero_ps();
991             fjz0             = _mm256_setzero_ps();
992
993             /**************************
994              * CALCULATE INTERACTIONS *
995              **************************/
996
997             if (gmx_mm256_any_lt(rsq10,rcutoff2))
998             {
999
1000             /* Compute parameters for interactions between i and j atoms */
1001             qq10             = _mm256_mul_ps(iq1,jq0);
1002
1003             /* REACTION-FIELD ELECTROSTATICS */
1004             felec            = _mm256_mul_ps(qq10,_mm256_sub_ps(_mm256_mul_ps(rinv10,rinvsq10),krf2));
1005
1006             cutoff_mask      = _mm256_cmp_ps(rsq10,rcutoff2,_CMP_LT_OQ);
1007
1008             fscal            = felec;
1009
1010             fscal            = _mm256_and_ps(fscal,cutoff_mask);
1011
1012             fscal            = _mm256_andnot_ps(dummy_mask,fscal);
1013
1014             /* Calculate temporary vectorial force */
1015             tx               = _mm256_mul_ps(fscal,dx10);
1016             ty               = _mm256_mul_ps(fscal,dy10);
1017             tz               = _mm256_mul_ps(fscal,dz10);
1018
1019             /* Update vectorial force */
1020             fix1             = _mm256_add_ps(fix1,tx);
1021             fiy1             = _mm256_add_ps(fiy1,ty);
1022             fiz1             = _mm256_add_ps(fiz1,tz);
1023
1024             fjx0             = _mm256_add_ps(fjx0,tx);
1025             fjy0             = _mm256_add_ps(fjy0,ty);
1026             fjz0             = _mm256_add_ps(fjz0,tz);
1027
1028             }
1029
1030             /**************************
1031              * CALCULATE INTERACTIONS *
1032              **************************/
1033
1034             if (gmx_mm256_any_lt(rsq20,rcutoff2))
1035             {
1036
1037             /* Compute parameters for interactions between i and j atoms */
1038             qq20             = _mm256_mul_ps(iq2,jq0);
1039
1040             /* REACTION-FIELD ELECTROSTATICS */
1041             felec            = _mm256_mul_ps(qq20,_mm256_sub_ps(_mm256_mul_ps(rinv20,rinvsq20),krf2));
1042
1043             cutoff_mask      = _mm256_cmp_ps(rsq20,rcutoff2,_CMP_LT_OQ);
1044
1045             fscal            = felec;
1046
1047             fscal            = _mm256_and_ps(fscal,cutoff_mask);
1048
1049             fscal            = _mm256_andnot_ps(dummy_mask,fscal);
1050
1051             /* Calculate temporary vectorial force */
1052             tx               = _mm256_mul_ps(fscal,dx20);
1053             ty               = _mm256_mul_ps(fscal,dy20);
1054             tz               = _mm256_mul_ps(fscal,dz20);
1055
1056             /* Update vectorial force */
1057             fix2             = _mm256_add_ps(fix2,tx);
1058             fiy2             = _mm256_add_ps(fiy2,ty);
1059             fiz2             = _mm256_add_ps(fiz2,tz);
1060
1061             fjx0             = _mm256_add_ps(fjx0,tx);
1062             fjy0             = _mm256_add_ps(fjy0,ty);
1063             fjz0             = _mm256_add_ps(fjz0,tz);
1064
1065             }
1066
1067             /**************************
1068              * CALCULATE INTERACTIONS *
1069              **************************/
1070
1071             if (gmx_mm256_any_lt(rsq30,rcutoff2))
1072             {
1073
1074             /* Compute parameters for interactions between i and j atoms */
1075             qq30             = _mm256_mul_ps(iq3,jq0);
1076
1077             /* REACTION-FIELD ELECTROSTATICS */
1078             felec            = _mm256_mul_ps(qq30,_mm256_sub_ps(_mm256_mul_ps(rinv30,rinvsq30),krf2));
1079
1080             cutoff_mask      = _mm256_cmp_ps(rsq30,rcutoff2,_CMP_LT_OQ);
1081
1082             fscal            = felec;
1083
1084             fscal            = _mm256_and_ps(fscal,cutoff_mask);
1085
1086             fscal            = _mm256_andnot_ps(dummy_mask,fscal);
1087
1088             /* Calculate temporary vectorial force */
1089             tx               = _mm256_mul_ps(fscal,dx30);
1090             ty               = _mm256_mul_ps(fscal,dy30);
1091             tz               = _mm256_mul_ps(fscal,dz30);
1092
1093             /* Update vectorial force */
1094             fix3             = _mm256_add_ps(fix3,tx);
1095             fiy3             = _mm256_add_ps(fiy3,ty);
1096             fiz3             = _mm256_add_ps(fiz3,tz);
1097
1098             fjx0             = _mm256_add_ps(fjx0,tx);
1099             fjy0             = _mm256_add_ps(fjy0,ty);
1100             fjz0             = _mm256_add_ps(fjz0,tz);
1101
1102             }
1103
1104             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
1105             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
1106             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
1107             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
1108             fjptrE             = (jnrlistE>=0) ? f+j_coord_offsetE : scratch;
1109             fjptrF             = (jnrlistF>=0) ? f+j_coord_offsetF : scratch;
1110             fjptrG             = (jnrlistG>=0) ? f+j_coord_offsetG : scratch;
1111             fjptrH             = (jnrlistH>=0) ? f+j_coord_offsetH : scratch;
1112
1113             gmx_mm256_decrement_1rvec_8ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,fjptrE,fjptrF,fjptrG,fjptrH,fjx0,fjy0,fjz0);
1114
1115             /* Inner loop uses 93 flops */
1116         }
1117
1118         /* End of innermost loop */
1119
1120         gmx_mm256_update_iforce_3atom_swizzle_ps(fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
1121                                                  f+i_coord_offset+DIM,fshift+i_shift_offset);
1122
1123         /* Increment number of inner iterations */
1124         inneriter                  += j_index_end - j_index_start;
1125
1126         /* Outer loop uses 18 flops */
1127     }
1128
1129     /* Increment number of outer iterations */
1130     outeriter        += nri;
1131
1132     /* Update outer/inner flops */
1133
1134     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_W4_F,outeriter*18 + inneriter*93);
1135 }