07584802d29a76dab72752531070ebf0160e9b77
[alexxy/gromacs.git] / src / gromacs / gmxlib / nonbonded / nb_kernel_avx_128_fma_double / nb_kernel_ElecCSTab_VdwNone_GeomW4P1_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 #ifdef HAVE_CONFIG_H
39 #include <config.h>
40 #endif
41
42 #include <math.h>
43
44 #include "../nb_kernel.h"
45 #include "types/simple.h"
46 #include "gromacs/math/vec.h"
47 #include "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_ElecCSTab_VdwNone_GeomW4P1_VF_avx_128_fma_double
54  * Electrostatics interaction: CubicSplineTable
55  * VdW interaction:            None
56  * Geometry:                   Water4-Particle
57  * Calculate force/pot:        PotentialAndForce
58  */
59 void
60 nb_kernel_ElecCSTab_VdwNone_GeomW4P1_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              vdwioffset1;
83     __m128d          ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
84     int              vdwioffset2;
85     __m128d          ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
86     int              vdwioffset3;
87     __m128d          ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
88     int              vdwjidx0A,vdwjidx0B;
89     __m128d          jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
90     __m128d          dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
91     __m128d          dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
92     __m128d          dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
93     __m128d          velec,felec,velecsum,facel,crf,krf,krf2;
94     real             *charge;
95     __m128i          vfitab;
96     __m128i          ifour       = _mm_set1_epi32(4);
97     __m128d          rt,vfeps,vftabscale,Y,F,G,H,Heps,Fp,VV,FF,twovfeps;
98     real             *vftab;
99     __m128d          dummy_mask,cutoff_mask;
100     __m128d          signbit   = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
101     __m128d          one     = _mm_set1_pd(1.0);
102     __m128d          two     = _mm_set1_pd(2.0);
103     x                = xx[0];
104     f                = ff[0];
105
106     nri              = nlist->nri;
107     iinr             = nlist->iinr;
108     jindex           = nlist->jindex;
109     jjnr             = nlist->jjnr;
110     shiftidx         = nlist->shift;
111     gid              = nlist->gid;
112     shiftvec         = fr->shift_vec[0];
113     fshift           = fr->fshift[0];
114     facel            = _mm_set1_pd(fr->epsfac);
115     charge           = mdatoms->chargeA;
116
117     vftab            = kernel_data->table_elec->data;
118     vftabscale       = _mm_set1_pd(kernel_data->table_elec->scale);
119
120     /* Setup water-specific parameters */
121     inr              = nlist->iinr[0];
122     iq1              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+1]));
123     iq2              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+2]));
124     iq3              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+3]));
125
126     /* Avoid stupid compiler warnings */
127     jnrA = jnrB = 0;
128     j_coord_offsetA = 0;
129     j_coord_offsetB = 0;
130
131     outeriter        = 0;
132     inneriter        = 0;
133
134     /* Start outer loop over neighborlists */
135     for(iidx=0; iidx<nri; iidx++)
136     {
137         /* Load shift vector for this list */
138         i_shift_offset   = DIM*shiftidx[iidx];
139
140         /* Load limits for loop over neighbors */
141         j_index_start    = jindex[iidx];
142         j_index_end      = jindex[iidx+1];
143
144         /* Get outer coordinate index */
145         inr              = iinr[iidx];
146         i_coord_offset   = DIM*inr;
147
148         /* Load i particle coords and add shift vector */
149         gmx_mm_load_shift_and_3rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset+DIM,
150                                                  &ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
151
152         fix1             = _mm_setzero_pd();
153         fiy1             = _mm_setzero_pd();
154         fiz1             = _mm_setzero_pd();
155         fix2             = _mm_setzero_pd();
156         fiy2             = _mm_setzero_pd();
157         fiz2             = _mm_setzero_pd();
158         fix3             = _mm_setzero_pd();
159         fiy3             = _mm_setzero_pd();
160         fiz3             = _mm_setzero_pd();
161
162         /* Reset potential sums */
163         velecsum         = _mm_setzero_pd();
164
165         /* Start inner kernel loop */
166         for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
167         {
168
169             /* Get j neighbor index, and coordinate index */
170             jnrA             = jjnr[jidx];
171             jnrB             = jjnr[jidx+1];
172             j_coord_offsetA  = DIM*jnrA;
173             j_coord_offsetB  = DIM*jnrB;
174
175             /* load j atom coordinates */
176             gmx_mm_load_1rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
177                                               &jx0,&jy0,&jz0);
178
179             /* Calculate displacement vector */
180             dx10             = _mm_sub_pd(ix1,jx0);
181             dy10             = _mm_sub_pd(iy1,jy0);
182             dz10             = _mm_sub_pd(iz1,jz0);
183             dx20             = _mm_sub_pd(ix2,jx0);
184             dy20             = _mm_sub_pd(iy2,jy0);
185             dz20             = _mm_sub_pd(iz2,jz0);
186             dx30             = _mm_sub_pd(ix3,jx0);
187             dy30             = _mm_sub_pd(iy3,jy0);
188             dz30             = _mm_sub_pd(iz3,jz0);
189
190             /* Calculate squared distance and things based on it */
191             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
192             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
193             rsq30            = gmx_mm_calc_rsq_pd(dx30,dy30,dz30);
194
195             rinv10           = gmx_mm_invsqrt_pd(rsq10);
196             rinv20           = gmx_mm_invsqrt_pd(rsq20);
197             rinv30           = gmx_mm_invsqrt_pd(rsq30);
198
199             /* Load parameters for j particles */
200             jq0              = gmx_mm_load_2real_swizzle_pd(charge+jnrA+0,charge+jnrB+0);
201
202             fjx0             = _mm_setzero_pd();
203             fjy0             = _mm_setzero_pd();
204             fjz0             = _mm_setzero_pd();
205
206             /**************************
207              * CALCULATE INTERACTIONS *
208              **************************/
209
210             r10              = _mm_mul_pd(rsq10,rinv10);
211
212             /* Compute parameters for interactions between i and j atoms */
213             qq10             = _mm_mul_pd(iq1,jq0);
214
215             /* Calculate table index by multiplying r with table scale and truncate to integer */
216             rt               = _mm_mul_pd(r10,vftabscale);
217             vfitab           = _mm_cvttpd_epi32(rt);
218 #ifdef __XOP__
219             vfeps            = _mm_frcz_pd(rt);
220 #else
221             vfeps            = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
222 #endif
223             twovfeps         = _mm_add_pd(vfeps,vfeps);
224             vfitab           = _mm_slli_epi32(vfitab,2);
225
226             /* CUBIC SPLINE TABLE ELECTROSTATICS */
227             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
228             F                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) );
229             GMX_MM_TRANSPOSE2_PD(Y,F);
230             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
231             H                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) +2);
232             GMX_MM_TRANSPOSE2_PD(G,H);
233             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(vfeps,H,G),F);
234             VV               = _mm_macc_pd(vfeps,Fp,Y);
235             velec            = _mm_mul_pd(qq10,VV);
236             FF               = _mm_macc_pd(_mm_macc_pd(twovfeps,H,G),vfeps,Fp);
237             felec            = _mm_xor_pd(signbit,_mm_mul_pd(_mm_mul_pd(qq10,FF),_mm_mul_pd(vftabscale,rinv10)));
238
239             /* Update potential sum for this i atom from the interaction with this j atom. */
240             velecsum         = _mm_add_pd(velecsum,velec);
241
242             fscal            = felec;
243
244             /* Update vectorial force */
245             fix1             = _mm_macc_pd(dx10,fscal,fix1);
246             fiy1             = _mm_macc_pd(dy10,fscal,fiy1);
247             fiz1             = _mm_macc_pd(dz10,fscal,fiz1);
248             
249             fjx0             = _mm_macc_pd(dx10,fscal,fjx0);
250             fjy0             = _mm_macc_pd(dy10,fscal,fjy0);
251             fjz0             = _mm_macc_pd(dz10,fscal,fjz0);
252
253             /**************************
254              * CALCULATE INTERACTIONS *
255              **************************/
256
257             r20              = _mm_mul_pd(rsq20,rinv20);
258
259             /* Compute parameters for interactions between i and j atoms */
260             qq20             = _mm_mul_pd(iq2,jq0);
261
262             /* Calculate table index by multiplying r with table scale and truncate to integer */
263             rt               = _mm_mul_pd(r20,vftabscale);
264             vfitab           = _mm_cvttpd_epi32(rt);
265 #ifdef __XOP__
266             vfeps            = _mm_frcz_pd(rt);
267 #else
268             vfeps            = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
269 #endif
270             twovfeps         = _mm_add_pd(vfeps,vfeps);
271             vfitab           = _mm_slli_epi32(vfitab,2);
272
273             /* CUBIC SPLINE TABLE ELECTROSTATICS */
274             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
275             F                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) );
276             GMX_MM_TRANSPOSE2_PD(Y,F);
277             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
278             H                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) +2);
279             GMX_MM_TRANSPOSE2_PD(G,H);
280             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(vfeps,H,G),F);
281             VV               = _mm_macc_pd(vfeps,Fp,Y);
282             velec            = _mm_mul_pd(qq20,VV);
283             FF               = _mm_macc_pd(_mm_macc_pd(twovfeps,H,G),vfeps,Fp);
284             felec            = _mm_xor_pd(signbit,_mm_mul_pd(_mm_mul_pd(qq20,FF),_mm_mul_pd(vftabscale,rinv20)));
285
286             /* Update potential sum for this i atom from the interaction with this j atom. */
287             velecsum         = _mm_add_pd(velecsum,velec);
288
289             fscal            = felec;
290
291             /* Update vectorial force */
292             fix2             = _mm_macc_pd(dx20,fscal,fix2);
293             fiy2             = _mm_macc_pd(dy20,fscal,fiy2);
294             fiz2             = _mm_macc_pd(dz20,fscal,fiz2);
295             
296             fjx0             = _mm_macc_pd(dx20,fscal,fjx0);
297             fjy0             = _mm_macc_pd(dy20,fscal,fjy0);
298             fjz0             = _mm_macc_pd(dz20,fscal,fjz0);
299
300             /**************************
301              * CALCULATE INTERACTIONS *
302              **************************/
303
304             r30              = _mm_mul_pd(rsq30,rinv30);
305
306             /* Compute parameters for interactions between i and j atoms */
307             qq30             = _mm_mul_pd(iq3,jq0);
308
309             /* Calculate table index by multiplying r with table scale and truncate to integer */
310             rt               = _mm_mul_pd(r30,vftabscale);
311             vfitab           = _mm_cvttpd_epi32(rt);
312 #ifdef __XOP__
313             vfeps            = _mm_frcz_pd(rt);
314 #else
315             vfeps            = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
316 #endif
317             twovfeps         = _mm_add_pd(vfeps,vfeps);
318             vfitab           = _mm_slli_epi32(vfitab,2);
319
320             /* CUBIC SPLINE TABLE ELECTROSTATICS */
321             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
322             F                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) );
323             GMX_MM_TRANSPOSE2_PD(Y,F);
324             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
325             H                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) +2);
326             GMX_MM_TRANSPOSE2_PD(G,H);
327             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(vfeps,H,G),F);
328             VV               = _mm_macc_pd(vfeps,Fp,Y);
329             velec            = _mm_mul_pd(qq30,VV);
330             FF               = _mm_macc_pd(_mm_macc_pd(twovfeps,H,G),vfeps,Fp);
331             felec            = _mm_xor_pd(signbit,_mm_mul_pd(_mm_mul_pd(qq30,FF),_mm_mul_pd(vftabscale,rinv30)));
332
333             /* Update potential sum for this i atom from the interaction with this j atom. */
334             velecsum         = _mm_add_pd(velecsum,velec);
335
336             fscal            = felec;
337
338             /* Update vectorial force */
339             fix3             = _mm_macc_pd(dx30,fscal,fix3);
340             fiy3             = _mm_macc_pd(dy30,fscal,fiy3);
341             fiz3             = _mm_macc_pd(dz30,fscal,fiz3);
342             
343             fjx0             = _mm_macc_pd(dx30,fscal,fjx0);
344             fjy0             = _mm_macc_pd(dy30,fscal,fjy0);
345             fjz0             = _mm_macc_pd(dz30,fscal,fjz0);
346
347             gmx_mm_decrement_1rvec_2ptr_swizzle_pd(f+j_coord_offsetA,f+j_coord_offsetB,fjx0,fjy0,fjz0);
348
349             /* Inner loop uses 141 flops */
350         }
351
352         if(jidx<j_index_end)
353         {
354
355             jnrA             = jjnr[jidx];
356             j_coord_offsetA  = DIM*jnrA;
357
358             /* load j atom coordinates */
359             gmx_mm_load_1rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
360                                               &jx0,&jy0,&jz0);
361
362             /* Calculate displacement vector */
363             dx10             = _mm_sub_pd(ix1,jx0);
364             dy10             = _mm_sub_pd(iy1,jy0);
365             dz10             = _mm_sub_pd(iz1,jz0);
366             dx20             = _mm_sub_pd(ix2,jx0);
367             dy20             = _mm_sub_pd(iy2,jy0);
368             dz20             = _mm_sub_pd(iz2,jz0);
369             dx30             = _mm_sub_pd(ix3,jx0);
370             dy30             = _mm_sub_pd(iy3,jy0);
371             dz30             = _mm_sub_pd(iz3,jz0);
372
373             /* Calculate squared distance and things based on it */
374             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
375             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
376             rsq30            = gmx_mm_calc_rsq_pd(dx30,dy30,dz30);
377
378             rinv10           = gmx_mm_invsqrt_pd(rsq10);
379             rinv20           = gmx_mm_invsqrt_pd(rsq20);
380             rinv30           = gmx_mm_invsqrt_pd(rsq30);
381
382             /* Load parameters for j particles */
383             jq0              = _mm_load_sd(charge+jnrA+0);
384
385             fjx0             = _mm_setzero_pd();
386             fjy0             = _mm_setzero_pd();
387             fjz0             = _mm_setzero_pd();
388
389             /**************************
390              * CALCULATE INTERACTIONS *
391              **************************/
392
393             r10              = _mm_mul_pd(rsq10,rinv10);
394
395             /* Compute parameters for interactions between i and j atoms */
396             qq10             = _mm_mul_pd(iq1,jq0);
397
398             /* Calculate table index by multiplying r with table scale and truncate to integer */
399             rt               = _mm_mul_pd(r10,vftabscale);
400             vfitab           = _mm_cvttpd_epi32(rt);
401 #ifdef __XOP__
402             vfeps            = _mm_frcz_pd(rt);
403 #else
404             vfeps            = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
405 #endif
406             twovfeps         = _mm_add_pd(vfeps,vfeps);
407             vfitab           = _mm_slli_epi32(vfitab,2);
408
409             /* CUBIC SPLINE TABLE ELECTROSTATICS */
410             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
411             F                = _mm_setzero_pd();
412             GMX_MM_TRANSPOSE2_PD(Y,F);
413             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
414             H                = _mm_setzero_pd();
415             GMX_MM_TRANSPOSE2_PD(G,H);
416             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(vfeps,H,G),F);
417             VV               = _mm_macc_pd(vfeps,Fp,Y);
418             velec            = _mm_mul_pd(qq10,VV);
419             FF               = _mm_macc_pd(_mm_macc_pd(twovfeps,H,G),vfeps,Fp);
420             felec            = _mm_xor_pd(signbit,_mm_mul_pd(_mm_mul_pd(qq10,FF),_mm_mul_pd(vftabscale,rinv10)));
421
422             /* Update potential sum for this i atom from the interaction with this j atom. */
423             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
424             velecsum         = _mm_add_pd(velecsum,velec);
425
426             fscal            = felec;
427
428             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
429
430             /* Update vectorial force */
431             fix1             = _mm_macc_pd(dx10,fscal,fix1);
432             fiy1             = _mm_macc_pd(dy10,fscal,fiy1);
433             fiz1             = _mm_macc_pd(dz10,fscal,fiz1);
434             
435             fjx0             = _mm_macc_pd(dx10,fscal,fjx0);
436             fjy0             = _mm_macc_pd(dy10,fscal,fjy0);
437             fjz0             = _mm_macc_pd(dz10,fscal,fjz0);
438
439             /**************************
440              * CALCULATE INTERACTIONS *
441              **************************/
442
443             r20              = _mm_mul_pd(rsq20,rinv20);
444
445             /* Compute parameters for interactions between i and j atoms */
446             qq20             = _mm_mul_pd(iq2,jq0);
447
448             /* Calculate table index by multiplying r with table scale and truncate to integer */
449             rt               = _mm_mul_pd(r20,vftabscale);
450             vfitab           = _mm_cvttpd_epi32(rt);
451 #ifdef __XOP__
452             vfeps            = _mm_frcz_pd(rt);
453 #else
454             vfeps            = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
455 #endif
456             twovfeps         = _mm_add_pd(vfeps,vfeps);
457             vfitab           = _mm_slli_epi32(vfitab,2);
458
459             /* CUBIC SPLINE TABLE ELECTROSTATICS */
460             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
461             F                = _mm_setzero_pd();
462             GMX_MM_TRANSPOSE2_PD(Y,F);
463             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
464             H                = _mm_setzero_pd();
465             GMX_MM_TRANSPOSE2_PD(G,H);
466             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(vfeps,H,G),F);
467             VV               = _mm_macc_pd(vfeps,Fp,Y);
468             velec            = _mm_mul_pd(qq20,VV);
469             FF               = _mm_macc_pd(_mm_macc_pd(twovfeps,H,G),vfeps,Fp);
470             felec            = _mm_xor_pd(signbit,_mm_mul_pd(_mm_mul_pd(qq20,FF),_mm_mul_pd(vftabscale,rinv20)));
471
472             /* Update potential sum for this i atom from the interaction with this j atom. */
473             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
474             velecsum         = _mm_add_pd(velecsum,velec);
475
476             fscal            = felec;
477
478             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
479
480             /* Update vectorial force */
481             fix2             = _mm_macc_pd(dx20,fscal,fix2);
482             fiy2             = _mm_macc_pd(dy20,fscal,fiy2);
483             fiz2             = _mm_macc_pd(dz20,fscal,fiz2);
484             
485             fjx0             = _mm_macc_pd(dx20,fscal,fjx0);
486             fjy0             = _mm_macc_pd(dy20,fscal,fjy0);
487             fjz0             = _mm_macc_pd(dz20,fscal,fjz0);
488
489             /**************************
490              * CALCULATE INTERACTIONS *
491              **************************/
492
493             r30              = _mm_mul_pd(rsq30,rinv30);
494
495             /* Compute parameters for interactions between i and j atoms */
496             qq30             = _mm_mul_pd(iq3,jq0);
497
498             /* Calculate table index by multiplying r with table scale and truncate to integer */
499             rt               = _mm_mul_pd(r30,vftabscale);
500             vfitab           = _mm_cvttpd_epi32(rt);
501 #ifdef __XOP__
502             vfeps            = _mm_frcz_pd(rt);
503 #else
504             vfeps            = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
505 #endif
506             twovfeps         = _mm_add_pd(vfeps,vfeps);
507             vfitab           = _mm_slli_epi32(vfitab,2);
508
509             /* CUBIC SPLINE TABLE ELECTROSTATICS */
510             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
511             F                = _mm_setzero_pd();
512             GMX_MM_TRANSPOSE2_PD(Y,F);
513             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
514             H                = _mm_setzero_pd();
515             GMX_MM_TRANSPOSE2_PD(G,H);
516             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(vfeps,H,G),F);
517             VV               = _mm_macc_pd(vfeps,Fp,Y);
518             velec            = _mm_mul_pd(qq30,VV);
519             FF               = _mm_macc_pd(_mm_macc_pd(twovfeps,H,G),vfeps,Fp);
520             felec            = _mm_xor_pd(signbit,_mm_mul_pd(_mm_mul_pd(qq30,FF),_mm_mul_pd(vftabscale,rinv30)));
521
522             /* Update potential sum for this i atom from the interaction with this j atom. */
523             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
524             velecsum         = _mm_add_pd(velecsum,velec);
525
526             fscal            = felec;
527
528             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
529
530             /* Update vectorial force */
531             fix3             = _mm_macc_pd(dx30,fscal,fix3);
532             fiy3             = _mm_macc_pd(dy30,fscal,fiy3);
533             fiz3             = _mm_macc_pd(dz30,fscal,fiz3);
534             
535             fjx0             = _mm_macc_pd(dx30,fscal,fjx0);
536             fjy0             = _mm_macc_pd(dy30,fscal,fjy0);
537             fjz0             = _mm_macc_pd(dz30,fscal,fjz0);
538
539             gmx_mm_decrement_1rvec_1ptr_swizzle_pd(f+j_coord_offsetA,fjx0,fjy0,fjz0);
540
541             /* Inner loop uses 141 flops */
542         }
543
544         /* End of innermost loop */
545
546         gmx_mm_update_iforce_3atom_swizzle_pd(fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
547                                               f+i_coord_offset+DIM,fshift+i_shift_offset);
548
549         ggid                        = gid[iidx];
550         /* Update potential energies */
551         gmx_mm_update_1pot_pd(velecsum,kernel_data->energygrp_elec+ggid);
552
553         /* Increment number of inner iterations */
554         inneriter                  += j_index_end - j_index_start;
555
556         /* Outer loop uses 19 flops */
557     }
558
559     /* Increment number of outer iterations */
560     outeriter        += nri;
561
562     /* Update outer/inner flops */
563
564     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_W4_VF,outeriter*19 + inneriter*141);
565 }
566 /*
567  * Gromacs nonbonded kernel:   nb_kernel_ElecCSTab_VdwNone_GeomW4P1_F_avx_128_fma_double
568  * Electrostatics interaction: CubicSplineTable
569  * VdW interaction:            None
570  * Geometry:                   Water4-Particle
571  * Calculate force/pot:        Force
572  */
573 void
574 nb_kernel_ElecCSTab_VdwNone_GeomW4P1_F_avx_128_fma_double
575                     (t_nblist                    * gmx_restrict       nlist,
576                      rvec                        * gmx_restrict          xx,
577                      rvec                        * gmx_restrict          ff,
578                      t_forcerec                  * gmx_restrict          fr,
579                      t_mdatoms                   * gmx_restrict     mdatoms,
580                      nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
581                      t_nrnb                      * gmx_restrict        nrnb)
582 {
583     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
584      * just 0 for non-waters.
585      * Suffixes A,B refer to j loop unrolling done with SSE double precision, e.g. for the two different
586      * jnr indices corresponding to data put in the four positions in the SIMD register.
587      */
588     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
589     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
590     int              jnrA,jnrB;
591     int              j_coord_offsetA,j_coord_offsetB;
592     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
593     real             rcutoff_scalar;
594     real             *shiftvec,*fshift,*x,*f;
595     __m128d          tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
596     int              vdwioffset1;
597     __m128d          ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
598     int              vdwioffset2;
599     __m128d          ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
600     int              vdwioffset3;
601     __m128d          ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
602     int              vdwjidx0A,vdwjidx0B;
603     __m128d          jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
604     __m128d          dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
605     __m128d          dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
606     __m128d          dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
607     __m128d          velec,felec,velecsum,facel,crf,krf,krf2;
608     real             *charge;
609     __m128i          vfitab;
610     __m128i          ifour       = _mm_set1_epi32(4);
611     __m128d          rt,vfeps,vftabscale,Y,F,G,H,Heps,Fp,VV,FF,twovfeps;
612     real             *vftab;
613     __m128d          dummy_mask,cutoff_mask;
614     __m128d          signbit   = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
615     __m128d          one     = _mm_set1_pd(1.0);
616     __m128d          two     = _mm_set1_pd(2.0);
617     x                = xx[0];
618     f                = ff[0];
619
620     nri              = nlist->nri;
621     iinr             = nlist->iinr;
622     jindex           = nlist->jindex;
623     jjnr             = nlist->jjnr;
624     shiftidx         = nlist->shift;
625     gid              = nlist->gid;
626     shiftvec         = fr->shift_vec[0];
627     fshift           = fr->fshift[0];
628     facel            = _mm_set1_pd(fr->epsfac);
629     charge           = mdatoms->chargeA;
630
631     vftab            = kernel_data->table_elec->data;
632     vftabscale       = _mm_set1_pd(kernel_data->table_elec->scale);
633
634     /* Setup water-specific parameters */
635     inr              = nlist->iinr[0];
636     iq1              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+1]));
637     iq2              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+2]));
638     iq3              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+3]));
639
640     /* Avoid stupid compiler warnings */
641     jnrA = jnrB = 0;
642     j_coord_offsetA = 0;
643     j_coord_offsetB = 0;
644
645     outeriter        = 0;
646     inneriter        = 0;
647
648     /* Start outer loop over neighborlists */
649     for(iidx=0; iidx<nri; iidx++)
650     {
651         /* Load shift vector for this list */
652         i_shift_offset   = DIM*shiftidx[iidx];
653
654         /* Load limits for loop over neighbors */
655         j_index_start    = jindex[iidx];
656         j_index_end      = jindex[iidx+1];
657
658         /* Get outer coordinate index */
659         inr              = iinr[iidx];
660         i_coord_offset   = DIM*inr;
661
662         /* Load i particle coords and add shift vector */
663         gmx_mm_load_shift_and_3rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset+DIM,
664                                                  &ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
665
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         fix3             = _mm_setzero_pd();
673         fiy3             = _mm_setzero_pd();
674         fiz3             = _mm_setzero_pd();
675
676         /* Start inner kernel loop */
677         for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
678         {
679
680             /* Get j neighbor index, and coordinate index */
681             jnrA             = jjnr[jidx];
682             jnrB             = jjnr[jidx+1];
683             j_coord_offsetA  = DIM*jnrA;
684             j_coord_offsetB  = DIM*jnrB;
685
686             /* load j atom coordinates */
687             gmx_mm_load_1rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
688                                               &jx0,&jy0,&jz0);
689
690             /* Calculate displacement vector */
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             dx30             = _mm_sub_pd(ix3,jx0);
698             dy30             = _mm_sub_pd(iy3,jy0);
699             dz30             = _mm_sub_pd(iz3,jz0);
700
701             /* Calculate squared distance and things based on it */
702             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
703             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
704             rsq30            = gmx_mm_calc_rsq_pd(dx30,dy30,dz30);
705
706             rinv10           = gmx_mm_invsqrt_pd(rsq10);
707             rinv20           = gmx_mm_invsqrt_pd(rsq20);
708             rinv30           = gmx_mm_invsqrt_pd(rsq30);
709
710             /* Load parameters for j particles */
711             jq0              = gmx_mm_load_2real_swizzle_pd(charge+jnrA+0,charge+jnrB+0);
712
713             fjx0             = _mm_setzero_pd();
714             fjy0             = _mm_setzero_pd();
715             fjz0             = _mm_setzero_pd();
716
717             /**************************
718              * CALCULATE INTERACTIONS *
719              **************************/
720
721             r10              = _mm_mul_pd(rsq10,rinv10);
722
723             /* Compute parameters for interactions between i and j atoms */
724             qq10             = _mm_mul_pd(iq1,jq0);
725
726             /* Calculate table index by multiplying r with table scale and truncate to integer */
727             rt               = _mm_mul_pd(r10,vftabscale);
728             vfitab           = _mm_cvttpd_epi32(rt);
729 #ifdef __XOP__
730             vfeps            = _mm_frcz_pd(rt);
731 #else
732             vfeps            = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
733 #endif
734             twovfeps         = _mm_add_pd(vfeps,vfeps);
735             vfitab           = _mm_slli_epi32(vfitab,2);
736
737             /* CUBIC SPLINE TABLE ELECTROSTATICS */
738             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
739             F                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) );
740             GMX_MM_TRANSPOSE2_PD(Y,F);
741             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
742             H                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) +2);
743             GMX_MM_TRANSPOSE2_PD(G,H);
744             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(vfeps,H,G),F);
745             FF               = _mm_macc_pd(_mm_macc_pd(twovfeps,H,G),vfeps,Fp);
746             felec            = _mm_xor_pd(signbit,_mm_mul_pd(_mm_mul_pd(qq10,FF),_mm_mul_pd(vftabscale,rinv10)));
747
748             fscal            = felec;
749
750             /* Update vectorial force */
751             fix1             = _mm_macc_pd(dx10,fscal,fix1);
752             fiy1             = _mm_macc_pd(dy10,fscal,fiy1);
753             fiz1             = _mm_macc_pd(dz10,fscal,fiz1);
754             
755             fjx0             = _mm_macc_pd(dx10,fscal,fjx0);
756             fjy0             = _mm_macc_pd(dy10,fscal,fjy0);
757             fjz0             = _mm_macc_pd(dz10,fscal,fjz0);
758
759             /**************************
760              * CALCULATE INTERACTIONS *
761              **************************/
762
763             r20              = _mm_mul_pd(rsq20,rinv20);
764
765             /* Compute parameters for interactions between i and j atoms */
766             qq20             = _mm_mul_pd(iq2,jq0);
767
768             /* Calculate table index by multiplying r with table scale and truncate to integer */
769             rt               = _mm_mul_pd(r20,vftabscale);
770             vfitab           = _mm_cvttpd_epi32(rt);
771 #ifdef __XOP__
772             vfeps            = _mm_frcz_pd(rt);
773 #else
774             vfeps            = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
775 #endif
776             twovfeps         = _mm_add_pd(vfeps,vfeps);
777             vfitab           = _mm_slli_epi32(vfitab,2);
778
779             /* CUBIC SPLINE TABLE ELECTROSTATICS */
780             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
781             F                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) );
782             GMX_MM_TRANSPOSE2_PD(Y,F);
783             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
784             H                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) +2);
785             GMX_MM_TRANSPOSE2_PD(G,H);
786             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(vfeps,H,G),F);
787             FF               = _mm_macc_pd(_mm_macc_pd(twovfeps,H,G),vfeps,Fp);
788             felec            = _mm_xor_pd(signbit,_mm_mul_pd(_mm_mul_pd(qq20,FF),_mm_mul_pd(vftabscale,rinv20)));
789
790             fscal            = felec;
791
792             /* Update vectorial force */
793             fix2             = _mm_macc_pd(dx20,fscal,fix2);
794             fiy2             = _mm_macc_pd(dy20,fscal,fiy2);
795             fiz2             = _mm_macc_pd(dz20,fscal,fiz2);
796             
797             fjx0             = _mm_macc_pd(dx20,fscal,fjx0);
798             fjy0             = _mm_macc_pd(dy20,fscal,fjy0);
799             fjz0             = _mm_macc_pd(dz20,fscal,fjz0);
800
801             /**************************
802              * CALCULATE INTERACTIONS *
803              **************************/
804
805             r30              = _mm_mul_pd(rsq30,rinv30);
806
807             /* Compute parameters for interactions between i and j atoms */
808             qq30             = _mm_mul_pd(iq3,jq0);
809
810             /* Calculate table index by multiplying r with table scale and truncate to integer */
811             rt               = _mm_mul_pd(r30,vftabscale);
812             vfitab           = _mm_cvttpd_epi32(rt);
813 #ifdef __XOP__
814             vfeps            = _mm_frcz_pd(rt);
815 #else
816             vfeps            = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
817 #endif
818             twovfeps         = _mm_add_pd(vfeps,vfeps);
819             vfitab           = _mm_slli_epi32(vfitab,2);
820
821             /* CUBIC SPLINE TABLE ELECTROSTATICS */
822             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
823             F                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) );
824             GMX_MM_TRANSPOSE2_PD(Y,F);
825             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
826             H                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,1) +2);
827             GMX_MM_TRANSPOSE2_PD(G,H);
828             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(vfeps,H,G),F);
829             FF               = _mm_macc_pd(_mm_macc_pd(twovfeps,H,G),vfeps,Fp);
830             felec            = _mm_xor_pd(signbit,_mm_mul_pd(_mm_mul_pd(qq30,FF),_mm_mul_pd(vftabscale,rinv30)));
831
832             fscal            = felec;
833
834             /* Update vectorial force */
835             fix3             = _mm_macc_pd(dx30,fscal,fix3);
836             fiy3             = _mm_macc_pd(dy30,fscal,fiy3);
837             fiz3             = _mm_macc_pd(dz30,fscal,fiz3);
838             
839             fjx0             = _mm_macc_pd(dx30,fscal,fjx0);
840             fjy0             = _mm_macc_pd(dy30,fscal,fjy0);
841             fjz0             = _mm_macc_pd(dz30,fscal,fjz0);
842
843             gmx_mm_decrement_1rvec_2ptr_swizzle_pd(f+j_coord_offsetA,f+j_coord_offsetB,fjx0,fjy0,fjz0);
844
845             /* Inner loop uses 129 flops */
846         }
847
848         if(jidx<j_index_end)
849         {
850
851             jnrA             = jjnr[jidx];
852             j_coord_offsetA  = DIM*jnrA;
853
854             /* load j atom coordinates */
855             gmx_mm_load_1rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
856                                               &jx0,&jy0,&jz0);
857
858             /* Calculate displacement vector */
859             dx10             = _mm_sub_pd(ix1,jx0);
860             dy10             = _mm_sub_pd(iy1,jy0);
861             dz10             = _mm_sub_pd(iz1,jz0);
862             dx20             = _mm_sub_pd(ix2,jx0);
863             dy20             = _mm_sub_pd(iy2,jy0);
864             dz20             = _mm_sub_pd(iz2,jz0);
865             dx30             = _mm_sub_pd(ix3,jx0);
866             dy30             = _mm_sub_pd(iy3,jy0);
867             dz30             = _mm_sub_pd(iz3,jz0);
868
869             /* Calculate squared distance and things based on it */
870             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
871             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
872             rsq30            = gmx_mm_calc_rsq_pd(dx30,dy30,dz30);
873
874             rinv10           = gmx_mm_invsqrt_pd(rsq10);
875             rinv20           = gmx_mm_invsqrt_pd(rsq20);
876             rinv30           = gmx_mm_invsqrt_pd(rsq30);
877
878             /* Load parameters for j particles */
879             jq0              = _mm_load_sd(charge+jnrA+0);
880
881             fjx0             = _mm_setzero_pd();
882             fjy0             = _mm_setzero_pd();
883             fjz0             = _mm_setzero_pd();
884
885             /**************************
886              * CALCULATE INTERACTIONS *
887              **************************/
888
889             r10              = _mm_mul_pd(rsq10,rinv10);
890
891             /* Compute parameters for interactions between i and j atoms */
892             qq10             = _mm_mul_pd(iq1,jq0);
893
894             /* Calculate table index by multiplying r with table scale and truncate to integer */
895             rt               = _mm_mul_pd(r10,vftabscale);
896             vfitab           = _mm_cvttpd_epi32(rt);
897 #ifdef __XOP__
898             vfeps            = _mm_frcz_pd(rt);
899 #else
900             vfeps            = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
901 #endif
902             twovfeps         = _mm_add_pd(vfeps,vfeps);
903             vfitab           = _mm_slli_epi32(vfitab,2);
904
905             /* CUBIC SPLINE TABLE ELECTROSTATICS */
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(vfeps,H,G),F);
913             FF               = _mm_macc_pd(_mm_macc_pd(twovfeps,H,G),vfeps,Fp);
914             felec            = _mm_xor_pd(signbit,_mm_mul_pd(_mm_mul_pd(qq10,FF),_mm_mul_pd(vftabscale,rinv10)));
915
916             fscal            = felec;
917
918             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
919
920             /* Update vectorial force */
921             fix1             = _mm_macc_pd(dx10,fscal,fix1);
922             fiy1             = _mm_macc_pd(dy10,fscal,fiy1);
923             fiz1             = _mm_macc_pd(dz10,fscal,fiz1);
924             
925             fjx0             = _mm_macc_pd(dx10,fscal,fjx0);
926             fjy0             = _mm_macc_pd(dy10,fscal,fjy0);
927             fjz0             = _mm_macc_pd(dz10,fscal,fjz0);
928
929             /**************************
930              * CALCULATE INTERACTIONS *
931              **************************/
932
933             r20              = _mm_mul_pd(rsq20,rinv20);
934
935             /* Compute parameters for interactions between i and j atoms */
936             qq20             = _mm_mul_pd(iq2,jq0);
937
938             /* Calculate table index by multiplying r with table scale and truncate to integer */
939             rt               = _mm_mul_pd(r20,vftabscale);
940             vfitab           = _mm_cvttpd_epi32(rt);
941 #ifdef __XOP__
942             vfeps            = _mm_frcz_pd(rt);
943 #else
944             vfeps            = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
945 #endif
946             twovfeps         = _mm_add_pd(vfeps,vfeps);
947             vfitab           = _mm_slli_epi32(vfitab,2);
948
949             /* CUBIC SPLINE TABLE ELECTROSTATICS */
950             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
951             F                = _mm_setzero_pd();
952             GMX_MM_TRANSPOSE2_PD(Y,F);
953             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
954             H                = _mm_setzero_pd();
955             GMX_MM_TRANSPOSE2_PD(G,H);
956             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(vfeps,H,G),F);
957             FF               = _mm_macc_pd(_mm_macc_pd(twovfeps,H,G),vfeps,Fp);
958             felec            = _mm_xor_pd(signbit,_mm_mul_pd(_mm_mul_pd(qq20,FF),_mm_mul_pd(vftabscale,rinv20)));
959
960             fscal            = felec;
961
962             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
963
964             /* Update vectorial force */
965             fix2             = _mm_macc_pd(dx20,fscal,fix2);
966             fiy2             = _mm_macc_pd(dy20,fscal,fiy2);
967             fiz2             = _mm_macc_pd(dz20,fscal,fiz2);
968             
969             fjx0             = _mm_macc_pd(dx20,fscal,fjx0);
970             fjy0             = _mm_macc_pd(dy20,fscal,fjy0);
971             fjz0             = _mm_macc_pd(dz20,fscal,fjz0);
972
973             /**************************
974              * CALCULATE INTERACTIONS *
975              **************************/
976
977             r30              = _mm_mul_pd(rsq30,rinv30);
978
979             /* Compute parameters for interactions between i and j atoms */
980             qq30             = _mm_mul_pd(iq3,jq0);
981
982             /* Calculate table index by multiplying r with table scale and truncate to integer */
983             rt               = _mm_mul_pd(r30,vftabscale);
984             vfitab           = _mm_cvttpd_epi32(rt);
985 #ifdef __XOP__
986             vfeps            = _mm_frcz_pd(rt);
987 #else
988             vfeps            = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
989 #endif
990             twovfeps         = _mm_add_pd(vfeps,vfeps);
991             vfitab           = _mm_slli_epi32(vfitab,2);
992
993             /* CUBIC SPLINE TABLE ELECTROSTATICS */
994             Y                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) );
995             F                = _mm_setzero_pd();
996             GMX_MM_TRANSPOSE2_PD(Y,F);
997             G                = _mm_load_pd( vftab + _mm_extract_epi32(vfitab,0) +2);
998             H                = _mm_setzero_pd();
999             GMX_MM_TRANSPOSE2_PD(G,H);
1000             Fp               = _mm_macc_pd(vfeps,_mm_macc_pd(vfeps,H,G),F);
1001             FF               = _mm_macc_pd(_mm_macc_pd(twovfeps,H,G),vfeps,Fp);
1002             felec            = _mm_xor_pd(signbit,_mm_mul_pd(_mm_mul_pd(qq30,FF),_mm_mul_pd(vftabscale,rinv30)));
1003
1004             fscal            = felec;
1005
1006             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1007
1008             /* Update vectorial force */
1009             fix3             = _mm_macc_pd(dx30,fscal,fix3);
1010             fiy3             = _mm_macc_pd(dy30,fscal,fiy3);
1011             fiz3             = _mm_macc_pd(dz30,fscal,fiz3);
1012             
1013             fjx0             = _mm_macc_pd(dx30,fscal,fjx0);
1014             fjy0             = _mm_macc_pd(dy30,fscal,fjy0);
1015             fjz0             = _mm_macc_pd(dz30,fscal,fjz0);
1016
1017             gmx_mm_decrement_1rvec_1ptr_swizzle_pd(f+j_coord_offsetA,fjx0,fjy0,fjz0);
1018
1019             /* Inner loop uses 129 flops */
1020         }
1021
1022         /* End of innermost loop */
1023
1024         gmx_mm_update_iforce_3atom_swizzle_pd(fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
1025                                               f+i_coord_offset+DIM,fshift+i_shift_offset);
1026
1027         /* Increment number of inner iterations */
1028         inneriter                  += j_index_end - j_index_start;
1029
1030         /* Outer loop uses 18 flops */
1031     }
1032
1033     /* Increment number of outer iterations */
1034     outeriter        += nri;
1035
1036     /* Update outer/inner flops */
1037
1038     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_W4_F,outeriter*18 + inneriter*129);
1039 }