Merge release-4-6 into master
[alexxy/gromacs.git] / src / gromacs / gmxlib / nonbonded / nb_kernel_avx_128_fma_single / nb_kernel_ElecCoul_VdwNone_GeomW4P1_avx_128_fma_single.c
1 /*
2  * Note: this file was generated by the Gromacs avx_128_fma_single kernel generator.
3  *
4  *                This source code is part of
5  *
6  *                 G   R   O   M   A   C   S
7  *
8  * Copyright (c) 2001-2012, The GROMACS Development Team
9  *
10  * Gromacs is a library for molecular simulation and trajectory analysis,
11  * written by Erik Lindahl, David van der Spoel, Berk Hess, and others - for
12  * a full list of developers and information, check out http://www.gromacs.org
13  *
14  * This program is free software; you can redistribute it and/or modify it under
15  * the terms of the GNU Lesser General Public License as published by the Free
16  * Software Foundation; either version 2 of the License, or (at your option) any
17  * later version.
18  *
19  * To help fund GROMACS development, we humbly ask that you cite
20  * the papers people have written on it - you can find them on the website.
21  */
22 #ifdef HAVE_CONFIG_H
23 #include <config.h>
24 #endif
25
26 #include <math.h>
27
28 #include "../nb_kernel.h"
29 #include "types/simple.h"
30 #include "vec.h"
31 #include "nrnb.h"
32
33 #include "gmx_math_x86_avx_128_fma_single.h"
34 #include "kernelutil_x86_avx_128_fma_single.h"
35
36 /*
37  * Gromacs nonbonded kernel:   nb_kernel_ElecCoul_VdwNone_GeomW4P1_VF_avx_128_fma_single
38  * Electrostatics interaction: Coulomb
39  * VdW interaction:            None
40  * Geometry:                   Water4-Particle
41  * Calculate force/pot:        PotentialAndForce
42  */
43 void
44 nb_kernel_ElecCoul_VdwNone_GeomW4P1_VF_avx_128_fma_single
45                     (t_nblist * gmx_restrict                nlist,
46                      rvec * gmx_restrict                    xx,
47                      rvec * gmx_restrict                    ff,
48                      t_forcerec * gmx_restrict              fr,
49                      t_mdatoms * gmx_restrict               mdatoms,
50                      nb_kernel_data_t * gmx_restrict        kernel_data,
51                      t_nrnb * gmx_restrict                  nrnb)
52 {
53     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
54      * just 0 for non-waters.
55      * Suffixes A,B,C,D refer to j loop unrolling done with AVX_128, e.g. for the four different
56      * jnr indices corresponding to data put in the four positions in the SIMD register.
57      */
58     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
59     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
60     int              jnrA,jnrB,jnrC,jnrD;
61     int              jnrlistA,jnrlistB,jnrlistC,jnrlistD;
62     int              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
63     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
64     real             rcutoff_scalar;
65     real             *shiftvec,*fshift,*x,*f;
66     real             *fjptrA,*fjptrB,*fjptrC,*fjptrD;
67     real             scratch[4*DIM];
68     __m128           fscal,rcutoff,rcutoff2,jidxall;
69     int              vdwioffset1;
70     __m128           ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
71     int              vdwioffset2;
72     __m128           ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
73     int              vdwioffset3;
74     __m128           ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
75     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
76     __m128           jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
77     __m128           dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
78     __m128           dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
79     __m128           dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
80     __m128           velec,felec,velecsum,facel,crf,krf,krf2;
81     real             *charge;
82     __m128           dummy_mask,cutoff_mask;
83     __m128           signbit = _mm_castsi128_ps( _mm_set1_epi32(0x80000000) );
84     __m128           one     = _mm_set1_ps(1.0);
85     __m128           two     = _mm_set1_ps(2.0);
86     x                = xx[0];
87     f                = ff[0];
88
89     nri              = nlist->nri;
90     iinr             = nlist->iinr;
91     jindex           = nlist->jindex;
92     jjnr             = nlist->jjnr;
93     shiftidx         = nlist->shift;
94     gid              = nlist->gid;
95     shiftvec         = fr->shift_vec[0];
96     fshift           = fr->fshift[0];
97     facel            = _mm_set1_ps(fr->epsfac);
98     charge           = mdatoms->chargeA;
99
100     /* Setup water-specific parameters */
101     inr              = nlist->iinr[0];
102     iq1              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+1]));
103     iq2              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+2]));
104     iq3              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+3]));
105
106     /* Avoid stupid compiler warnings */
107     jnrA = jnrB = jnrC = jnrD = 0;
108     j_coord_offsetA = 0;
109     j_coord_offsetB = 0;
110     j_coord_offsetC = 0;
111     j_coord_offsetD = 0;
112
113     outeriter        = 0;
114     inneriter        = 0;
115
116     for(iidx=0;iidx<4*DIM;iidx++)
117     {
118         scratch[iidx] = 0.0;
119     }
120
121     /* Start outer loop over neighborlists */
122     for(iidx=0; iidx<nri; iidx++)
123     {
124         /* Load shift vector for this list */
125         i_shift_offset   = DIM*shiftidx[iidx];
126
127         /* Load limits for loop over neighbors */
128         j_index_start    = jindex[iidx];
129         j_index_end      = jindex[iidx+1];
130
131         /* Get outer coordinate index */
132         inr              = iinr[iidx];
133         i_coord_offset   = DIM*inr;
134
135         /* Load i particle coords and add shift vector */
136         gmx_mm_load_shift_and_3rvec_broadcast_ps(shiftvec+i_shift_offset,x+i_coord_offset+DIM,
137                                                  &ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
138
139         fix1             = _mm_setzero_ps();
140         fiy1             = _mm_setzero_ps();
141         fiz1             = _mm_setzero_ps();
142         fix2             = _mm_setzero_ps();
143         fiy2             = _mm_setzero_ps();
144         fiz2             = _mm_setzero_ps();
145         fix3             = _mm_setzero_ps();
146         fiy3             = _mm_setzero_ps();
147         fiz3             = _mm_setzero_ps();
148
149         /* Reset potential sums */
150         velecsum         = _mm_setzero_ps();
151
152         /* Start inner kernel loop */
153         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
154         {
155
156             /* Get j neighbor index, and coordinate index */
157             jnrA             = jjnr[jidx];
158             jnrB             = jjnr[jidx+1];
159             jnrC             = jjnr[jidx+2];
160             jnrD             = jjnr[jidx+3];
161             j_coord_offsetA  = DIM*jnrA;
162             j_coord_offsetB  = DIM*jnrB;
163             j_coord_offsetC  = DIM*jnrC;
164             j_coord_offsetD  = DIM*jnrD;
165
166             /* load j atom coordinates */
167             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
168                                               x+j_coord_offsetC,x+j_coord_offsetD,
169                                               &jx0,&jy0,&jz0);
170
171             /* Calculate displacement vector */
172             dx10             = _mm_sub_ps(ix1,jx0);
173             dy10             = _mm_sub_ps(iy1,jy0);
174             dz10             = _mm_sub_ps(iz1,jz0);
175             dx20             = _mm_sub_ps(ix2,jx0);
176             dy20             = _mm_sub_ps(iy2,jy0);
177             dz20             = _mm_sub_ps(iz2,jz0);
178             dx30             = _mm_sub_ps(ix3,jx0);
179             dy30             = _mm_sub_ps(iy3,jy0);
180             dz30             = _mm_sub_ps(iz3,jz0);
181
182             /* Calculate squared distance and things based on it */
183             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
184             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
185             rsq30            = gmx_mm_calc_rsq_ps(dx30,dy30,dz30);
186
187             rinv10           = gmx_mm_invsqrt_ps(rsq10);
188             rinv20           = gmx_mm_invsqrt_ps(rsq20);
189             rinv30           = gmx_mm_invsqrt_ps(rsq30);
190
191             rinvsq10         = _mm_mul_ps(rinv10,rinv10);
192             rinvsq20         = _mm_mul_ps(rinv20,rinv20);
193             rinvsq30         = _mm_mul_ps(rinv30,rinv30);
194
195             /* Load parameters for j particles */
196             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
197                                                               charge+jnrC+0,charge+jnrD+0);
198
199             fjx0             = _mm_setzero_ps();
200             fjy0             = _mm_setzero_ps();
201             fjz0             = _mm_setzero_ps();
202
203             /**************************
204              * CALCULATE INTERACTIONS *
205              **************************/
206
207             /* Compute parameters for interactions between i and j atoms */
208             qq10             = _mm_mul_ps(iq1,jq0);
209
210             /* COULOMB ELECTROSTATICS */
211             velec            = _mm_mul_ps(qq10,rinv10);
212             felec            = _mm_mul_ps(velec,rinvsq10);
213
214             /* Update potential sum for this i atom from the interaction with this j atom. */
215             velecsum         = _mm_add_ps(velecsum,velec);
216
217             fscal            = felec;
218
219              /* Update vectorial force */
220             fix1             = _mm_macc_ps(dx10,fscal,fix1);
221             fiy1             = _mm_macc_ps(dy10,fscal,fiy1);
222             fiz1             = _mm_macc_ps(dz10,fscal,fiz1);
223
224             fjx0             = _mm_macc_ps(dx10,fscal,fjx0);
225             fjy0             = _mm_macc_ps(dy10,fscal,fjy0);
226             fjz0             = _mm_macc_ps(dz10,fscal,fjz0);
227
228             /**************************
229              * CALCULATE INTERACTIONS *
230              **************************/
231
232             /* Compute parameters for interactions between i and j atoms */
233             qq20             = _mm_mul_ps(iq2,jq0);
234
235             /* COULOMB ELECTROSTATICS */
236             velec            = _mm_mul_ps(qq20,rinv20);
237             felec            = _mm_mul_ps(velec,rinvsq20);
238
239             /* Update potential sum for this i atom from the interaction with this j atom. */
240             velecsum         = _mm_add_ps(velecsum,velec);
241
242             fscal            = felec;
243
244              /* Update vectorial force */
245             fix2             = _mm_macc_ps(dx20,fscal,fix2);
246             fiy2             = _mm_macc_ps(dy20,fscal,fiy2);
247             fiz2             = _mm_macc_ps(dz20,fscal,fiz2);
248
249             fjx0             = _mm_macc_ps(dx20,fscal,fjx0);
250             fjy0             = _mm_macc_ps(dy20,fscal,fjy0);
251             fjz0             = _mm_macc_ps(dz20,fscal,fjz0);
252
253             /**************************
254              * CALCULATE INTERACTIONS *
255              **************************/
256
257             /* Compute parameters for interactions between i and j atoms */
258             qq30             = _mm_mul_ps(iq3,jq0);
259
260             /* COULOMB ELECTROSTATICS */
261             velec            = _mm_mul_ps(qq30,rinv30);
262             felec            = _mm_mul_ps(velec,rinvsq30);
263
264             /* Update potential sum for this i atom from the interaction with this j atom. */
265             velecsum         = _mm_add_ps(velecsum,velec);
266
267             fscal            = felec;
268
269              /* Update vectorial force */
270             fix3             = _mm_macc_ps(dx30,fscal,fix3);
271             fiy3             = _mm_macc_ps(dy30,fscal,fiy3);
272             fiz3             = _mm_macc_ps(dz30,fscal,fiz3);
273
274             fjx0             = _mm_macc_ps(dx30,fscal,fjx0);
275             fjy0             = _mm_macc_ps(dy30,fscal,fjy0);
276             fjz0             = _mm_macc_ps(dz30,fscal,fjz0);
277
278             fjptrA             = f+j_coord_offsetA;
279             fjptrB             = f+j_coord_offsetB;
280             fjptrC             = f+j_coord_offsetC;
281             fjptrD             = f+j_coord_offsetD;
282
283             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,fjx0,fjy0,fjz0);
284
285             /* Inner loop uses 93 flops */
286         }
287
288         if(jidx<j_index_end)
289         {
290
291             /* Get j neighbor index, and coordinate index */
292             jnrlistA         = jjnr[jidx];
293             jnrlistB         = jjnr[jidx+1];
294             jnrlistC         = jjnr[jidx+2];
295             jnrlistD         = jjnr[jidx+3];
296             /* Sign of each element will be negative for non-real atoms.
297              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
298              * so use it as val = _mm_andnot_ps(mask,val) to clear dummy entries.
299              */
300             dummy_mask = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
301             jnrA       = (jnrlistA>=0) ? jnrlistA : 0;
302             jnrB       = (jnrlistB>=0) ? jnrlistB : 0;
303             jnrC       = (jnrlistC>=0) ? jnrlistC : 0;
304             jnrD       = (jnrlistD>=0) ? jnrlistD : 0;
305             j_coord_offsetA  = DIM*jnrA;
306             j_coord_offsetB  = DIM*jnrB;
307             j_coord_offsetC  = DIM*jnrC;
308             j_coord_offsetD  = DIM*jnrD;
309
310             /* load j atom coordinates */
311             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
312                                               x+j_coord_offsetC,x+j_coord_offsetD,
313                                               &jx0,&jy0,&jz0);
314
315             /* Calculate displacement vector */
316             dx10             = _mm_sub_ps(ix1,jx0);
317             dy10             = _mm_sub_ps(iy1,jy0);
318             dz10             = _mm_sub_ps(iz1,jz0);
319             dx20             = _mm_sub_ps(ix2,jx0);
320             dy20             = _mm_sub_ps(iy2,jy0);
321             dz20             = _mm_sub_ps(iz2,jz0);
322             dx30             = _mm_sub_ps(ix3,jx0);
323             dy30             = _mm_sub_ps(iy3,jy0);
324             dz30             = _mm_sub_ps(iz3,jz0);
325
326             /* Calculate squared distance and things based on it */
327             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
328             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
329             rsq30            = gmx_mm_calc_rsq_ps(dx30,dy30,dz30);
330
331             rinv10           = gmx_mm_invsqrt_ps(rsq10);
332             rinv20           = gmx_mm_invsqrt_ps(rsq20);
333             rinv30           = gmx_mm_invsqrt_ps(rsq30);
334
335             rinvsq10         = _mm_mul_ps(rinv10,rinv10);
336             rinvsq20         = _mm_mul_ps(rinv20,rinv20);
337             rinvsq30         = _mm_mul_ps(rinv30,rinv30);
338
339             /* Load parameters for j particles */
340             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
341                                                               charge+jnrC+0,charge+jnrD+0);
342
343             fjx0             = _mm_setzero_ps();
344             fjy0             = _mm_setzero_ps();
345             fjz0             = _mm_setzero_ps();
346
347             /**************************
348              * CALCULATE INTERACTIONS *
349              **************************/
350
351             /* Compute parameters for interactions between i and j atoms */
352             qq10             = _mm_mul_ps(iq1,jq0);
353
354             /* COULOMB ELECTROSTATICS */
355             velec            = _mm_mul_ps(qq10,rinv10);
356             felec            = _mm_mul_ps(velec,rinvsq10);
357
358             /* Update potential sum for this i atom from the interaction with this j atom. */
359             velec            = _mm_andnot_ps(dummy_mask,velec);
360             velecsum         = _mm_add_ps(velecsum,velec);
361
362             fscal            = felec;
363
364             fscal            = _mm_andnot_ps(dummy_mask,fscal);
365
366              /* Update vectorial force */
367             fix1             = _mm_macc_ps(dx10,fscal,fix1);
368             fiy1             = _mm_macc_ps(dy10,fscal,fiy1);
369             fiz1             = _mm_macc_ps(dz10,fscal,fiz1);
370
371             fjx0             = _mm_macc_ps(dx10,fscal,fjx0);
372             fjy0             = _mm_macc_ps(dy10,fscal,fjy0);
373             fjz0             = _mm_macc_ps(dz10,fscal,fjz0);
374
375             /**************************
376              * CALCULATE INTERACTIONS *
377              **************************/
378
379             /* Compute parameters for interactions between i and j atoms */
380             qq20             = _mm_mul_ps(iq2,jq0);
381
382             /* COULOMB ELECTROSTATICS */
383             velec            = _mm_mul_ps(qq20,rinv20);
384             felec            = _mm_mul_ps(velec,rinvsq20);
385
386             /* Update potential sum for this i atom from the interaction with this j atom. */
387             velec            = _mm_andnot_ps(dummy_mask,velec);
388             velecsum         = _mm_add_ps(velecsum,velec);
389
390             fscal            = felec;
391
392             fscal            = _mm_andnot_ps(dummy_mask,fscal);
393
394              /* Update vectorial force */
395             fix2             = _mm_macc_ps(dx20,fscal,fix2);
396             fiy2             = _mm_macc_ps(dy20,fscal,fiy2);
397             fiz2             = _mm_macc_ps(dz20,fscal,fiz2);
398
399             fjx0             = _mm_macc_ps(dx20,fscal,fjx0);
400             fjy0             = _mm_macc_ps(dy20,fscal,fjy0);
401             fjz0             = _mm_macc_ps(dz20,fscal,fjz0);
402
403             /**************************
404              * CALCULATE INTERACTIONS *
405              **************************/
406
407             /* Compute parameters for interactions between i and j atoms */
408             qq30             = _mm_mul_ps(iq3,jq0);
409
410             /* COULOMB ELECTROSTATICS */
411             velec            = _mm_mul_ps(qq30,rinv30);
412             felec            = _mm_mul_ps(velec,rinvsq30);
413
414             /* Update potential sum for this i atom from the interaction with this j atom. */
415             velec            = _mm_andnot_ps(dummy_mask,velec);
416             velecsum         = _mm_add_ps(velecsum,velec);
417
418             fscal            = felec;
419
420             fscal            = _mm_andnot_ps(dummy_mask,fscal);
421
422              /* Update vectorial force */
423             fix3             = _mm_macc_ps(dx30,fscal,fix3);
424             fiy3             = _mm_macc_ps(dy30,fscal,fiy3);
425             fiz3             = _mm_macc_ps(dz30,fscal,fiz3);
426
427             fjx0             = _mm_macc_ps(dx30,fscal,fjx0);
428             fjy0             = _mm_macc_ps(dy30,fscal,fjy0);
429             fjz0             = _mm_macc_ps(dz30,fscal,fjz0);
430
431             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
432             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
433             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
434             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
435
436             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,fjx0,fjy0,fjz0);
437
438             /* Inner loop uses 93 flops */
439         }
440
441         /* End of innermost loop */
442
443         gmx_mm_update_iforce_3atom_swizzle_ps(fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
444                                               f+i_coord_offset+DIM,fshift+i_shift_offset);
445
446         ggid                        = gid[iidx];
447         /* Update potential energies */
448         gmx_mm_update_1pot_ps(velecsum,kernel_data->energygrp_elec+ggid);
449
450         /* Increment number of inner iterations */
451         inneriter                  += j_index_end - j_index_start;
452
453         /* Outer loop uses 19 flops */
454     }
455
456     /* Increment number of outer iterations */
457     outeriter        += nri;
458
459     /* Update outer/inner flops */
460
461     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_W4_VF,outeriter*19 + inneriter*93);
462 }
463 /*
464  * Gromacs nonbonded kernel:   nb_kernel_ElecCoul_VdwNone_GeomW4P1_F_avx_128_fma_single
465  * Electrostatics interaction: Coulomb
466  * VdW interaction:            None
467  * Geometry:                   Water4-Particle
468  * Calculate force/pot:        Force
469  */
470 void
471 nb_kernel_ElecCoul_VdwNone_GeomW4P1_F_avx_128_fma_single
472                     (t_nblist * gmx_restrict                nlist,
473                      rvec * gmx_restrict                    xx,
474                      rvec * gmx_restrict                    ff,
475                      t_forcerec * gmx_restrict              fr,
476                      t_mdatoms * gmx_restrict               mdatoms,
477                      nb_kernel_data_t * gmx_restrict        kernel_data,
478                      t_nrnb * gmx_restrict                  nrnb)
479 {
480     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
481      * just 0 for non-waters.
482      * Suffixes A,B,C,D refer to j loop unrolling done with AVX_128, e.g. for the four different
483      * jnr indices corresponding to data put in the four positions in the SIMD register.
484      */
485     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
486     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
487     int              jnrA,jnrB,jnrC,jnrD;
488     int              jnrlistA,jnrlistB,jnrlistC,jnrlistD;
489     int              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
490     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
491     real             rcutoff_scalar;
492     real             *shiftvec,*fshift,*x,*f;
493     real             *fjptrA,*fjptrB,*fjptrC,*fjptrD;
494     real             scratch[4*DIM];
495     __m128           fscal,rcutoff,rcutoff2,jidxall;
496     int              vdwioffset1;
497     __m128           ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
498     int              vdwioffset2;
499     __m128           ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
500     int              vdwioffset3;
501     __m128           ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
502     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
503     __m128           jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
504     __m128           dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
505     __m128           dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
506     __m128           dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
507     __m128           velec,felec,velecsum,facel,crf,krf,krf2;
508     real             *charge;
509     __m128           dummy_mask,cutoff_mask;
510     __m128           signbit = _mm_castsi128_ps( _mm_set1_epi32(0x80000000) );
511     __m128           one     = _mm_set1_ps(1.0);
512     __m128           two     = _mm_set1_ps(2.0);
513     x                = xx[0];
514     f                = ff[0];
515
516     nri              = nlist->nri;
517     iinr             = nlist->iinr;
518     jindex           = nlist->jindex;
519     jjnr             = nlist->jjnr;
520     shiftidx         = nlist->shift;
521     gid              = nlist->gid;
522     shiftvec         = fr->shift_vec[0];
523     fshift           = fr->fshift[0];
524     facel            = _mm_set1_ps(fr->epsfac);
525     charge           = mdatoms->chargeA;
526
527     /* Setup water-specific parameters */
528     inr              = nlist->iinr[0];
529     iq1              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+1]));
530     iq2              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+2]));
531     iq3              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+3]));
532
533     /* Avoid stupid compiler warnings */
534     jnrA = jnrB = jnrC = jnrD = 0;
535     j_coord_offsetA = 0;
536     j_coord_offsetB = 0;
537     j_coord_offsetC = 0;
538     j_coord_offsetD = 0;
539
540     outeriter        = 0;
541     inneriter        = 0;
542
543     for(iidx=0;iidx<4*DIM;iidx++)
544     {
545         scratch[iidx] = 0.0;
546     }
547
548     /* Start outer loop over neighborlists */
549     for(iidx=0; iidx<nri; iidx++)
550     {
551         /* Load shift vector for this list */
552         i_shift_offset   = DIM*shiftidx[iidx];
553
554         /* Load limits for loop over neighbors */
555         j_index_start    = jindex[iidx];
556         j_index_end      = jindex[iidx+1];
557
558         /* Get outer coordinate index */
559         inr              = iinr[iidx];
560         i_coord_offset   = DIM*inr;
561
562         /* Load i particle coords and add shift vector */
563         gmx_mm_load_shift_and_3rvec_broadcast_ps(shiftvec+i_shift_offset,x+i_coord_offset+DIM,
564                                                  &ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
565
566         fix1             = _mm_setzero_ps();
567         fiy1             = _mm_setzero_ps();
568         fiz1             = _mm_setzero_ps();
569         fix2             = _mm_setzero_ps();
570         fiy2             = _mm_setzero_ps();
571         fiz2             = _mm_setzero_ps();
572         fix3             = _mm_setzero_ps();
573         fiy3             = _mm_setzero_ps();
574         fiz3             = _mm_setzero_ps();
575
576         /* Start inner kernel loop */
577         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
578         {
579
580             /* Get j neighbor index, and coordinate index */
581             jnrA             = jjnr[jidx];
582             jnrB             = jjnr[jidx+1];
583             jnrC             = jjnr[jidx+2];
584             jnrD             = jjnr[jidx+3];
585             j_coord_offsetA  = DIM*jnrA;
586             j_coord_offsetB  = DIM*jnrB;
587             j_coord_offsetC  = DIM*jnrC;
588             j_coord_offsetD  = DIM*jnrD;
589
590             /* load j atom coordinates */
591             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
592                                               x+j_coord_offsetC,x+j_coord_offsetD,
593                                               &jx0,&jy0,&jz0);
594
595             /* Calculate displacement vector */
596             dx10             = _mm_sub_ps(ix1,jx0);
597             dy10             = _mm_sub_ps(iy1,jy0);
598             dz10             = _mm_sub_ps(iz1,jz0);
599             dx20             = _mm_sub_ps(ix2,jx0);
600             dy20             = _mm_sub_ps(iy2,jy0);
601             dz20             = _mm_sub_ps(iz2,jz0);
602             dx30             = _mm_sub_ps(ix3,jx0);
603             dy30             = _mm_sub_ps(iy3,jy0);
604             dz30             = _mm_sub_ps(iz3,jz0);
605
606             /* Calculate squared distance and things based on it */
607             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
608             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
609             rsq30            = gmx_mm_calc_rsq_ps(dx30,dy30,dz30);
610
611             rinv10           = gmx_mm_invsqrt_ps(rsq10);
612             rinv20           = gmx_mm_invsqrt_ps(rsq20);
613             rinv30           = gmx_mm_invsqrt_ps(rsq30);
614
615             rinvsq10         = _mm_mul_ps(rinv10,rinv10);
616             rinvsq20         = _mm_mul_ps(rinv20,rinv20);
617             rinvsq30         = _mm_mul_ps(rinv30,rinv30);
618
619             /* Load parameters for j particles */
620             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
621                                                               charge+jnrC+0,charge+jnrD+0);
622
623             fjx0             = _mm_setzero_ps();
624             fjy0             = _mm_setzero_ps();
625             fjz0             = _mm_setzero_ps();
626
627             /**************************
628              * CALCULATE INTERACTIONS *
629              **************************/
630
631             /* Compute parameters for interactions between i and j atoms */
632             qq10             = _mm_mul_ps(iq1,jq0);
633
634             /* COULOMB ELECTROSTATICS */
635             velec            = _mm_mul_ps(qq10,rinv10);
636             felec            = _mm_mul_ps(velec,rinvsq10);
637
638             fscal            = felec;
639
640              /* Update vectorial force */
641             fix1             = _mm_macc_ps(dx10,fscal,fix1);
642             fiy1             = _mm_macc_ps(dy10,fscal,fiy1);
643             fiz1             = _mm_macc_ps(dz10,fscal,fiz1);
644
645             fjx0             = _mm_macc_ps(dx10,fscal,fjx0);
646             fjy0             = _mm_macc_ps(dy10,fscal,fjy0);
647             fjz0             = _mm_macc_ps(dz10,fscal,fjz0);
648
649             /**************************
650              * CALCULATE INTERACTIONS *
651              **************************/
652
653             /* Compute parameters for interactions between i and j atoms */
654             qq20             = _mm_mul_ps(iq2,jq0);
655
656             /* COULOMB ELECTROSTATICS */
657             velec            = _mm_mul_ps(qq20,rinv20);
658             felec            = _mm_mul_ps(velec,rinvsq20);
659
660             fscal            = felec;
661
662              /* Update vectorial force */
663             fix2             = _mm_macc_ps(dx20,fscal,fix2);
664             fiy2             = _mm_macc_ps(dy20,fscal,fiy2);
665             fiz2             = _mm_macc_ps(dz20,fscal,fiz2);
666
667             fjx0             = _mm_macc_ps(dx20,fscal,fjx0);
668             fjy0             = _mm_macc_ps(dy20,fscal,fjy0);
669             fjz0             = _mm_macc_ps(dz20,fscal,fjz0);
670
671             /**************************
672              * CALCULATE INTERACTIONS *
673              **************************/
674
675             /* Compute parameters for interactions between i and j atoms */
676             qq30             = _mm_mul_ps(iq3,jq0);
677
678             /* COULOMB ELECTROSTATICS */
679             velec            = _mm_mul_ps(qq30,rinv30);
680             felec            = _mm_mul_ps(velec,rinvsq30);
681
682             fscal            = felec;
683
684              /* Update vectorial force */
685             fix3             = _mm_macc_ps(dx30,fscal,fix3);
686             fiy3             = _mm_macc_ps(dy30,fscal,fiy3);
687             fiz3             = _mm_macc_ps(dz30,fscal,fiz3);
688
689             fjx0             = _mm_macc_ps(dx30,fscal,fjx0);
690             fjy0             = _mm_macc_ps(dy30,fscal,fjy0);
691             fjz0             = _mm_macc_ps(dz30,fscal,fjz0);
692
693             fjptrA             = f+j_coord_offsetA;
694             fjptrB             = f+j_coord_offsetB;
695             fjptrC             = f+j_coord_offsetC;
696             fjptrD             = f+j_coord_offsetD;
697
698             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,fjx0,fjy0,fjz0);
699
700             /* Inner loop uses 90 flops */
701         }
702
703         if(jidx<j_index_end)
704         {
705
706             /* Get j neighbor index, and coordinate index */
707             jnrlistA         = jjnr[jidx];
708             jnrlistB         = jjnr[jidx+1];
709             jnrlistC         = jjnr[jidx+2];
710             jnrlistD         = jjnr[jidx+3];
711             /* Sign of each element will be negative for non-real atoms.
712              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
713              * so use it as val = _mm_andnot_ps(mask,val) to clear dummy entries.
714              */
715             dummy_mask = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
716             jnrA       = (jnrlistA>=0) ? jnrlistA : 0;
717             jnrB       = (jnrlistB>=0) ? jnrlistB : 0;
718             jnrC       = (jnrlistC>=0) ? jnrlistC : 0;
719             jnrD       = (jnrlistD>=0) ? jnrlistD : 0;
720             j_coord_offsetA  = DIM*jnrA;
721             j_coord_offsetB  = DIM*jnrB;
722             j_coord_offsetC  = DIM*jnrC;
723             j_coord_offsetD  = DIM*jnrD;
724
725             /* load j atom coordinates */
726             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
727                                               x+j_coord_offsetC,x+j_coord_offsetD,
728                                               &jx0,&jy0,&jz0);
729
730             /* Calculate displacement vector */
731             dx10             = _mm_sub_ps(ix1,jx0);
732             dy10             = _mm_sub_ps(iy1,jy0);
733             dz10             = _mm_sub_ps(iz1,jz0);
734             dx20             = _mm_sub_ps(ix2,jx0);
735             dy20             = _mm_sub_ps(iy2,jy0);
736             dz20             = _mm_sub_ps(iz2,jz0);
737             dx30             = _mm_sub_ps(ix3,jx0);
738             dy30             = _mm_sub_ps(iy3,jy0);
739             dz30             = _mm_sub_ps(iz3,jz0);
740
741             /* Calculate squared distance and things based on it */
742             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
743             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
744             rsq30            = gmx_mm_calc_rsq_ps(dx30,dy30,dz30);
745
746             rinv10           = gmx_mm_invsqrt_ps(rsq10);
747             rinv20           = gmx_mm_invsqrt_ps(rsq20);
748             rinv30           = gmx_mm_invsqrt_ps(rsq30);
749
750             rinvsq10         = _mm_mul_ps(rinv10,rinv10);
751             rinvsq20         = _mm_mul_ps(rinv20,rinv20);
752             rinvsq30         = _mm_mul_ps(rinv30,rinv30);
753
754             /* Load parameters for j particles */
755             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
756                                                               charge+jnrC+0,charge+jnrD+0);
757
758             fjx0             = _mm_setzero_ps();
759             fjy0             = _mm_setzero_ps();
760             fjz0             = _mm_setzero_ps();
761
762             /**************************
763              * CALCULATE INTERACTIONS *
764              **************************/
765
766             /* Compute parameters for interactions between i and j atoms */
767             qq10             = _mm_mul_ps(iq1,jq0);
768
769             /* COULOMB ELECTROSTATICS */
770             velec            = _mm_mul_ps(qq10,rinv10);
771             felec            = _mm_mul_ps(velec,rinvsq10);
772
773             fscal            = felec;
774
775             fscal            = _mm_andnot_ps(dummy_mask,fscal);
776
777              /* Update vectorial force */
778             fix1             = _mm_macc_ps(dx10,fscal,fix1);
779             fiy1             = _mm_macc_ps(dy10,fscal,fiy1);
780             fiz1             = _mm_macc_ps(dz10,fscal,fiz1);
781
782             fjx0             = _mm_macc_ps(dx10,fscal,fjx0);
783             fjy0             = _mm_macc_ps(dy10,fscal,fjy0);
784             fjz0             = _mm_macc_ps(dz10,fscal,fjz0);
785
786             /**************************
787              * CALCULATE INTERACTIONS *
788              **************************/
789
790             /* Compute parameters for interactions between i and j atoms */
791             qq20             = _mm_mul_ps(iq2,jq0);
792
793             /* COULOMB ELECTROSTATICS */
794             velec            = _mm_mul_ps(qq20,rinv20);
795             felec            = _mm_mul_ps(velec,rinvsq20);
796
797             fscal            = felec;
798
799             fscal            = _mm_andnot_ps(dummy_mask,fscal);
800
801              /* Update vectorial force */
802             fix2             = _mm_macc_ps(dx20,fscal,fix2);
803             fiy2             = _mm_macc_ps(dy20,fscal,fiy2);
804             fiz2             = _mm_macc_ps(dz20,fscal,fiz2);
805
806             fjx0             = _mm_macc_ps(dx20,fscal,fjx0);
807             fjy0             = _mm_macc_ps(dy20,fscal,fjy0);
808             fjz0             = _mm_macc_ps(dz20,fscal,fjz0);
809
810             /**************************
811              * CALCULATE INTERACTIONS *
812              **************************/
813
814             /* Compute parameters for interactions between i and j atoms */
815             qq30             = _mm_mul_ps(iq3,jq0);
816
817             /* COULOMB ELECTROSTATICS */
818             velec            = _mm_mul_ps(qq30,rinv30);
819             felec            = _mm_mul_ps(velec,rinvsq30);
820
821             fscal            = felec;
822
823             fscal            = _mm_andnot_ps(dummy_mask,fscal);
824
825              /* Update vectorial force */
826             fix3             = _mm_macc_ps(dx30,fscal,fix3);
827             fiy3             = _mm_macc_ps(dy30,fscal,fiy3);
828             fiz3             = _mm_macc_ps(dz30,fscal,fiz3);
829
830             fjx0             = _mm_macc_ps(dx30,fscal,fjx0);
831             fjy0             = _mm_macc_ps(dy30,fscal,fjy0);
832             fjz0             = _mm_macc_ps(dz30,fscal,fjz0);
833
834             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
835             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
836             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
837             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
838
839             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,fjx0,fjy0,fjz0);
840
841             /* Inner loop uses 90 flops */
842         }
843
844         /* End of innermost loop */
845
846         gmx_mm_update_iforce_3atom_swizzle_ps(fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
847                                               f+i_coord_offset+DIM,fshift+i_shift_offset);
848
849         /* Increment number of inner iterations */
850         inneriter                  += j_index_end - j_index_start;
851
852         /* Outer loop uses 18 flops */
853     }
854
855     /* Increment number of outer iterations */
856     outeriter        += nri;
857
858     /* Update outer/inner flops */
859
860     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_W4_F,outeriter*18 + inneriter*90);
861 }