Merge release-4-6 into master
[alexxy/gromacs.git] / src / gromacs / gmxlib / nonbonded / nb_kernel_sse2_single / nb_kernel_ElecCSTab_VdwCSTab_GeomW3P1_sse2_single.c
1 /*
2  * Note: this file was generated by the Gromacs sse2_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_sse2_single.h"
34 #include "kernelutil_x86_sse2_single.h"
35
36 /*
37  * Gromacs nonbonded kernel:   nb_kernel_ElecCSTab_VdwCSTab_GeomW3P1_VF_sse2_single
38  * Electrostatics interaction: CubicSplineTable
39  * VdW interaction:            CubicSplineTable
40  * Geometry:                   Water3-Particle
41  * Calculate force/pot:        PotentialAndForce
42  */
43 void
44 nb_kernel_ElecCSTab_VdwCSTab_GeomW3P1_VF_sse2_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 refer to j loop unrolling done with SSE, e.g. for the four 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              jnrlistA,jnrlistB,jnrlistC,jnrlistD;
62     int              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
63     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
64     real             rcutoff_scalar;
65     real             *shiftvec,*fshift,*x,*f;
66     real             *fjptrA,*fjptrB,*fjptrC,*fjptrD;
67     real             scratch[4*DIM];
68     __m128           tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
69     int              vdwioffset0;
70     __m128           ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
71     int              vdwioffset1;
72     __m128           ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
73     int              vdwioffset2;
74     __m128           ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
75     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
76     __m128           jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
77     __m128           dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
78     __m128           dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
79     __m128           dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
80     __m128           velec,felec,velecsum,facel,crf,krf,krf2;
81     real             *charge;
82     int              nvdwtype;
83     __m128           rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
84     int              *vdwtype;
85     real             *vdwparam;
86     __m128           one_sixth   = _mm_set1_ps(1.0/6.0);
87     __m128           one_twelfth = _mm_set1_ps(1.0/12.0);
88     __m128i          vfitab;
89     __m128i          ifour       = _mm_set1_epi32(4);
90     __m128           rt,vfeps,vftabscale,Y,F,G,H,Heps,Fp,VV,FF;
91     real             *vftab;
92     __m128           dummy_mask,cutoff_mask;
93     __m128           signbit = _mm_castsi128_ps( _mm_set1_epi32(0x80000000) );
94     __m128           one     = _mm_set1_ps(1.0);
95     __m128           two     = _mm_set1_ps(2.0);
96     x                = xx[0];
97     f                = ff[0];
98
99     nri              = nlist->nri;
100     iinr             = nlist->iinr;
101     jindex           = nlist->jindex;
102     jjnr             = nlist->jjnr;
103     shiftidx         = nlist->shift;
104     gid              = nlist->gid;
105     shiftvec         = fr->shift_vec[0];
106     fshift           = fr->fshift[0];
107     facel            = _mm_set1_ps(fr->epsfac);
108     charge           = mdatoms->chargeA;
109     nvdwtype         = fr->ntype;
110     vdwparam         = fr->nbfp;
111     vdwtype          = mdatoms->typeA;
112
113     vftab            = kernel_data->table_elec_vdw->data;
114     vftabscale       = _mm_set1_ps(kernel_data->table_elec_vdw->scale);
115
116     /* Setup water-specific parameters */
117     inr              = nlist->iinr[0];
118     iq0              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+0]));
119     iq1              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+1]));
120     iq2              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+2]));
121     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
122
123     /* Avoid stupid compiler warnings */
124     jnrA = jnrB = jnrC = jnrD = 0;
125     j_coord_offsetA = 0;
126     j_coord_offsetB = 0;
127     j_coord_offsetC = 0;
128     j_coord_offsetD = 0;
129
130     outeriter        = 0;
131     inneriter        = 0;
132
133     for(iidx=0;iidx<4*DIM;iidx++)
134     {
135         scratch[iidx] = 0.0;
136     }  
137
138     /* Start outer loop over neighborlists */
139     for(iidx=0; iidx<nri; iidx++)
140     {
141         /* Load shift vector for this list */
142         i_shift_offset   = DIM*shiftidx[iidx];
143
144         /* Load limits for loop over neighbors */
145         j_index_start    = jindex[iidx];
146         j_index_end      = jindex[iidx+1];
147
148         /* Get outer coordinate index */
149         inr              = iinr[iidx];
150         i_coord_offset   = DIM*inr;
151
152         /* Load i particle coords and add shift vector */
153         gmx_mm_load_shift_and_3rvec_broadcast_ps(shiftvec+i_shift_offset,x+i_coord_offset,
154                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2);
155         
156         fix0             = _mm_setzero_ps();
157         fiy0             = _mm_setzero_ps();
158         fiz0             = _mm_setzero_ps();
159         fix1             = _mm_setzero_ps();
160         fiy1             = _mm_setzero_ps();
161         fiz1             = _mm_setzero_ps();
162         fix2             = _mm_setzero_ps();
163         fiy2             = _mm_setzero_ps();
164         fiz2             = _mm_setzero_ps();
165
166         /* Reset potential sums */
167         velecsum         = _mm_setzero_ps();
168         vvdwsum          = _mm_setzero_ps();
169
170         /* Start inner kernel loop */
171         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
172         {
173
174             /* Get j neighbor index, and coordinate index */
175             jnrA             = jjnr[jidx];
176             jnrB             = jjnr[jidx+1];
177             jnrC             = jjnr[jidx+2];
178             jnrD             = jjnr[jidx+3];
179             j_coord_offsetA  = DIM*jnrA;
180             j_coord_offsetB  = DIM*jnrB;
181             j_coord_offsetC  = DIM*jnrC;
182             j_coord_offsetD  = DIM*jnrD;
183
184             /* load j atom coordinates */
185             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
186                                               x+j_coord_offsetC,x+j_coord_offsetD,
187                                               &jx0,&jy0,&jz0);
188
189             /* Calculate displacement vector */
190             dx00             = _mm_sub_ps(ix0,jx0);
191             dy00             = _mm_sub_ps(iy0,jy0);
192             dz00             = _mm_sub_ps(iz0,jz0);
193             dx10             = _mm_sub_ps(ix1,jx0);
194             dy10             = _mm_sub_ps(iy1,jy0);
195             dz10             = _mm_sub_ps(iz1,jz0);
196             dx20             = _mm_sub_ps(ix2,jx0);
197             dy20             = _mm_sub_ps(iy2,jy0);
198             dz20             = _mm_sub_ps(iz2,jz0);
199
200             /* Calculate squared distance and things based on it */
201             rsq00            = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
202             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
203             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
204
205             rinv00           = gmx_mm_invsqrt_ps(rsq00);
206             rinv10           = gmx_mm_invsqrt_ps(rsq10);
207             rinv20           = gmx_mm_invsqrt_ps(rsq20);
208
209             /* Load parameters for j particles */
210             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
211                                                               charge+jnrC+0,charge+jnrD+0);
212             vdwjidx0A        = 2*vdwtype[jnrA+0];
213             vdwjidx0B        = 2*vdwtype[jnrB+0];
214             vdwjidx0C        = 2*vdwtype[jnrC+0];
215             vdwjidx0D        = 2*vdwtype[jnrD+0];
216
217             /**************************
218              * CALCULATE INTERACTIONS *
219              **************************/
220
221             r00              = _mm_mul_ps(rsq00,rinv00);
222
223             /* Compute parameters for interactions between i and j atoms */
224             qq00             = _mm_mul_ps(iq0,jq0);
225             gmx_mm_load_4pair_swizzle_ps(vdwparam+vdwioffset0+vdwjidx0A,
226                                          vdwparam+vdwioffset0+vdwjidx0B,
227                                          vdwparam+vdwioffset0+vdwjidx0C,
228                                          vdwparam+vdwioffset0+vdwjidx0D,
229                                          &c6_00,&c12_00);
230
231             /* Calculate table index by multiplying r with table scale and truncate to integer */
232             rt               = _mm_mul_ps(r00,vftabscale);
233             vfitab           = _mm_cvttps_epi32(rt);
234             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
235             vfitab           = _mm_slli_epi32(_mm_add_epi32(vfitab,_mm_slli_epi32(vfitab,1)),2);
236
237             /* CUBIC SPLINE TABLE ELECTROSTATICS */
238             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
239             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
240             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
241             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
242             _MM_TRANSPOSE4_PS(Y,F,G,H);
243             Heps             = _mm_mul_ps(vfeps,H);
244             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
245             VV               = _mm_add_ps(Y,_mm_mul_ps(vfeps,Fp));
246             velec            = _mm_mul_ps(qq00,VV);
247             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
248             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq00,FF),_mm_mul_ps(vftabscale,rinv00)));
249
250             /* CUBIC SPLINE TABLE DISPERSION */
251             vfitab           = _mm_add_epi32(vfitab,ifour);
252             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
253             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
254             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
255             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
256             _MM_TRANSPOSE4_PS(Y,F,G,H);
257             Heps             = _mm_mul_ps(vfeps,H);
258             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
259             VV               = _mm_add_ps(Y,_mm_mul_ps(vfeps,Fp));
260             vvdw6            = _mm_mul_ps(c6_00,VV);
261             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
262             fvdw6            = _mm_mul_ps(c6_00,FF);
263
264             /* CUBIC SPLINE TABLE REPULSION */
265             vfitab           = _mm_add_epi32(vfitab,ifour);
266             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
267             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
268             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
269             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
270             _MM_TRANSPOSE4_PS(Y,F,G,H);
271             Heps             = _mm_mul_ps(vfeps,H);
272             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
273             VV               = _mm_add_ps(Y,_mm_mul_ps(vfeps,Fp));
274             vvdw12           = _mm_mul_ps(c12_00,VV);
275             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
276             fvdw12           = _mm_mul_ps(c12_00,FF);
277             vvdw             = _mm_add_ps(vvdw12,vvdw6);
278             fvdw             = _mm_xor_ps(signbit,_mm_mul_ps(_mm_add_ps(fvdw6,fvdw12),_mm_mul_ps(vftabscale,rinv00)));
279
280             /* Update potential sum for this i atom from the interaction with this j atom. */
281             velecsum         = _mm_add_ps(velecsum,velec);
282             vvdwsum          = _mm_add_ps(vvdwsum,vvdw);
283
284             fscal            = _mm_add_ps(felec,fvdw);
285
286             /* Calculate temporary vectorial force */
287             tx               = _mm_mul_ps(fscal,dx00);
288             ty               = _mm_mul_ps(fscal,dy00);
289             tz               = _mm_mul_ps(fscal,dz00);
290
291             /* Update vectorial force */
292             fix0             = _mm_add_ps(fix0,tx);
293             fiy0             = _mm_add_ps(fiy0,ty);
294             fiz0             = _mm_add_ps(fiz0,tz);
295
296             fjptrA             = f+j_coord_offsetA;
297             fjptrB             = f+j_coord_offsetB;
298             fjptrC             = f+j_coord_offsetC;
299             fjptrD             = f+j_coord_offsetD;
300             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
301             
302             /**************************
303              * CALCULATE INTERACTIONS *
304              **************************/
305
306             r10              = _mm_mul_ps(rsq10,rinv10);
307
308             /* Compute parameters for interactions between i and j atoms */
309             qq10             = _mm_mul_ps(iq1,jq0);
310
311             /* Calculate table index by multiplying r with table scale and truncate to integer */
312             rt               = _mm_mul_ps(r10,vftabscale);
313             vfitab           = _mm_cvttps_epi32(rt);
314             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
315             vfitab           = _mm_slli_epi32(_mm_add_epi32(vfitab,_mm_slli_epi32(vfitab,1)),2);
316
317             /* CUBIC SPLINE TABLE ELECTROSTATICS */
318             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
319             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
320             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
321             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
322             _MM_TRANSPOSE4_PS(Y,F,G,H);
323             Heps             = _mm_mul_ps(vfeps,H);
324             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
325             VV               = _mm_add_ps(Y,_mm_mul_ps(vfeps,Fp));
326             velec            = _mm_mul_ps(qq10,VV);
327             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
328             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq10,FF),_mm_mul_ps(vftabscale,rinv10)));
329
330             /* Update potential sum for this i atom from the interaction with this j atom. */
331             velecsum         = _mm_add_ps(velecsum,velec);
332
333             fscal            = felec;
334
335             /* Calculate temporary vectorial force */
336             tx               = _mm_mul_ps(fscal,dx10);
337             ty               = _mm_mul_ps(fscal,dy10);
338             tz               = _mm_mul_ps(fscal,dz10);
339
340             /* Update vectorial force */
341             fix1             = _mm_add_ps(fix1,tx);
342             fiy1             = _mm_add_ps(fiy1,ty);
343             fiz1             = _mm_add_ps(fiz1,tz);
344
345             fjptrA             = f+j_coord_offsetA;
346             fjptrB             = f+j_coord_offsetB;
347             fjptrC             = f+j_coord_offsetC;
348             fjptrD             = f+j_coord_offsetD;
349             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
350             
351             /**************************
352              * CALCULATE INTERACTIONS *
353              **************************/
354
355             r20              = _mm_mul_ps(rsq20,rinv20);
356
357             /* Compute parameters for interactions between i and j atoms */
358             qq20             = _mm_mul_ps(iq2,jq0);
359
360             /* Calculate table index by multiplying r with table scale and truncate to integer */
361             rt               = _mm_mul_ps(r20,vftabscale);
362             vfitab           = _mm_cvttps_epi32(rt);
363             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
364             vfitab           = _mm_slli_epi32(_mm_add_epi32(vfitab,_mm_slli_epi32(vfitab,1)),2);
365
366             /* CUBIC SPLINE TABLE ELECTROSTATICS */
367             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
368             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
369             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
370             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
371             _MM_TRANSPOSE4_PS(Y,F,G,H);
372             Heps             = _mm_mul_ps(vfeps,H);
373             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
374             VV               = _mm_add_ps(Y,_mm_mul_ps(vfeps,Fp));
375             velec            = _mm_mul_ps(qq20,VV);
376             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
377             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq20,FF),_mm_mul_ps(vftabscale,rinv20)));
378
379             /* Update potential sum for this i atom from the interaction with this j atom. */
380             velecsum         = _mm_add_ps(velecsum,velec);
381
382             fscal            = felec;
383
384             /* Calculate temporary vectorial force */
385             tx               = _mm_mul_ps(fscal,dx20);
386             ty               = _mm_mul_ps(fscal,dy20);
387             tz               = _mm_mul_ps(fscal,dz20);
388
389             /* Update vectorial force */
390             fix2             = _mm_add_ps(fix2,tx);
391             fiy2             = _mm_add_ps(fiy2,ty);
392             fiz2             = _mm_add_ps(fiz2,tz);
393
394             fjptrA             = f+j_coord_offsetA;
395             fjptrB             = f+j_coord_offsetB;
396             fjptrC             = f+j_coord_offsetC;
397             fjptrD             = f+j_coord_offsetD;
398             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
399             
400             /* Inner loop uses 159 flops */
401         }
402
403         if(jidx<j_index_end)
404         {
405
406             /* Get j neighbor index, and coordinate index */
407             jnrlistA         = jjnr[jidx];
408             jnrlistB         = jjnr[jidx+1];
409             jnrlistC         = jjnr[jidx+2];
410             jnrlistD         = jjnr[jidx+3];
411             /* Sign of each element will be negative for non-real atoms.
412              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
413              * so use it as val = _mm_andnot_ps(mask,val) to clear dummy entries.
414              */
415             dummy_mask = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
416             jnrA       = (jnrlistA>=0) ? jnrlistA : 0;
417             jnrB       = (jnrlistB>=0) ? jnrlistB : 0;
418             jnrC       = (jnrlistC>=0) ? jnrlistC : 0;
419             jnrD       = (jnrlistD>=0) ? jnrlistD : 0;
420             j_coord_offsetA  = DIM*jnrA;
421             j_coord_offsetB  = DIM*jnrB;
422             j_coord_offsetC  = DIM*jnrC;
423             j_coord_offsetD  = DIM*jnrD;
424
425             /* load j atom coordinates */
426             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
427                                               x+j_coord_offsetC,x+j_coord_offsetD,
428                                               &jx0,&jy0,&jz0);
429
430             /* Calculate displacement vector */
431             dx00             = _mm_sub_ps(ix0,jx0);
432             dy00             = _mm_sub_ps(iy0,jy0);
433             dz00             = _mm_sub_ps(iz0,jz0);
434             dx10             = _mm_sub_ps(ix1,jx0);
435             dy10             = _mm_sub_ps(iy1,jy0);
436             dz10             = _mm_sub_ps(iz1,jz0);
437             dx20             = _mm_sub_ps(ix2,jx0);
438             dy20             = _mm_sub_ps(iy2,jy0);
439             dz20             = _mm_sub_ps(iz2,jz0);
440
441             /* Calculate squared distance and things based on it */
442             rsq00            = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
443             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
444             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
445
446             rinv00           = gmx_mm_invsqrt_ps(rsq00);
447             rinv10           = gmx_mm_invsqrt_ps(rsq10);
448             rinv20           = gmx_mm_invsqrt_ps(rsq20);
449
450             /* Load parameters for j particles */
451             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
452                                                               charge+jnrC+0,charge+jnrD+0);
453             vdwjidx0A        = 2*vdwtype[jnrA+0];
454             vdwjidx0B        = 2*vdwtype[jnrB+0];
455             vdwjidx0C        = 2*vdwtype[jnrC+0];
456             vdwjidx0D        = 2*vdwtype[jnrD+0];
457
458             /**************************
459              * CALCULATE INTERACTIONS *
460              **************************/
461
462             r00              = _mm_mul_ps(rsq00,rinv00);
463             r00              = _mm_andnot_ps(dummy_mask,r00);
464
465             /* Compute parameters for interactions between i and j atoms */
466             qq00             = _mm_mul_ps(iq0,jq0);
467             gmx_mm_load_4pair_swizzle_ps(vdwparam+vdwioffset0+vdwjidx0A,
468                                          vdwparam+vdwioffset0+vdwjidx0B,
469                                          vdwparam+vdwioffset0+vdwjidx0C,
470                                          vdwparam+vdwioffset0+vdwjidx0D,
471                                          &c6_00,&c12_00);
472
473             /* Calculate table index by multiplying r with table scale and truncate to integer */
474             rt               = _mm_mul_ps(r00,vftabscale);
475             vfitab           = _mm_cvttps_epi32(rt);
476             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
477             vfitab           = _mm_slli_epi32(_mm_add_epi32(vfitab,_mm_slli_epi32(vfitab,1)),2);
478
479             /* CUBIC SPLINE TABLE ELECTROSTATICS */
480             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
481             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
482             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
483             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
484             _MM_TRANSPOSE4_PS(Y,F,G,H);
485             Heps             = _mm_mul_ps(vfeps,H);
486             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
487             VV               = _mm_add_ps(Y,_mm_mul_ps(vfeps,Fp));
488             velec            = _mm_mul_ps(qq00,VV);
489             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
490             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq00,FF),_mm_mul_ps(vftabscale,rinv00)));
491
492             /* CUBIC SPLINE TABLE DISPERSION */
493             vfitab           = _mm_add_epi32(vfitab,ifour);
494             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
495             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
496             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
497             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
498             _MM_TRANSPOSE4_PS(Y,F,G,H);
499             Heps             = _mm_mul_ps(vfeps,H);
500             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
501             VV               = _mm_add_ps(Y,_mm_mul_ps(vfeps,Fp));
502             vvdw6            = _mm_mul_ps(c6_00,VV);
503             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
504             fvdw6            = _mm_mul_ps(c6_00,FF);
505
506             /* CUBIC SPLINE TABLE REPULSION */
507             vfitab           = _mm_add_epi32(vfitab,ifour);
508             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
509             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
510             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
511             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
512             _MM_TRANSPOSE4_PS(Y,F,G,H);
513             Heps             = _mm_mul_ps(vfeps,H);
514             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
515             VV               = _mm_add_ps(Y,_mm_mul_ps(vfeps,Fp));
516             vvdw12           = _mm_mul_ps(c12_00,VV);
517             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
518             fvdw12           = _mm_mul_ps(c12_00,FF);
519             vvdw             = _mm_add_ps(vvdw12,vvdw6);
520             fvdw             = _mm_xor_ps(signbit,_mm_mul_ps(_mm_add_ps(fvdw6,fvdw12),_mm_mul_ps(vftabscale,rinv00)));
521
522             /* Update potential sum for this i atom from the interaction with this j atom. */
523             velec            = _mm_andnot_ps(dummy_mask,velec);
524             velecsum         = _mm_add_ps(velecsum,velec);
525             vvdw             = _mm_andnot_ps(dummy_mask,vvdw);
526             vvdwsum          = _mm_add_ps(vvdwsum,vvdw);
527
528             fscal            = _mm_add_ps(felec,fvdw);
529
530             fscal            = _mm_andnot_ps(dummy_mask,fscal);
531
532             /* Calculate temporary vectorial force */
533             tx               = _mm_mul_ps(fscal,dx00);
534             ty               = _mm_mul_ps(fscal,dy00);
535             tz               = _mm_mul_ps(fscal,dz00);
536
537             /* Update vectorial force */
538             fix0             = _mm_add_ps(fix0,tx);
539             fiy0             = _mm_add_ps(fiy0,ty);
540             fiz0             = _mm_add_ps(fiz0,tz);
541
542             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
543             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
544             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
545             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
546             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
547             
548             /**************************
549              * CALCULATE INTERACTIONS *
550              **************************/
551
552             r10              = _mm_mul_ps(rsq10,rinv10);
553             r10              = _mm_andnot_ps(dummy_mask,r10);
554
555             /* Compute parameters for interactions between i and j atoms */
556             qq10             = _mm_mul_ps(iq1,jq0);
557
558             /* Calculate table index by multiplying r with table scale and truncate to integer */
559             rt               = _mm_mul_ps(r10,vftabscale);
560             vfitab           = _mm_cvttps_epi32(rt);
561             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
562             vfitab           = _mm_slli_epi32(_mm_add_epi32(vfitab,_mm_slli_epi32(vfitab,1)),2);
563
564             /* CUBIC SPLINE TABLE ELECTROSTATICS */
565             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
566             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
567             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
568             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
569             _MM_TRANSPOSE4_PS(Y,F,G,H);
570             Heps             = _mm_mul_ps(vfeps,H);
571             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
572             VV               = _mm_add_ps(Y,_mm_mul_ps(vfeps,Fp));
573             velec            = _mm_mul_ps(qq10,VV);
574             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
575             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq10,FF),_mm_mul_ps(vftabscale,rinv10)));
576
577             /* Update potential sum for this i atom from the interaction with this j atom. */
578             velec            = _mm_andnot_ps(dummy_mask,velec);
579             velecsum         = _mm_add_ps(velecsum,velec);
580
581             fscal            = felec;
582
583             fscal            = _mm_andnot_ps(dummy_mask,fscal);
584
585             /* Calculate temporary vectorial force */
586             tx               = _mm_mul_ps(fscal,dx10);
587             ty               = _mm_mul_ps(fscal,dy10);
588             tz               = _mm_mul_ps(fscal,dz10);
589
590             /* Update vectorial force */
591             fix1             = _mm_add_ps(fix1,tx);
592             fiy1             = _mm_add_ps(fiy1,ty);
593             fiz1             = _mm_add_ps(fiz1,tz);
594
595             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
596             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
597             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
598             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
599             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
600             
601             /**************************
602              * CALCULATE INTERACTIONS *
603              **************************/
604
605             r20              = _mm_mul_ps(rsq20,rinv20);
606             r20              = _mm_andnot_ps(dummy_mask,r20);
607
608             /* Compute parameters for interactions between i and j atoms */
609             qq20             = _mm_mul_ps(iq2,jq0);
610
611             /* Calculate table index by multiplying r with table scale and truncate to integer */
612             rt               = _mm_mul_ps(r20,vftabscale);
613             vfitab           = _mm_cvttps_epi32(rt);
614             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
615             vfitab           = _mm_slli_epi32(_mm_add_epi32(vfitab,_mm_slli_epi32(vfitab,1)),2);
616
617             /* CUBIC SPLINE TABLE ELECTROSTATICS */
618             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
619             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
620             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
621             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
622             _MM_TRANSPOSE4_PS(Y,F,G,H);
623             Heps             = _mm_mul_ps(vfeps,H);
624             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
625             VV               = _mm_add_ps(Y,_mm_mul_ps(vfeps,Fp));
626             velec            = _mm_mul_ps(qq20,VV);
627             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
628             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq20,FF),_mm_mul_ps(vftabscale,rinv20)));
629
630             /* Update potential sum for this i atom from the interaction with this j atom. */
631             velec            = _mm_andnot_ps(dummy_mask,velec);
632             velecsum         = _mm_add_ps(velecsum,velec);
633
634             fscal            = felec;
635
636             fscal            = _mm_andnot_ps(dummy_mask,fscal);
637
638             /* Calculate temporary vectorial force */
639             tx               = _mm_mul_ps(fscal,dx20);
640             ty               = _mm_mul_ps(fscal,dy20);
641             tz               = _mm_mul_ps(fscal,dz20);
642
643             /* Update vectorial force */
644             fix2             = _mm_add_ps(fix2,tx);
645             fiy2             = _mm_add_ps(fiy2,ty);
646             fiz2             = _mm_add_ps(fiz2,tz);
647
648             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
649             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
650             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
651             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
652             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
653             
654             /* Inner loop uses 162 flops */
655         }
656
657         /* End of innermost loop */
658
659         gmx_mm_update_iforce_3atom_swizzle_ps(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,
660                                               f+i_coord_offset,fshift+i_shift_offset);
661
662         ggid                        = gid[iidx];
663         /* Update potential energies */
664         gmx_mm_update_1pot_ps(velecsum,kernel_data->energygrp_elec+ggid);
665         gmx_mm_update_1pot_ps(vvdwsum,kernel_data->energygrp_vdw+ggid);
666
667         /* Increment number of inner iterations */
668         inneriter                  += j_index_end - j_index_start;
669
670         /* Outer loop uses 20 flops */
671     }
672
673     /* Increment number of outer iterations */
674     outeriter        += nri;
675
676     /* Update outer/inner flops */
677
678     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W3_VF,outeriter*20 + inneriter*162);
679 }
680 /*
681  * Gromacs nonbonded kernel:   nb_kernel_ElecCSTab_VdwCSTab_GeomW3P1_F_sse2_single
682  * Electrostatics interaction: CubicSplineTable
683  * VdW interaction:            CubicSplineTable
684  * Geometry:                   Water3-Particle
685  * Calculate force/pot:        Force
686  */
687 void
688 nb_kernel_ElecCSTab_VdwCSTab_GeomW3P1_F_sse2_single
689                     (t_nblist * gmx_restrict                nlist,
690                      rvec * gmx_restrict                    xx,
691                      rvec * gmx_restrict                    ff,
692                      t_forcerec * gmx_restrict              fr,
693                      t_mdatoms * gmx_restrict               mdatoms,
694                      nb_kernel_data_t * gmx_restrict        kernel_data,
695                      t_nrnb * gmx_restrict                  nrnb)
696 {
697     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or 
698      * just 0 for non-waters.
699      * Suffixes A,B,C,D refer to j loop unrolling done with SSE, e.g. for the four different
700      * jnr indices corresponding to data put in the four positions in the SIMD register.
701      */
702     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
703     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
704     int              jnrA,jnrB,jnrC,jnrD;
705     int              jnrlistA,jnrlistB,jnrlistC,jnrlistD;
706     int              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
707     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
708     real             rcutoff_scalar;
709     real             *shiftvec,*fshift,*x,*f;
710     real             *fjptrA,*fjptrB,*fjptrC,*fjptrD;
711     real             scratch[4*DIM];
712     __m128           tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
713     int              vdwioffset0;
714     __m128           ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
715     int              vdwioffset1;
716     __m128           ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
717     int              vdwioffset2;
718     __m128           ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
719     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
720     __m128           jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
721     __m128           dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
722     __m128           dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
723     __m128           dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
724     __m128           velec,felec,velecsum,facel,crf,krf,krf2;
725     real             *charge;
726     int              nvdwtype;
727     __m128           rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
728     int              *vdwtype;
729     real             *vdwparam;
730     __m128           one_sixth   = _mm_set1_ps(1.0/6.0);
731     __m128           one_twelfth = _mm_set1_ps(1.0/12.0);
732     __m128i          vfitab;
733     __m128i          ifour       = _mm_set1_epi32(4);
734     __m128           rt,vfeps,vftabscale,Y,F,G,H,Heps,Fp,VV,FF;
735     real             *vftab;
736     __m128           dummy_mask,cutoff_mask;
737     __m128           signbit = _mm_castsi128_ps( _mm_set1_epi32(0x80000000) );
738     __m128           one     = _mm_set1_ps(1.0);
739     __m128           two     = _mm_set1_ps(2.0);
740     x                = xx[0];
741     f                = ff[0];
742
743     nri              = nlist->nri;
744     iinr             = nlist->iinr;
745     jindex           = nlist->jindex;
746     jjnr             = nlist->jjnr;
747     shiftidx         = nlist->shift;
748     gid              = nlist->gid;
749     shiftvec         = fr->shift_vec[0];
750     fshift           = fr->fshift[0];
751     facel            = _mm_set1_ps(fr->epsfac);
752     charge           = mdatoms->chargeA;
753     nvdwtype         = fr->ntype;
754     vdwparam         = fr->nbfp;
755     vdwtype          = mdatoms->typeA;
756
757     vftab            = kernel_data->table_elec_vdw->data;
758     vftabscale       = _mm_set1_ps(kernel_data->table_elec_vdw->scale);
759
760     /* Setup water-specific parameters */
761     inr              = nlist->iinr[0];
762     iq0              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+0]));
763     iq1              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+1]));
764     iq2              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+2]));
765     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
766
767     /* Avoid stupid compiler warnings */
768     jnrA = jnrB = jnrC = jnrD = 0;
769     j_coord_offsetA = 0;
770     j_coord_offsetB = 0;
771     j_coord_offsetC = 0;
772     j_coord_offsetD = 0;
773
774     outeriter        = 0;
775     inneriter        = 0;
776
777     for(iidx=0;iidx<4*DIM;iidx++)
778     {
779         scratch[iidx] = 0.0;
780     }  
781
782     /* Start outer loop over neighborlists */
783     for(iidx=0; iidx<nri; iidx++)
784     {
785         /* Load shift vector for this list */
786         i_shift_offset   = DIM*shiftidx[iidx];
787
788         /* Load limits for loop over neighbors */
789         j_index_start    = jindex[iidx];
790         j_index_end      = jindex[iidx+1];
791
792         /* Get outer coordinate index */
793         inr              = iinr[iidx];
794         i_coord_offset   = DIM*inr;
795
796         /* Load i particle coords and add shift vector */
797         gmx_mm_load_shift_and_3rvec_broadcast_ps(shiftvec+i_shift_offset,x+i_coord_offset,
798                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2);
799         
800         fix0             = _mm_setzero_ps();
801         fiy0             = _mm_setzero_ps();
802         fiz0             = _mm_setzero_ps();
803         fix1             = _mm_setzero_ps();
804         fiy1             = _mm_setzero_ps();
805         fiz1             = _mm_setzero_ps();
806         fix2             = _mm_setzero_ps();
807         fiy2             = _mm_setzero_ps();
808         fiz2             = _mm_setzero_ps();
809
810         /* Start inner kernel loop */
811         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
812         {
813
814             /* Get j neighbor index, and coordinate index */
815             jnrA             = jjnr[jidx];
816             jnrB             = jjnr[jidx+1];
817             jnrC             = jjnr[jidx+2];
818             jnrD             = jjnr[jidx+3];
819             j_coord_offsetA  = DIM*jnrA;
820             j_coord_offsetB  = DIM*jnrB;
821             j_coord_offsetC  = DIM*jnrC;
822             j_coord_offsetD  = DIM*jnrD;
823
824             /* load j atom coordinates */
825             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
826                                               x+j_coord_offsetC,x+j_coord_offsetD,
827                                               &jx0,&jy0,&jz0);
828
829             /* Calculate displacement vector */
830             dx00             = _mm_sub_ps(ix0,jx0);
831             dy00             = _mm_sub_ps(iy0,jy0);
832             dz00             = _mm_sub_ps(iz0,jz0);
833             dx10             = _mm_sub_ps(ix1,jx0);
834             dy10             = _mm_sub_ps(iy1,jy0);
835             dz10             = _mm_sub_ps(iz1,jz0);
836             dx20             = _mm_sub_ps(ix2,jx0);
837             dy20             = _mm_sub_ps(iy2,jy0);
838             dz20             = _mm_sub_ps(iz2,jz0);
839
840             /* Calculate squared distance and things based on it */
841             rsq00            = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
842             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
843             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
844
845             rinv00           = gmx_mm_invsqrt_ps(rsq00);
846             rinv10           = gmx_mm_invsqrt_ps(rsq10);
847             rinv20           = gmx_mm_invsqrt_ps(rsq20);
848
849             /* Load parameters for j particles */
850             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
851                                                               charge+jnrC+0,charge+jnrD+0);
852             vdwjidx0A        = 2*vdwtype[jnrA+0];
853             vdwjidx0B        = 2*vdwtype[jnrB+0];
854             vdwjidx0C        = 2*vdwtype[jnrC+0];
855             vdwjidx0D        = 2*vdwtype[jnrD+0];
856
857             /**************************
858              * CALCULATE INTERACTIONS *
859              **************************/
860
861             r00              = _mm_mul_ps(rsq00,rinv00);
862
863             /* Compute parameters for interactions between i and j atoms */
864             qq00             = _mm_mul_ps(iq0,jq0);
865             gmx_mm_load_4pair_swizzle_ps(vdwparam+vdwioffset0+vdwjidx0A,
866                                          vdwparam+vdwioffset0+vdwjidx0B,
867                                          vdwparam+vdwioffset0+vdwjidx0C,
868                                          vdwparam+vdwioffset0+vdwjidx0D,
869                                          &c6_00,&c12_00);
870
871             /* Calculate table index by multiplying r with table scale and truncate to integer */
872             rt               = _mm_mul_ps(r00,vftabscale);
873             vfitab           = _mm_cvttps_epi32(rt);
874             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
875             vfitab           = _mm_slli_epi32(_mm_add_epi32(vfitab,_mm_slli_epi32(vfitab,1)),2);
876
877             /* CUBIC SPLINE TABLE ELECTROSTATICS */
878             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
879             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
880             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
881             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
882             _MM_TRANSPOSE4_PS(Y,F,G,H);
883             Heps             = _mm_mul_ps(vfeps,H);
884             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
885             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
886             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq00,FF),_mm_mul_ps(vftabscale,rinv00)));
887
888             /* CUBIC SPLINE TABLE DISPERSION */
889             vfitab           = _mm_add_epi32(vfitab,ifour);
890             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
891             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
892             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
893             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
894             _MM_TRANSPOSE4_PS(Y,F,G,H);
895             Heps             = _mm_mul_ps(vfeps,H);
896             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
897             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
898             fvdw6            = _mm_mul_ps(c6_00,FF);
899
900             /* CUBIC SPLINE TABLE REPULSION */
901             vfitab           = _mm_add_epi32(vfitab,ifour);
902             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
903             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
904             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
905             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
906             _MM_TRANSPOSE4_PS(Y,F,G,H);
907             Heps             = _mm_mul_ps(vfeps,H);
908             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
909             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
910             fvdw12           = _mm_mul_ps(c12_00,FF);
911             fvdw             = _mm_xor_ps(signbit,_mm_mul_ps(_mm_add_ps(fvdw6,fvdw12),_mm_mul_ps(vftabscale,rinv00)));
912
913             fscal            = _mm_add_ps(felec,fvdw);
914
915             /* Calculate temporary vectorial force */
916             tx               = _mm_mul_ps(fscal,dx00);
917             ty               = _mm_mul_ps(fscal,dy00);
918             tz               = _mm_mul_ps(fscal,dz00);
919
920             /* Update vectorial force */
921             fix0             = _mm_add_ps(fix0,tx);
922             fiy0             = _mm_add_ps(fiy0,ty);
923             fiz0             = _mm_add_ps(fiz0,tz);
924
925             fjptrA             = f+j_coord_offsetA;
926             fjptrB             = f+j_coord_offsetB;
927             fjptrC             = f+j_coord_offsetC;
928             fjptrD             = f+j_coord_offsetD;
929             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
930             
931             /**************************
932              * CALCULATE INTERACTIONS *
933              **************************/
934
935             r10              = _mm_mul_ps(rsq10,rinv10);
936
937             /* Compute parameters for interactions between i and j atoms */
938             qq10             = _mm_mul_ps(iq1,jq0);
939
940             /* Calculate table index by multiplying r with table scale and truncate to integer */
941             rt               = _mm_mul_ps(r10,vftabscale);
942             vfitab           = _mm_cvttps_epi32(rt);
943             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
944             vfitab           = _mm_slli_epi32(_mm_add_epi32(vfitab,_mm_slli_epi32(vfitab,1)),2);
945
946             /* CUBIC SPLINE TABLE ELECTROSTATICS */
947             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
948             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
949             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
950             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
951             _MM_TRANSPOSE4_PS(Y,F,G,H);
952             Heps             = _mm_mul_ps(vfeps,H);
953             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
954             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
955             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq10,FF),_mm_mul_ps(vftabscale,rinv10)));
956
957             fscal            = felec;
958
959             /* Calculate temporary vectorial force */
960             tx               = _mm_mul_ps(fscal,dx10);
961             ty               = _mm_mul_ps(fscal,dy10);
962             tz               = _mm_mul_ps(fscal,dz10);
963
964             /* Update vectorial force */
965             fix1             = _mm_add_ps(fix1,tx);
966             fiy1             = _mm_add_ps(fiy1,ty);
967             fiz1             = _mm_add_ps(fiz1,tz);
968
969             fjptrA             = f+j_coord_offsetA;
970             fjptrB             = f+j_coord_offsetB;
971             fjptrC             = f+j_coord_offsetC;
972             fjptrD             = f+j_coord_offsetD;
973             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
974             
975             /**************************
976              * CALCULATE INTERACTIONS *
977              **************************/
978
979             r20              = _mm_mul_ps(rsq20,rinv20);
980
981             /* Compute parameters for interactions between i and j atoms */
982             qq20             = _mm_mul_ps(iq2,jq0);
983
984             /* Calculate table index by multiplying r with table scale and truncate to integer */
985             rt               = _mm_mul_ps(r20,vftabscale);
986             vfitab           = _mm_cvttps_epi32(rt);
987             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
988             vfitab           = _mm_slli_epi32(_mm_add_epi32(vfitab,_mm_slli_epi32(vfitab,1)),2);
989
990             /* CUBIC SPLINE TABLE ELECTROSTATICS */
991             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
992             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
993             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
994             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
995             _MM_TRANSPOSE4_PS(Y,F,G,H);
996             Heps             = _mm_mul_ps(vfeps,H);
997             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
998             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
999             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq20,FF),_mm_mul_ps(vftabscale,rinv20)));
1000
1001             fscal            = felec;
1002
1003             /* Calculate temporary vectorial force */
1004             tx               = _mm_mul_ps(fscal,dx20);
1005             ty               = _mm_mul_ps(fscal,dy20);
1006             tz               = _mm_mul_ps(fscal,dz20);
1007
1008             /* Update vectorial force */
1009             fix2             = _mm_add_ps(fix2,tx);
1010             fiy2             = _mm_add_ps(fiy2,ty);
1011             fiz2             = _mm_add_ps(fiz2,tz);
1012
1013             fjptrA             = f+j_coord_offsetA;
1014             fjptrB             = f+j_coord_offsetB;
1015             fjptrC             = f+j_coord_offsetC;
1016             fjptrD             = f+j_coord_offsetD;
1017             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
1018             
1019             /* Inner loop uses 139 flops */
1020         }
1021
1022         if(jidx<j_index_end)
1023         {
1024
1025             /* Get j neighbor index, and coordinate index */
1026             jnrlistA         = jjnr[jidx];
1027             jnrlistB         = jjnr[jidx+1];
1028             jnrlistC         = jjnr[jidx+2];
1029             jnrlistD         = jjnr[jidx+3];
1030             /* Sign of each element will be negative for non-real atoms.
1031              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
1032              * so use it as val = _mm_andnot_ps(mask,val) to clear dummy entries.
1033              */
1034             dummy_mask = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
1035             jnrA       = (jnrlistA>=0) ? jnrlistA : 0;
1036             jnrB       = (jnrlistB>=0) ? jnrlistB : 0;
1037             jnrC       = (jnrlistC>=0) ? jnrlistC : 0;
1038             jnrD       = (jnrlistD>=0) ? jnrlistD : 0;
1039             j_coord_offsetA  = DIM*jnrA;
1040             j_coord_offsetB  = DIM*jnrB;
1041             j_coord_offsetC  = DIM*jnrC;
1042             j_coord_offsetD  = DIM*jnrD;
1043
1044             /* load j atom coordinates */
1045             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
1046                                               x+j_coord_offsetC,x+j_coord_offsetD,
1047                                               &jx0,&jy0,&jz0);
1048
1049             /* Calculate displacement vector */
1050             dx00             = _mm_sub_ps(ix0,jx0);
1051             dy00             = _mm_sub_ps(iy0,jy0);
1052             dz00             = _mm_sub_ps(iz0,jz0);
1053             dx10             = _mm_sub_ps(ix1,jx0);
1054             dy10             = _mm_sub_ps(iy1,jy0);
1055             dz10             = _mm_sub_ps(iz1,jz0);
1056             dx20             = _mm_sub_ps(ix2,jx0);
1057             dy20             = _mm_sub_ps(iy2,jy0);
1058             dz20             = _mm_sub_ps(iz2,jz0);
1059
1060             /* Calculate squared distance and things based on it */
1061             rsq00            = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
1062             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
1063             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
1064
1065             rinv00           = gmx_mm_invsqrt_ps(rsq00);
1066             rinv10           = gmx_mm_invsqrt_ps(rsq10);
1067             rinv20           = gmx_mm_invsqrt_ps(rsq20);
1068
1069             /* Load parameters for j particles */
1070             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
1071                                                               charge+jnrC+0,charge+jnrD+0);
1072             vdwjidx0A        = 2*vdwtype[jnrA+0];
1073             vdwjidx0B        = 2*vdwtype[jnrB+0];
1074             vdwjidx0C        = 2*vdwtype[jnrC+0];
1075             vdwjidx0D        = 2*vdwtype[jnrD+0];
1076
1077             /**************************
1078              * CALCULATE INTERACTIONS *
1079              **************************/
1080
1081             r00              = _mm_mul_ps(rsq00,rinv00);
1082             r00              = _mm_andnot_ps(dummy_mask,r00);
1083
1084             /* Compute parameters for interactions between i and j atoms */
1085             qq00             = _mm_mul_ps(iq0,jq0);
1086             gmx_mm_load_4pair_swizzle_ps(vdwparam+vdwioffset0+vdwjidx0A,
1087                                          vdwparam+vdwioffset0+vdwjidx0B,
1088                                          vdwparam+vdwioffset0+vdwjidx0C,
1089                                          vdwparam+vdwioffset0+vdwjidx0D,
1090                                          &c6_00,&c12_00);
1091
1092             /* Calculate table index by multiplying r with table scale and truncate to integer */
1093             rt               = _mm_mul_ps(r00,vftabscale);
1094             vfitab           = _mm_cvttps_epi32(rt);
1095             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
1096             vfitab           = _mm_slli_epi32(_mm_add_epi32(vfitab,_mm_slli_epi32(vfitab,1)),2);
1097
1098             /* CUBIC SPLINE TABLE ELECTROSTATICS */
1099             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
1100             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
1101             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
1102             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
1103             _MM_TRANSPOSE4_PS(Y,F,G,H);
1104             Heps             = _mm_mul_ps(vfeps,H);
1105             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
1106             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
1107             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq00,FF),_mm_mul_ps(vftabscale,rinv00)));
1108
1109             /* CUBIC SPLINE TABLE DISPERSION */
1110             vfitab           = _mm_add_epi32(vfitab,ifour);
1111             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
1112             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
1113             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
1114             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
1115             _MM_TRANSPOSE4_PS(Y,F,G,H);
1116             Heps             = _mm_mul_ps(vfeps,H);
1117             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
1118             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
1119             fvdw6            = _mm_mul_ps(c6_00,FF);
1120
1121             /* CUBIC SPLINE TABLE REPULSION */
1122             vfitab           = _mm_add_epi32(vfitab,ifour);
1123             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
1124             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
1125             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
1126             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
1127             _MM_TRANSPOSE4_PS(Y,F,G,H);
1128             Heps             = _mm_mul_ps(vfeps,H);
1129             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
1130             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
1131             fvdw12           = _mm_mul_ps(c12_00,FF);
1132             fvdw             = _mm_xor_ps(signbit,_mm_mul_ps(_mm_add_ps(fvdw6,fvdw12),_mm_mul_ps(vftabscale,rinv00)));
1133
1134             fscal            = _mm_add_ps(felec,fvdw);
1135
1136             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1137
1138             /* Calculate temporary vectorial force */
1139             tx               = _mm_mul_ps(fscal,dx00);
1140             ty               = _mm_mul_ps(fscal,dy00);
1141             tz               = _mm_mul_ps(fscal,dz00);
1142
1143             /* Update vectorial force */
1144             fix0             = _mm_add_ps(fix0,tx);
1145             fiy0             = _mm_add_ps(fiy0,ty);
1146             fiz0             = _mm_add_ps(fiz0,tz);
1147
1148             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
1149             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
1150             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
1151             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
1152             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
1153             
1154             /**************************
1155              * CALCULATE INTERACTIONS *
1156              **************************/
1157
1158             r10              = _mm_mul_ps(rsq10,rinv10);
1159             r10              = _mm_andnot_ps(dummy_mask,r10);
1160
1161             /* Compute parameters for interactions between i and j atoms */
1162             qq10             = _mm_mul_ps(iq1,jq0);
1163
1164             /* Calculate table index by multiplying r with table scale and truncate to integer */
1165             rt               = _mm_mul_ps(r10,vftabscale);
1166             vfitab           = _mm_cvttps_epi32(rt);
1167             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
1168             vfitab           = _mm_slli_epi32(_mm_add_epi32(vfitab,_mm_slli_epi32(vfitab,1)),2);
1169
1170             /* CUBIC SPLINE TABLE ELECTROSTATICS */
1171             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
1172             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
1173             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
1174             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
1175             _MM_TRANSPOSE4_PS(Y,F,G,H);
1176             Heps             = _mm_mul_ps(vfeps,H);
1177             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
1178             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
1179             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq10,FF),_mm_mul_ps(vftabscale,rinv10)));
1180
1181             fscal            = felec;
1182
1183             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1184
1185             /* Calculate temporary vectorial force */
1186             tx               = _mm_mul_ps(fscal,dx10);
1187             ty               = _mm_mul_ps(fscal,dy10);
1188             tz               = _mm_mul_ps(fscal,dz10);
1189
1190             /* Update vectorial force */
1191             fix1             = _mm_add_ps(fix1,tx);
1192             fiy1             = _mm_add_ps(fiy1,ty);
1193             fiz1             = _mm_add_ps(fiz1,tz);
1194
1195             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
1196             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
1197             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
1198             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
1199             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
1200             
1201             /**************************
1202              * CALCULATE INTERACTIONS *
1203              **************************/
1204
1205             r20              = _mm_mul_ps(rsq20,rinv20);
1206             r20              = _mm_andnot_ps(dummy_mask,r20);
1207
1208             /* Compute parameters for interactions between i and j atoms */
1209             qq20             = _mm_mul_ps(iq2,jq0);
1210
1211             /* Calculate table index by multiplying r with table scale and truncate to integer */
1212             rt               = _mm_mul_ps(r20,vftabscale);
1213             vfitab           = _mm_cvttps_epi32(rt);
1214             vfeps            = _mm_sub_ps(rt,_mm_cvtepi32_ps(vfitab));
1215             vfitab           = _mm_slli_epi32(_mm_add_epi32(vfitab,_mm_slli_epi32(vfitab,1)),2);
1216
1217             /* CUBIC SPLINE TABLE ELECTROSTATICS */
1218             Y                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,0) );
1219             F                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,1) );
1220             G                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,2) );
1221             H                = _mm_load_ps( vftab + gmx_mm_extract_epi32(vfitab,3) );
1222             _MM_TRANSPOSE4_PS(Y,F,G,H);
1223             Heps             = _mm_mul_ps(vfeps,H);
1224             Fp               = _mm_add_ps(F,_mm_mul_ps(vfeps,_mm_add_ps(G,Heps)));
1225             FF               = _mm_add_ps(Fp,_mm_mul_ps(vfeps,_mm_add_ps(G,_mm_add_ps(Heps,Heps))));
1226             felec            = _mm_xor_ps(signbit,_mm_mul_ps(_mm_mul_ps(qq20,FF),_mm_mul_ps(vftabscale,rinv20)));
1227
1228             fscal            = felec;
1229
1230             fscal            = _mm_andnot_ps(dummy_mask,fscal);
1231
1232             /* Calculate temporary vectorial force */
1233             tx               = _mm_mul_ps(fscal,dx20);
1234             ty               = _mm_mul_ps(fscal,dy20);
1235             tz               = _mm_mul_ps(fscal,dz20);
1236
1237             /* Update vectorial force */
1238             fix2             = _mm_add_ps(fix2,tx);
1239             fiy2             = _mm_add_ps(fiy2,ty);
1240             fiz2             = _mm_add_ps(fiz2,tz);
1241
1242             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
1243             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
1244             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
1245             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
1246             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
1247             
1248             /* Inner loop uses 142 flops */
1249         }
1250
1251         /* End of innermost loop */
1252
1253         gmx_mm_update_iforce_3atom_swizzle_ps(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,
1254                                               f+i_coord_offset,fshift+i_shift_offset);
1255
1256         /* Increment number of inner iterations */
1257         inneriter                  += j_index_end - j_index_start;
1258
1259         /* Outer loop uses 18 flops */
1260     }
1261
1262     /* Increment number of outer iterations */
1263     outeriter        += nri;
1264
1265     /* Update outer/inner flops */
1266
1267     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W3_F,outeriter*18 + inneriter*142);
1268 }