Merge release-5-0 into master
[alexxy/gromacs.git] / src / gromacs / gmxlib / nonbonded / nb_kernel_avx_128_fma_double / nb_kernel_ElecRF_VdwCSTab_GeomW3P1_avx_128_fma_double.c
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6  * and including many others, as listed in the AUTHORS file in the
7  * top-level source directory and at http://www.gromacs.org.
8  *
9  * GROMACS is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public License
11  * as published by the Free Software Foundation; either version 2.1
12  * of the License, or (at your option) any later version.
13  *
14  * GROMACS is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17  * Lesser General Public License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with GROMACS; if not, see
21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
23  *
24  * If you want to redistribute modifications to GROMACS, please
25  * consider that scientific software is very special. Version
26  * control is crucial - bugs must be traceable. We will be happy to
27  * consider code for inclusion in the official distribution, but
28  * derived work must not be called official GROMACS. Details are found
29  * in the README & COPYING files - if they are missing, get the
30  * official version at http://www.gromacs.org.
31  *
32  * To help us fund GROMACS development, we humbly ask that you cite
33  * the research papers on the package. Check out http://www.gromacs.org.
34  */
35 /*
36  * Note: this file was generated by the GROMACS avx_128_fma_double kernel generator.
37  */
38 #include "gmxpre.h"
39
40 #include "config.h"
41
42 #include <math.h>
43
44 #include "../nb_kernel.h"
45 #include "gromacs/legacyheaders/types/simple.h"
46 #include "gromacs/math/vec.h"
47 #include "gromacs/legacyheaders/nrnb.h"
48
49 #include "gromacs/simd/math_x86_avx_128_fma_double.h"
50 #include "kernelutil_x86_avx_128_fma_double.h"
51
52 /*
53  * Gromacs nonbonded kernel:   nb_kernel_ElecRF_VdwCSTab_GeomW3P1_VF_avx_128_fma_double
54  * Electrostatics interaction: ReactionField
55  * VdW interaction:            CubicSplineTable
56  * Geometry:                   Water3-Particle
57  * Calculate force/pot:        PotentialAndForce
58  */
59 void
60 nb_kernel_ElecRF_VdwCSTab_GeomW3P1_VF_avx_128_fma_double
61                     (t_nblist                    * gmx_restrict       nlist,
62                      rvec                        * gmx_restrict          xx,
63                      rvec                        * gmx_restrict          ff,
64                      t_forcerec                  * gmx_restrict          fr,
65                      t_mdatoms                   * gmx_restrict     mdatoms,
66                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
67                      t_nrnb                      * gmx_restrict        nrnb)
68 {
69     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
70      * just 0 for non-waters.
71      * Suffixes A,B refer to j loop unrolling done with SSE double precision, e.g. for the two different
72      * jnr indices corresponding to data put in the four positions in the SIMD register.
73      */
74     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
75     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
76     int              jnrA,jnrB;
77     int              j_coord_offsetA,j_coord_offsetB;
78     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
79     real             rcutoff_scalar;
80     real             *shiftvec,*fshift,*x,*f;
81     __m128d          tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
82     int              vdwioffset0;
83     __m128d          ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
84     int              vdwioffset1;
85     __m128d          ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
86     int              vdwioffset2;
87     __m128d          ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
88     int              vdwjidx0A,vdwjidx0B;
89     __m128d          jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
90     __m128d          dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
91     __m128d          dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
92     __m128d          dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
93     __m128d          velec,felec,velecsum,facel,crf,krf,krf2;
94     real             *charge;
95     int              nvdwtype;
96     __m128d          rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
97     int              *vdwtype;
98     real             *vdwparam;
99     __m128d          one_sixth   = _mm_set1_pd(1.0/6.0);
100     __m128d          one_twelfth = _mm_set1_pd(1.0/12.0);
101     __m128i          vfitab;
102     __m128i          ifour       = _mm_set1_epi32(4);
103     __m128d          rt,vfeps,vftabscale,Y,F,G,H,Heps,Fp,VV,FF,twovfeps;
104     real             *vftab;
105     __m128d          dummy_mask,cutoff_mask;
106     __m128d          signbit   = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
107     __m128d          one     = _mm_set1_pd(1.0);
108     __m128d          two     = _mm_set1_pd(2.0);
109     x                = xx[0];
110     f                = ff[0];
111
112     nri              = nlist->nri;
113     iinr             = nlist->iinr;
114     jindex           = nlist->jindex;
115     jjnr             = nlist->jjnr;
116     shiftidx         = nlist->shift;
117     gid              = nlist->gid;
118     shiftvec         = fr->shift_vec[0];
119     fshift           = fr->fshift[0];
120     facel            = _mm_set1_pd(fr->epsfac);
121     charge           = mdatoms->chargeA;
122     krf              = _mm_set1_pd(fr->ic->k_rf);
123     krf2             = _mm_set1_pd(fr->ic->k_rf*2.0);
124     crf              = _mm_set1_pd(fr->ic->c_rf);
125     nvdwtype         = fr->ntype;
126     vdwparam         = fr->nbfp;
127     vdwtype          = mdatoms->typeA;
128
129     vftab            = kernel_data->table_vdw->data;
130     vftabscale       = _mm_set1_pd(kernel_data->table_vdw->scale);
131
132     /* Setup water-specific parameters */
133     inr              = nlist->iinr[0];
134     iq0              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+0]));
135     iq1              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+1]));
136     iq2              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+2]));
137     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
138
139     /* Avoid stupid compiler warnings */
140     jnrA = jnrB = 0;
141     j_coord_offsetA = 0;
142     j_coord_offsetB = 0;
143
144     outeriter        = 0;
145     inneriter        = 0;
146
147     /* Start outer loop over neighborlists */
148     for(iidx=0; iidx<nri; iidx++)
149     {
150         /* Load shift vector for this list */
151         i_shift_offset   = DIM*shiftidx[iidx];
152
153         /* Load limits for loop over neighbors */
154         j_index_start    = jindex[iidx];
155         j_index_end      = jindex[iidx+1];
156
157         /* Get outer coordinate index */
158         inr              = iinr[iidx];
159         i_coord_offset   = DIM*inr;
160
161         /* Load i particle coords and add shift vector */
162         gmx_mm_load_shift_and_3rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,
163                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2);
164
165         fix0             = _mm_setzero_pd();
166         fiy0             = _mm_setzero_pd();
167         fiz0             = _mm_setzero_pd();
168         fix1             = _mm_setzero_pd();
169         fiy1             = _mm_setzero_pd();
170         fiz1             = _mm_setzero_pd();
171         fix2             = _mm_setzero_pd();
172         fiy2             = _mm_setzero_pd();
173         fiz2             = _mm_setzero_pd();
174
175         /* Reset potential sums */
176         velecsum         = _mm_setzero_pd();
177         vvdwsum          = _mm_setzero_pd();
178
179         /* Start inner kernel loop */
180         for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
181         {
182
183             /* Get j neighbor index, and coordinate index */
184             jnrA             = jjnr[jidx];
185             jnrB             = jjnr[jidx+1];
186             j_coord_offsetA  = DIM*jnrA;
187             j_coord_offsetB  = DIM*jnrB;
188
189             /* load j atom coordinates */
190             gmx_mm_load_1rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
191                                               &jx0,&jy0,&jz0);
192
193             /* Calculate displacement vector */
194             dx00             = _mm_sub_pd(ix0,jx0);
195             dy00             = _mm_sub_pd(iy0,jy0);
196             dz00             = _mm_sub_pd(iz0,jz0);
197             dx10             = _mm_sub_pd(ix1,jx0);
198             dy10             = _mm_sub_pd(iy1,jy0);
199             dz10             = _mm_sub_pd(iz1,jz0);
200             dx20             = _mm_sub_pd(ix2,jx0);
201             dy20             = _mm_sub_pd(iy2,jy0);
202             dz20             = _mm_sub_pd(iz2,jz0);
203
204             /* Calculate squared distance and things based on it */
205             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
206             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
207             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
208
209             rinv00           = gmx_mm_invsqrt_pd(rsq00);
210             rinv10           = gmx_mm_invsqrt_pd(rsq10);
211             rinv20           = gmx_mm_invsqrt_pd(rsq20);
212
213             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
214             rinvsq10         = _mm_mul_pd(rinv10,rinv10);
215             rinvsq20         = _mm_mul_pd(rinv20,rinv20);
216
217             /* Load parameters for j particles */
218             jq0              = gmx_mm_load_2real_swizzle_pd(charge+jnrA+0,charge+jnrB+0);
219             vdwjidx0A        = 2*vdwtype[jnrA+0];
220             vdwjidx0B        = 2*vdwtype[jnrB+0];
221
222             fjx0             = _mm_setzero_pd();
223             fjy0             = _mm_setzero_pd();
224             fjz0             = _mm_setzero_pd();
225
226             /**************************
227              * CALCULATE INTERACTIONS *
228              **************************/
229
230             r00              = _mm_mul_pd(rsq00,rinv00);
231
232             /* Compute parameters for interactions between i and j atoms */
233             qq00             = _mm_mul_pd(iq0,jq0);
234             gmx_mm_load_2pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,
235                                          vdwparam+vdwioffset0+vdwjidx0B,&c6_00,&c12_00);
236
237             /* Calculate table index by multiplying r with table scale and truncate to integer */
238             rt               = _mm_mul_pd(r00,vftabscale);
239             vfitab           = _mm_cvttpd_epi32(rt);
240 #ifdef __XOP__
241             vfeps            = _mm_frcz_pd(rt);
242 #else
243             vfeps            = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
244 #endif
245             twovfeps         = _mm_add_pd(vfeps,vfeps);
246             vfitab           = _mm_slli_epi32(vfitab,3);
247
248             /* REACTION-FIELD ELECTROSTATICS */
249             velec            = _mm_mul_pd(qq00,_mm_sub_pd(_mm_macc_pd(krf,rsq00,rinv00),crf));
250             felec            = _mm_mul_pd(qq00,_mm_msub_pd(rinv00,rinvsq00,krf2));
251
252             /* CUBIC SPLINE TABLE DISPERSION */
253             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
254             F                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) );
255             GMX_MM_TRANSPOSE2_PD(Y,F);
256             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
257             H                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) +2);
258             GMX_MM_TRANSPOSE2_PD(G,H);
259             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(H,vfeps,G),F);
260             VV               = _mm_macc_pd(vfeps,Fp,Y);
261             vvdw6            = _mm_mul_pd(c6_00,VV);
262             FF               = _mm_macc_pd(vfeps,_mm_macc_pd(twovfeps,H,G),Fp);
263             fvdw6            = _mm_mul_pd(c6_00,FF);
264
265             /* CUBIC SPLINE TABLE REPULSION */
266             vfitab           = _mm_add_epi32(vfitab,ifour);
267             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
268             F                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) );
269             GMX_MM_TRANSPOSE2_PD(Y,F);
270             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
271             H                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) +2);
272             GMX_MM_TRANSPOSE2_PD(G,H);
273             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(H,vfeps,G),F);
274             VV               = _mm_macc_pd(vfeps,Fp,Y);
275             vvdw12           = _mm_mul_pd(c12_00,VV);
276             FF               = _mm_macc_pd(vfeps,_mm_macc_pd(twovfeps,H,G),Fp);
277             fvdw12           = _mm_mul_pd(c12_00,FF);
278             vvdw             = _mm_add_pd(vvdw12,vvdw6);
279             fvdw             = _mm_xor_pd(signbit,_mm_mul_pd(_mm_add_pd(fvdw6,fvdw12),_mm_mul_pd(vftabscale,rinv00)));
280
281             /* Update potential sum for this i atom from the interaction with this j atom. */
282             velecsum         = _mm_add_pd(velecsum,velec);
283             vvdwsum          = _mm_add_pd(vvdwsum,vvdw);
284
285             fscal            = _mm_add_pd(felec,fvdw);
286
287             /* Update vectorial force */
288             fix0             = _mm_macc_pd(dx00,fscal,fix0);
289             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
290             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
291             
292             fjx0             = _mm_macc_pd(dx00,fscal,fjx0);
293             fjy0             = _mm_macc_pd(dy00,fscal,fjy0);
294             fjz0             = _mm_macc_pd(dz00,fscal,fjz0);
295
296             /**************************
297              * CALCULATE INTERACTIONS *
298              **************************/
299
300             /* Compute parameters for interactions between i and j atoms */
301             qq10             = _mm_mul_pd(iq1,jq0);
302
303             /* REACTION-FIELD ELECTROSTATICS */
304             velec            = _mm_mul_pd(qq10,_mm_sub_pd(_mm_macc_pd(krf,rsq10,rinv10),crf));
305             felec            = _mm_mul_pd(qq10,_mm_msub_pd(rinv10,rinvsq10,krf2));
306
307             /* Update potential sum for this i atom from the interaction with this j atom. */
308             velecsum         = _mm_add_pd(velecsum,velec);
309
310             fscal            = felec;
311
312             /* Update vectorial force */
313             fix1             = _mm_macc_pd(dx10,fscal,fix1);
314             fiy1             = _mm_macc_pd(dy10,fscal,fiy1);
315             fiz1             = _mm_macc_pd(dz10,fscal,fiz1);
316             
317             fjx0             = _mm_macc_pd(dx10,fscal,fjx0);
318             fjy0             = _mm_macc_pd(dy10,fscal,fjy0);
319             fjz0             = _mm_macc_pd(dz10,fscal,fjz0);
320
321             /**************************
322              * CALCULATE INTERACTIONS *
323              **************************/
324
325             /* Compute parameters for interactions between i and j atoms */
326             qq20             = _mm_mul_pd(iq2,jq0);
327
328             /* REACTION-FIELD ELECTROSTATICS */
329             velec            = _mm_mul_pd(qq20,_mm_sub_pd(_mm_macc_pd(krf,rsq20,rinv20),crf));
330             felec            = _mm_mul_pd(qq20,_mm_msub_pd(rinv20,rinvsq20,krf2));
331
332             /* Update potential sum for this i atom from the interaction with this j atom. */
333             velecsum         = _mm_add_pd(velecsum,velec);
334
335             fscal            = felec;
336
337             /* Update vectorial force */
338             fix2             = _mm_macc_pd(dx20,fscal,fix2);
339             fiy2             = _mm_macc_pd(dy20,fscal,fiy2);
340             fiz2             = _mm_macc_pd(dz20,fscal,fiz2);
341             
342             fjx0             = _mm_macc_pd(dx20,fscal,fjx0);
343             fjy0             = _mm_macc_pd(dy20,fscal,fjy0);
344             fjz0             = _mm_macc_pd(dz20,fscal,fjz0);
345
346             gmx_mm_decrement_1rvec_2ptr_swizzle_pd(f+j_coord_offsetA,f+j_coord_offsetB,fjx0,fjy0,fjz0);
347
348             /* Inner loop uses 143 flops */
349         }
350
351         if(jidx<j_index_end)
352         {
353
354             jnrA             = jjnr[jidx];
355             j_coord_offsetA  = DIM*jnrA;
356
357             /* load j atom coordinates */
358             gmx_mm_load_1rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
359                                               &jx0,&jy0,&jz0);
360
361             /* Calculate displacement vector */
362             dx00             = _mm_sub_pd(ix0,jx0);
363             dy00             = _mm_sub_pd(iy0,jy0);
364             dz00             = _mm_sub_pd(iz0,jz0);
365             dx10             = _mm_sub_pd(ix1,jx0);
366             dy10             = _mm_sub_pd(iy1,jy0);
367             dz10             = _mm_sub_pd(iz1,jz0);
368             dx20             = _mm_sub_pd(ix2,jx0);
369             dy20             = _mm_sub_pd(iy2,jy0);
370             dz20             = _mm_sub_pd(iz2,jz0);
371
372             /* Calculate squared distance and things based on it */
373             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
374             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
375             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
376
377             rinv00           = gmx_mm_invsqrt_pd(rsq00);
378             rinv10           = gmx_mm_invsqrt_pd(rsq10);
379             rinv20           = gmx_mm_invsqrt_pd(rsq20);
380
381             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
382             rinvsq10         = _mm_mul_pd(rinv10,rinv10);
383             rinvsq20         = _mm_mul_pd(rinv20,rinv20);
384
385             /* Load parameters for j particles */
386             jq0              = _mm_load_sd(charge+jnrA+0);
387             vdwjidx0A        = 2*vdwtype[jnrA+0];
388
389             fjx0             = _mm_setzero_pd();
390             fjy0             = _mm_setzero_pd();
391             fjz0             = _mm_setzero_pd();
392
393             /**************************
394              * CALCULATE INTERACTIONS *
395              **************************/
396
397             r00              = _mm_mul_pd(rsq00,rinv00);
398
399             /* Compute parameters for interactions between i and j atoms */
400             qq00             = _mm_mul_pd(iq0,jq0);
401             gmx_mm_load_1pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,&c6_00,&c12_00);
402
403             /* Calculate table index by multiplying r with table scale and truncate to integer */
404             rt               = _mm_mul_pd(r00,vftabscale);
405             vfitab           = _mm_cvttpd_epi32(rt);
406 #ifdef __XOP__
407             vfeps            = _mm_frcz_pd(rt);
408 #else
409             vfeps            = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
410 #endif
411             twovfeps         = _mm_add_pd(vfeps,vfeps);
412             vfitab           = _mm_slli_epi32(vfitab,3);
413
414             /* REACTION-FIELD ELECTROSTATICS */
415             velec            = _mm_mul_pd(qq00,_mm_sub_pd(_mm_macc_pd(krf,rsq00,rinv00),crf));
416             felec            = _mm_mul_pd(qq00,_mm_msub_pd(rinv00,rinvsq00,krf2));
417
418             /* CUBIC SPLINE TABLE DISPERSION */
419             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
420             F                = _mm_setzero_pd();
421             GMX_MM_TRANSPOSE2_PD(Y,F);
422             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
423             H                = _mm_setzero_pd();
424             GMX_MM_TRANSPOSE2_PD(G,H);
425             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(H,vfeps,G),F);
426             VV               = _mm_macc_pd(vfeps,Fp,Y);
427             vvdw6            = _mm_mul_pd(c6_00,VV);
428             FF               = _mm_macc_pd(vfeps,_mm_macc_pd(twovfeps,H,G),Fp);
429             fvdw6            = _mm_mul_pd(c6_00,FF);
430
431             /* CUBIC SPLINE TABLE REPULSION */
432             vfitab           = _mm_add_epi32(vfitab,ifour);
433             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
434             F                = _mm_setzero_pd();
435             GMX_MM_TRANSPOSE2_PD(Y,F);
436             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
437             H                = _mm_setzero_pd();
438             GMX_MM_TRANSPOSE2_PD(G,H);
439             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(H,vfeps,G),F);
440             VV               = _mm_macc_pd(vfeps,Fp,Y);
441             vvdw12           = _mm_mul_pd(c12_00,VV);
442             FF               = _mm_macc_pd(vfeps,_mm_macc_pd(twovfeps,H,G),Fp);
443             fvdw12           = _mm_mul_pd(c12_00,FF);
444             vvdw             = _mm_add_pd(vvdw12,vvdw6);
445             fvdw             = _mm_xor_pd(signbit,_mm_mul_pd(_mm_add_pd(fvdw6,fvdw12),_mm_mul_pd(vftabscale,rinv00)));
446
447             /* Update potential sum for this i atom from the interaction with this j atom. */
448             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
449             velecsum         = _mm_add_pd(velecsum,velec);
450             vvdw             = _mm_unpacklo_pd(vvdw,_mm_setzero_pd());
451             vvdwsum          = _mm_add_pd(vvdwsum,vvdw);
452
453             fscal            = _mm_add_pd(felec,fvdw);
454
455             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
456
457             /* Update vectorial force */
458             fix0             = _mm_macc_pd(dx00,fscal,fix0);
459             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
460             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
461             
462             fjx0             = _mm_macc_pd(dx00,fscal,fjx0);
463             fjy0             = _mm_macc_pd(dy00,fscal,fjy0);
464             fjz0             = _mm_macc_pd(dz00,fscal,fjz0);
465
466             /**************************
467              * CALCULATE INTERACTIONS *
468              **************************/
469
470             /* Compute parameters for interactions between i and j atoms */
471             qq10             = _mm_mul_pd(iq1,jq0);
472
473             /* REACTION-FIELD ELECTROSTATICS */
474             velec            = _mm_mul_pd(qq10,_mm_sub_pd(_mm_macc_pd(krf,rsq10,rinv10),crf));
475             felec            = _mm_mul_pd(qq10,_mm_msub_pd(rinv10,rinvsq10,krf2));
476
477             /* Update potential sum for this i atom from the interaction with this j atom. */
478             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
479             velecsum         = _mm_add_pd(velecsum,velec);
480
481             fscal            = felec;
482
483             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
484
485             /* Update vectorial force */
486             fix1             = _mm_macc_pd(dx10,fscal,fix1);
487             fiy1             = _mm_macc_pd(dy10,fscal,fiy1);
488             fiz1             = _mm_macc_pd(dz10,fscal,fiz1);
489             
490             fjx0             = _mm_macc_pd(dx10,fscal,fjx0);
491             fjy0             = _mm_macc_pd(dy10,fscal,fjy0);
492             fjz0             = _mm_macc_pd(dz10,fscal,fjz0);
493
494             /**************************
495              * CALCULATE INTERACTIONS *
496              **************************/
497
498             /* Compute parameters for interactions between i and j atoms */
499             qq20             = _mm_mul_pd(iq2,jq0);
500
501             /* REACTION-FIELD ELECTROSTATICS */
502             velec            = _mm_mul_pd(qq20,_mm_sub_pd(_mm_macc_pd(krf,rsq20,rinv20),crf));
503             felec            = _mm_mul_pd(qq20,_mm_msub_pd(rinv20,rinvsq20,krf2));
504
505             /* Update potential sum for this i atom from the interaction with this j atom. */
506             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
507             velecsum         = _mm_add_pd(velecsum,velec);
508
509             fscal            = felec;
510
511             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
512
513             /* Update vectorial force */
514             fix2             = _mm_macc_pd(dx20,fscal,fix2);
515             fiy2             = _mm_macc_pd(dy20,fscal,fiy2);
516             fiz2             = _mm_macc_pd(dz20,fscal,fiz2);
517             
518             fjx0             = _mm_macc_pd(dx20,fscal,fjx0);
519             fjy0             = _mm_macc_pd(dy20,fscal,fjy0);
520             fjz0             = _mm_macc_pd(dz20,fscal,fjz0);
521
522             gmx_mm_decrement_1rvec_1ptr_swizzle_pd(f+j_coord_offsetA,fjx0,fjy0,fjz0);
523
524             /* Inner loop uses 143 flops */
525         }
526
527         /* End of innermost loop */
528
529         gmx_mm_update_iforce_3atom_swizzle_pd(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,
530                                               f+i_coord_offset,fshift+i_shift_offset);
531
532         ggid                        = gid[iidx];
533         /* Update potential energies */
534         gmx_mm_update_1pot_pd(velecsum,kernel_data->energygrp_elec+ggid);
535         gmx_mm_update_1pot_pd(vvdwsum,kernel_data->energygrp_vdw+ggid);
536
537         /* Increment number of inner iterations */
538         inneriter                  += j_index_end - j_index_start;
539
540         /* Outer loop uses 20 flops */
541     }
542
543     /* Increment number of outer iterations */
544     outeriter        += nri;
545
546     /* Update outer/inner flops */
547
548     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W3_VF,outeriter*20 + inneriter*143);
549 }
550 /*
551  * Gromacs nonbonded kernel:   nb_kernel_ElecRF_VdwCSTab_GeomW3P1_F_avx_128_fma_double
552  * Electrostatics interaction: ReactionField
553  * VdW interaction:            CubicSplineTable
554  * Geometry:                   Water3-Particle
555  * Calculate force/pot:        Force
556  */
557 void
558 nb_kernel_ElecRF_VdwCSTab_GeomW3P1_F_avx_128_fma_double
559                     (t_nblist                    * gmx_restrict       nlist,
560                      rvec                        * gmx_restrict          xx,
561                      rvec                        * gmx_restrict          ff,
562                      t_forcerec                  * gmx_restrict          fr,
563                      t_mdatoms                   * gmx_restrict     mdatoms,
564                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
565                      t_nrnb                      * gmx_restrict        nrnb)
566 {
567     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
568      * just 0 for non-waters.
569      * Suffixes A,B refer to j loop unrolling done with SSE double precision, e.g. for the two different
570      * jnr indices corresponding to data put in the four positions in the SIMD register.
571      */
572     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
573     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
574     int              jnrA,jnrB;
575     int              j_coord_offsetA,j_coord_offsetB;
576     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
577     real             rcutoff_scalar;
578     real             *shiftvec,*fshift,*x,*f;
579     __m128d          tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
580     int              vdwioffset0;
581     __m128d          ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
582     int              vdwioffset1;
583     __m128d          ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
584     int              vdwioffset2;
585     __m128d          ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
586     int              vdwjidx0A,vdwjidx0B;
587     __m128d          jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
588     __m128d          dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
589     __m128d          dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
590     __m128d          dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
591     __m128d          velec,felec,velecsum,facel,crf,krf,krf2;
592     real             *charge;
593     int              nvdwtype;
594     __m128d          rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
595     int              *vdwtype;
596     real             *vdwparam;
597     __m128d          one_sixth   = _mm_set1_pd(1.0/6.0);
598     __m128d          one_twelfth = _mm_set1_pd(1.0/12.0);
599     __m128i          vfitab;
600     __m128i          ifour       = _mm_set1_epi32(4);
601     __m128d          rt,vfeps,vftabscale,Y,F,G,H,Heps,Fp,VV,FF,twovfeps;
602     real             *vftab;
603     __m128d          dummy_mask,cutoff_mask;
604     __m128d          signbit   = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
605     __m128d          one     = _mm_set1_pd(1.0);
606     __m128d          two     = _mm_set1_pd(2.0);
607     x                = xx[0];
608     f                = ff[0];
609
610     nri              = nlist->nri;
611     iinr             = nlist->iinr;
612     jindex           = nlist->jindex;
613     jjnr             = nlist->jjnr;
614     shiftidx         = nlist->shift;
615     gid              = nlist->gid;
616     shiftvec         = fr->shift_vec[0];
617     fshift           = fr->fshift[0];
618     facel            = _mm_set1_pd(fr->epsfac);
619     charge           = mdatoms->chargeA;
620     krf              = _mm_set1_pd(fr->ic->k_rf);
621     krf2             = _mm_set1_pd(fr->ic->k_rf*2.0);
622     crf              = _mm_set1_pd(fr->ic->c_rf);
623     nvdwtype         = fr->ntype;
624     vdwparam         = fr->nbfp;
625     vdwtype          = mdatoms->typeA;
626
627     vftab            = kernel_data->table_vdw->data;
628     vftabscale       = _mm_set1_pd(kernel_data->table_vdw->scale);
629
630     /* Setup water-specific parameters */
631     inr              = nlist->iinr[0];
632     iq0              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+0]));
633     iq1              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+1]));
634     iq2              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+2]));
635     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
636
637     /* Avoid stupid compiler warnings */
638     jnrA = jnrB = 0;
639     j_coord_offsetA = 0;
640     j_coord_offsetB = 0;
641
642     outeriter        = 0;
643     inneriter        = 0;
644
645     /* Start outer loop over neighborlists */
646     for(iidx=0; iidx<nri; iidx++)
647     {
648         /* Load shift vector for this list */
649         i_shift_offset   = DIM*shiftidx[iidx];
650
651         /* Load limits for loop over neighbors */
652         j_index_start    = jindex[iidx];
653         j_index_end      = jindex[iidx+1];
654
655         /* Get outer coordinate index */
656         inr              = iinr[iidx];
657         i_coord_offset   = DIM*inr;
658
659         /* Load i particle coords and add shift vector */
660         gmx_mm_load_shift_and_3rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,
661                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2);
662
663         fix0             = _mm_setzero_pd();
664         fiy0             = _mm_setzero_pd();
665         fiz0             = _mm_setzero_pd();
666         fix1             = _mm_setzero_pd();
667         fiy1             = _mm_setzero_pd();
668         fiz1             = _mm_setzero_pd();
669         fix2             = _mm_setzero_pd();
670         fiy2             = _mm_setzero_pd();
671         fiz2             = _mm_setzero_pd();
672
673         /* Start inner kernel loop */
674         for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
675         {
676
677             /* Get j neighbor index, and coordinate index */
678             jnrA             = jjnr[jidx];
679             jnrB             = jjnr[jidx+1];
680             j_coord_offsetA  = DIM*jnrA;
681             j_coord_offsetB  = DIM*jnrB;
682
683             /* load j atom coordinates */
684             gmx_mm_load_1rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
685                                               &jx0,&jy0,&jz0);
686
687             /* Calculate displacement vector */
688             dx00             = _mm_sub_pd(ix0,jx0);
689             dy00             = _mm_sub_pd(iy0,jy0);
690             dz00             = _mm_sub_pd(iz0,jz0);
691             dx10             = _mm_sub_pd(ix1,jx0);
692             dy10             = _mm_sub_pd(iy1,jy0);
693             dz10             = _mm_sub_pd(iz1,jz0);
694             dx20             = _mm_sub_pd(ix2,jx0);
695             dy20             = _mm_sub_pd(iy2,jy0);
696             dz20             = _mm_sub_pd(iz2,jz0);
697
698             /* Calculate squared distance and things based on it */
699             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
700             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
701             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
702
703             rinv00           = gmx_mm_invsqrt_pd(rsq00);
704             rinv10           = gmx_mm_invsqrt_pd(rsq10);
705             rinv20           = gmx_mm_invsqrt_pd(rsq20);
706
707             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
708             rinvsq10         = _mm_mul_pd(rinv10,rinv10);
709             rinvsq20         = _mm_mul_pd(rinv20,rinv20);
710
711             /* Load parameters for j particles */
712             jq0              = gmx_mm_load_2real_swizzle_pd(charge+jnrA+0,charge+jnrB+0);
713             vdwjidx0A        = 2*vdwtype[jnrA+0];
714             vdwjidx0B        = 2*vdwtype[jnrB+0];
715
716             fjx0             = _mm_setzero_pd();
717             fjy0             = _mm_setzero_pd();
718             fjz0             = _mm_setzero_pd();
719
720             /**************************
721              * CALCULATE INTERACTIONS *
722              **************************/
723
724             r00              = _mm_mul_pd(rsq00,rinv00);
725
726             /* Compute parameters for interactions between i and j atoms */
727             qq00             = _mm_mul_pd(iq0,jq0);
728             gmx_mm_load_2pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,
729                                          vdwparam+vdwioffset0+vdwjidx0B,&c6_00,&c12_00);
730
731             /* Calculate table index by multiplying r with table scale and truncate to integer */
732             rt               = _mm_mul_pd(r00,vftabscale);
733             vfitab           = _mm_cvttpd_epi32(rt);
734 #ifdef __XOP__
735             vfeps            = _mm_frcz_pd(rt);
736 #else
737             vfeps            = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
738 #endif
739             twovfeps         = _mm_add_pd(vfeps,vfeps);
740             vfitab           = _mm_slli_epi32(vfitab,3);
741
742             /* REACTION-FIELD ELECTROSTATICS */
743             felec            = _mm_mul_pd(qq00,_mm_msub_pd(rinv00,rinvsq00,krf2));
744
745             /* CUBIC SPLINE TABLE DISPERSION */
746             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
747             F                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) );
748             GMX_MM_TRANSPOSE2_PD(Y,F);
749             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
750             H                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) +2);
751             GMX_MM_TRANSPOSE2_PD(G,H);
752             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(H,vfeps,G),F);
753             FF               = _mm_macc_pd(vfeps,_mm_macc_pd(twovfeps,H,G),Fp);
754             fvdw6            = _mm_mul_pd(c6_00,FF);
755
756             /* CUBIC SPLINE TABLE REPULSION */
757             vfitab           = _mm_add_epi32(vfitab,ifour);
758             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
759             F                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) );
760             GMX_MM_TRANSPOSE2_PD(Y,F);
761             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
762             H                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) +2);
763             GMX_MM_TRANSPOSE2_PD(G,H);
764             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(H,vfeps,G),F);
765             FF               = _mm_macc_pd(vfeps,_mm_macc_pd(twovfeps,H,G),Fp);
766             fvdw12           = _mm_mul_pd(c12_00,FF);
767             fvdw             = _mm_xor_pd(signbit,_mm_mul_pd(_mm_add_pd(fvdw6,fvdw12),_mm_mul_pd(vftabscale,rinv00)));
768
769             fscal            = _mm_add_pd(felec,fvdw);
770
771             /* Update vectorial force */
772             fix0             = _mm_macc_pd(dx00,fscal,fix0);
773             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
774             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
775             
776             fjx0             = _mm_macc_pd(dx00,fscal,fjx0);
777             fjy0             = _mm_macc_pd(dy00,fscal,fjy0);
778             fjz0             = _mm_macc_pd(dz00,fscal,fjz0);
779
780             /**************************
781              * CALCULATE INTERACTIONS *
782              **************************/
783
784             /* Compute parameters for interactions between i and j atoms */
785             qq10             = _mm_mul_pd(iq1,jq0);
786
787             /* REACTION-FIELD ELECTROSTATICS */
788             felec            = _mm_mul_pd(qq10,_mm_msub_pd(rinv10,rinvsq10,krf2));
789
790             fscal            = felec;
791
792             /* Update vectorial force */
793             fix1             = _mm_macc_pd(dx10,fscal,fix1);
794             fiy1             = _mm_macc_pd(dy10,fscal,fiy1);
795             fiz1             = _mm_macc_pd(dz10,fscal,fiz1);
796             
797             fjx0             = _mm_macc_pd(dx10,fscal,fjx0);
798             fjy0             = _mm_macc_pd(dy10,fscal,fjy0);
799             fjz0             = _mm_macc_pd(dz10,fscal,fjz0);
800
801             /**************************
802              * CALCULATE INTERACTIONS *
803              **************************/
804
805             /* Compute parameters for interactions between i and j atoms */
806             qq20             = _mm_mul_pd(iq2,jq0);
807
808             /* REACTION-FIELD ELECTROSTATICS */
809             felec            = _mm_mul_pd(qq20,_mm_msub_pd(rinv20,rinvsq20,krf2));
810
811             fscal            = felec;
812
813             /* Update vectorial force */
814             fix2             = _mm_macc_pd(dx20,fscal,fix2);
815             fiy2             = _mm_macc_pd(dy20,fscal,fiy2);
816             fiz2             = _mm_macc_pd(dz20,fscal,fiz2);
817             
818             fjx0             = _mm_macc_pd(dx20,fscal,fjx0);
819             fjy0             = _mm_macc_pd(dy20,fscal,fjy0);
820             fjz0             = _mm_macc_pd(dz20,fscal,fjz0);
821
822             gmx_mm_decrement_1rvec_2ptr_swizzle_pd(f+j_coord_offsetA,f+j_coord_offsetB,fjx0,fjy0,fjz0);
823
824             /* Inner loop uses 120 flops */
825         }
826
827         if(jidx<j_index_end)
828         {
829
830             jnrA             = jjnr[jidx];
831             j_coord_offsetA  = DIM*jnrA;
832
833             /* load j atom coordinates */
834             gmx_mm_load_1rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
835                                               &jx0,&jy0,&jz0);
836
837             /* Calculate displacement vector */
838             dx00             = _mm_sub_pd(ix0,jx0);
839             dy00             = _mm_sub_pd(iy0,jy0);
840             dz00             = _mm_sub_pd(iz0,jz0);
841             dx10             = _mm_sub_pd(ix1,jx0);
842             dy10             = _mm_sub_pd(iy1,jy0);
843             dz10             = _mm_sub_pd(iz1,jz0);
844             dx20             = _mm_sub_pd(ix2,jx0);
845             dy20             = _mm_sub_pd(iy2,jy0);
846             dz20             = _mm_sub_pd(iz2,jz0);
847
848             /* Calculate squared distance and things based on it */
849             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
850             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
851             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
852
853             rinv00           = gmx_mm_invsqrt_pd(rsq00);
854             rinv10           = gmx_mm_invsqrt_pd(rsq10);
855             rinv20           = gmx_mm_invsqrt_pd(rsq20);
856
857             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
858             rinvsq10         = _mm_mul_pd(rinv10,rinv10);
859             rinvsq20         = _mm_mul_pd(rinv20,rinv20);
860
861             /* Load parameters for j particles */
862             jq0              = _mm_load_sd(charge+jnrA+0);
863             vdwjidx0A        = 2*vdwtype[jnrA+0];
864
865             fjx0             = _mm_setzero_pd();
866             fjy0             = _mm_setzero_pd();
867             fjz0             = _mm_setzero_pd();
868
869             /**************************
870              * CALCULATE INTERACTIONS *
871              **************************/
872
873             r00              = _mm_mul_pd(rsq00,rinv00);
874
875             /* Compute parameters for interactions between i and j atoms */
876             qq00             = _mm_mul_pd(iq0,jq0);
877             gmx_mm_load_1pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,&c6_00,&c12_00);
878
879             /* Calculate table index by multiplying r with table scale and truncate to integer */
880             rt               = _mm_mul_pd(r00,vftabscale);
881             vfitab           = _mm_cvttpd_epi32(rt);
882 #ifdef __XOP__
883             vfeps            = _mm_frcz_pd(rt);
884 #else
885             vfeps            = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
886 #endif
887             twovfeps         = _mm_add_pd(vfeps,vfeps);
888             vfitab           = _mm_slli_epi32(vfitab,3);
889
890             /* REACTION-FIELD ELECTROSTATICS */
891             felec            = _mm_mul_pd(qq00,_mm_msub_pd(rinv00,rinvsq00,krf2));
892
893             /* CUBIC SPLINE TABLE DISPERSION */
894             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
895             F                = _mm_setzero_pd();
896             GMX_MM_TRANSPOSE2_PD(Y,F);
897             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
898             H                = _mm_setzero_pd();
899             GMX_MM_TRANSPOSE2_PD(G,H);
900             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(H,vfeps,G),F);
901             FF               = _mm_macc_pd(vfeps,_mm_macc_pd(twovfeps,H,G),Fp);
902             fvdw6            = _mm_mul_pd(c6_00,FF);
903
904             /* CUBIC SPLINE TABLE REPULSION */
905             vfitab           = _mm_add_epi32(vfitab,ifour);
906             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
907             F                = _mm_setzero_pd();
908             GMX_MM_TRANSPOSE2_PD(Y,F);
909             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
910             H                = _mm_setzero_pd();
911             GMX_MM_TRANSPOSE2_PD(G,H);
912             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(H,vfeps,G),F);
913             FF               = _mm_macc_pd(vfeps,_mm_macc_pd(twovfeps,H,G),Fp);
914             fvdw12           = _mm_mul_pd(c12_00,FF);
915             fvdw             = _mm_xor_pd(signbit,_mm_mul_pd(_mm_add_pd(fvdw6,fvdw12),_mm_mul_pd(vftabscale,rinv00)));
916
917             fscal            = _mm_add_pd(felec,fvdw);
918
919             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
920
921             /* Update vectorial force */
922             fix0             = _mm_macc_pd(dx00,fscal,fix0);
923             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
924             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
925             
926             fjx0             = _mm_macc_pd(dx00,fscal,fjx0);
927             fjy0             = _mm_macc_pd(dy00,fscal,fjy0);
928             fjz0             = _mm_macc_pd(dz00,fscal,fjz0);
929
930             /**************************
931              * CALCULATE INTERACTIONS *
932              **************************/
933
934             /* Compute parameters for interactions between i and j atoms */
935             qq10             = _mm_mul_pd(iq1,jq0);
936
937             /* REACTION-FIELD ELECTROSTATICS */
938             felec            = _mm_mul_pd(qq10,_mm_msub_pd(rinv10,rinvsq10,krf2));
939
940             fscal            = felec;
941
942             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
943
944             /* Update vectorial force */
945             fix1             = _mm_macc_pd(dx10,fscal,fix1);
946             fiy1             = _mm_macc_pd(dy10,fscal,fiy1);
947             fiz1             = _mm_macc_pd(dz10,fscal,fiz1);
948             
949             fjx0             = _mm_macc_pd(dx10,fscal,fjx0);
950             fjy0             = _mm_macc_pd(dy10,fscal,fjy0);
951             fjz0             = _mm_macc_pd(dz10,fscal,fjz0);
952
953             /**************************
954              * CALCULATE INTERACTIONS *
955              **************************/
956
957             /* Compute parameters for interactions between i and j atoms */
958             qq20             = _mm_mul_pd(iq2,jq0);
959
960             /* REACTION-FIELD ELECTROSTATICS */
961             felec            = _mm_mul_pd(qq20,_mm_msub_pd(rinv20,rinvsq20,krf2));
962
963             fscal            = felec;
964
965             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
966
967             /* Update vectorial force */
968             fix2             = _mm_macc_pd(dx20,fscal,fix2);
969             fiy2             = _mm_macc_pd(dy20,fscal,fiy2);
970             fiz2             = _mm_macc_pd(dz20,fscal,fiz2);
971             
972             fjx0             = _mm_macc_pd(dx20,fscal,fjx0);
973             fjy0             = _mm_macc_pd(dy20,fscal,fjy0);
974             fjz0             = _mm_macc_pd(dz20,fscal,fjz0);
975
976             gmx_mm_decrement_1rvec_1ptr_swizzle_pd(f+j_coord_offsetA,fjx0,fjy0,fjz0);
977
978             /* Inner loop uses 120 flops */
979         }
980
981         /* End of innermost loop */
982
983         gmx_mm_update_iforce_3atom_swizzle_pd(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,
984                                               f+i_coord_offset,fshift+i_shift_offset);
985
986         /* Increment number of inner iterations */
987         inneriter                  += j_index_end - j_index_start;
988
989         /* Outer loop uses 18 flops */
990     }
991
992     /* Increment number of outer iterations */
993     outeriter        += nri;
994
995     /* Update outer/inner flops */
996
997     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W3_F,outeriter*18 + inneriter*120);
998 }