2f8f18a44e5236f38b3db18338080d269d56c24d
[alexxy/gromacs.git] / src / gmxlib / nonbonded / nb_kernel_sse4_1_single / nb_kernel_ElecRFCut_VdwLJSh_GeomP1P1_sse4_1_single.c
1 /*
2  * Note: this file was generated by the Gromacs sse4_1_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_sse4_1_single.h"
34 #include "kernelutil_x86_sse4_1_single.h"
35
36 /*
37  * Gromacs nonbonded kernel:   nb_kernel_ElecRFCut_VdwLJSh_GeomP1P1_VF_sse4_1_single
38  * Electrostatics interaction: ReactionField
39  * VdW interaction:            LennardJones
40  * Geometry:                   Particle-Particle
41  * Calculate force/pot:        PotentialAndForce
42  */
43 void
44 nb_kernel_ElecRFCut_VdwLJSh_GeomP1P1_VF_sse4_1_single
45                     (t_nblist * gmx_restrict                nlist,
46                      rvec * gmx_restrict                    xx,
47                      rvec * gmx_restrict                    ff,
48                      t_forcerec * gmx_restrict              fr,
49                      t_mdatoms * gmx_restrict               mdatoms,
50                      nb_kernel_data_t * gmx_restrict        kernel_data,
51                      t_nrnb * gmx_restrict                  nrnb)
52 {
53     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or 
54      * just 0 for non-waters.
55      * Suffixes A,B,C,D refer to j loop unrolling done with SSE, e.g. for the four different
56      * jnr indices corresponding to data put in the four positions in the SIMD register.
57      */
58     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
59     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
60     int              jnrA,jnrB,jnrC,jnrD;
61     int              jnrlistA,jnrlistB,jnrlistC,jnrlistD;
62     int              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
63     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
64     real             rcutoff_scalar;
65     real             *shiftvec,*fshift,*x,*f;
66     real             *fjptrA,*fjptrB,*fjptrC,*fjptrD;
67     real             scratch[4*DIM];
68     __m128           tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
69     int              vdwioffset0;
70     __m128           ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
71     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
72     __m128           jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
73     __m128           dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
74     __m128           velec,felec,velecsum,facel,crf,krf,krf2;
75     real             *charge;
76     int              nvdwtype;
77     __m128           rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
78     int              *vdwtype;
79     real             *vdwparam;
80     __m128           one_sixth   = _mm_set1_ps(1.0/6.0);
81     __m128           one_twelfth = _mm_set1_ps(1.0/12.0);
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     krf              = _mm_set1_ps(fr->ic->k_rf);
100     krf2             = _mm_set1_ps(fr->ic->k_rf*2.0);
101     crf              = _mm_set1_ps(fr->ic->c_rf);
102     nvdwtype         = fr->ntype;
103     vdwparam         = fr->nbfp;
104     vdwtype          = mdatoms->typeA;
105
106     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
107     rcutoff_scalar   = fr->rcoulomb;
108     rcutoff          = _mm_set1_ps(rcutoff_scalar);
109     rcutoff2         = _mm_mul_ps(rcutoff,rcutoff);
110
111     sh_vdw_invrcut6  = _mm_set1_ps(fr->ic->sh_invrc6);
112     rvdw             = _mm_set1_ps(fr->rvdw);
113
114     /* Avoid stupid compiler warnings */
115     jnrA = jnrB = jnrC = jnrD = 0;
116     j_coord_offsetA = 0;
117     j_coord_offsetB = 0;
118     j_coord_offsetC = 0;
119     j_coord_offsetD = 0;
120
121     outeriter        = 0;
122     inneriter        = 0;
123
124     for(iidx=0;iidx<4*DIM;iidx++)
125     {
126         scratch[iidx] = 0.0;
127     }
128
129     /* Start outer loop over neighborlists */
130     for(iidx=0; iidx<nri; iidx++)
131     {
132         /* Load shift vector for this list */
133         i_shift_offset   = DIM*shiftidx[iidx];
134
135         /* Load limits for loop over neighbors */
136         j_index_start    = jindex[iidx];
137         j_index_end      = jindex[iidx+1];
138
139         /* Get outer coordinate index */
140         inr              = iinr[iidx];
141         i_coord_offset   = DIM*inr;
142
143         /* Load i particle coords and add shift vector */
144         gmx_mm_load_shift_and_1rvec_broadcast_ps(shiftvec+i_shift_offset,x+i_coord_offset,&ix0,&iy0,&iz0);
145
146         fix0             = _mm_setzero_ps();
147         fiy0             = _mm_setzero_ps();
148         fiz0             = _mm_setzero_ps();
149
150         /* Load parameters for i particles */
151         iq0              = _mm_mul_ps(facel,_mm_load1_ps(charge+inr+0));
152         vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
153
154         /* Reset potential sums */
155         velecsum         = _mm_setzero_ps();
156         vvdwsum          = _mm_setzero_ps();
157
158         /* Start inner kernel loop */
159         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
160         {
161
162             /* Get j neighbor index, and coordinate index */
163             jnrA             = jjnr[jidx];
164             jnrB             = jjnr[jidx+1];
165             jnrC             = jjnr[jidx+2];
166             jnrD             = jjnr[jidx+3];
167             j_coord_offsetA  = DIM*jnrA;
168             j_coord_offsetB  = DIM*jnrB;
169             j_coord_offsetC  = DIM*jnrC;
170             j_coord_offsetD  = DIM*jnrD;
171
172             /* load j atom coordinates */
173             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
174                                               x+j_coord_offsetC,x+j_coord_offsetD,
175                                               &jx0,&jy0,&jz0);
176
177             /* Calculate displacement vector */
178             dx00             = _mm_sub_ps(ix0,jx0);
179             dy00             = _mm_sub_ps(iy0,jy0);
180             dz00             = _mm_sub_ps(iz0,jz0);
181
182             /* Calculate squared distance and things based on it */
183             rsq00            = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
184
185             rinv00           = gmx_mm_invsqrt_ps(rsq00);
186
187             rinvsq00         = _mm_mul_ps(rinv00,rinv00);
188
189             /* Load parameters for j particles */
190             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
191                                                               charge+jnrC+0,charge+jnrD+0);
192             vdwjidx0A        = 2*vdwtype[jnrA+0];
193             vdwjidx0B        = 2*vdwtype[jnrB+0];
194             vdwjidx0C        = 2*vdwtype[jnrC+0];
195             vdwjidx0D        = 2*vdwtype[jnrD+0];
196
197             /**************************
198              * CALCULATE INTERACTIONS *
199              **************************/
200
201             if (gmx_mm_any_lt(rsq00,rcutoff2))
202             {
203
204             /* Compute parameters for interactions between i and j atoms */
205             qq00             = _mm_mul_ps(iq0,jq0);
206             gmx_mm_load_4pair_swizzle_ps(vdwparam+vdwioffset0+vdwjidx0A,
207                                          vdwparam+vdwioffset0+vdwjidx0B,
208                                          vdwparam+vdwioffset0+vdwjidx0C,
209                                          vdwparam+vdwioffset0+vdwjidx0D,
210                                          &c6_00,&c12_00);
211
212             /* REACTION-FIELD ELECTROSTATICS */
213             velec            = _mm_mul_ps(qq00,_mm_sub_ps(_mm_add_ps(rinv00,_mm_mul_ps(krf,rsq00)),crf));
214             felec            = _mm_mul_ps(qq00,_mm_sub_ps(_mm_mul_ps(rinv00,rinvsq00),krf2));
215
216             /* LENNARD-JONES DISPERSION/REPULSION */
217
218             rinvsix          = _mm_mul_ps(_mm_mul_ps(rinvsq00,rinvsq00),rinvsq00);
219             vvdw6            = _mm_mul_ps(c6_00,rinvsix);
220             vvdw12           = _mm_mul_ps(c12_00,_mm_mul_ps(rinvsix,rinvsix));
221             vvdw             = _mm_sub_ps(_mm_mul_ps( _mm_sub_ps(vvdw12 , _mm_mul_ps(c12_00,_mm_mul_ps(sh_vdw_invrcut6,sh_vdw_invrcut6))), one_twelfth) ,
222                                           _mm_mul_ps( _mm_sub_ps(vvdw6,_mm_mul_ps(c6_00,sh_vdw_invrcut6)),one_sixth));
223             fvdw             = _mm_mul_ps(_mm_sub_ps(vvdw12,vvdw6),rinvsq00);
224
225             cutoff_mask      = _mm_cmplt_ps(rsq00,rcutoff2);
226
227             /* Update potential sum for this i atom from the interaction with this j atom. */
228             velec            = _mm_and_ps(velec,cutoff_mask);
229             velecsum         = _mm_add_ps(velecsum,velec);
230             vvdw             = _mm_and_ps(vvdw,cutoff_mask);
231             vvdwsum          = _mm_add_ps(vvdwsum,vvdw);
232
233             fscal            = _mm_add_ps(felec,fvdw);
234
235             fscal            = _mm_and_ps(fscal,cutoff_mask);
236
237             /* Calculate temporary vectorial force */
238             tx               = _mm_mul_ps(fscal,dx00);
239             ty               = _mm_mul_ps(fscal,dy00);
240             tz               = _mm_mul_ps(fscal,dz00);
241
242             /* Update vectorial force */
243             fix0             = _mm_add_ps(fix0,tx);
244             fiy0             = _mm_add_ps(fiy0,ty);
245             fiz0             = _mm_add_ps(fiz0,tz);
246
247             fjptrA             = f+j_coord_offsetA;
248             fjptrB             = f+j_coord_offsetB;
249             fjptrC             = f+j_coord_offsetC;
250             fjptrD             = f+j_coord_offsetD;
251             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
252
253             }
254
255             /* Inner loop uses 54 flops */
256         }
257
258         if(jidx<j_index_end)
259         {
260
261             /* Get j neighbor index, and coordinate index */
262             jnrlistA         = jjnr[jidx];
263             jnrlistB         = jjnr[jidx+1];
264             jnrlistC         = jjnr[jidx+2];
265             jnrlistD         = jjnr[jidx+3];
266             /* Sign of each element will be negative for non-real atoms.
267              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
268              * so use it as val = _mm_andnot_ps(mask,val) to clear dummy entries.
269              */
270             dummy_mask = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
271             jnrA       = (jnrlistA>=0) ? jnrlistA : 0;
272             jnrB       = (jnrlistB>=0) ? jnrlistB : 0;
273             jnrC       = (jnrlistC>=0) ? jnrlistC : 0;
274             jnrD       = (jnrlistD>=0) ? jnrlistD : 0;
275             j_coord_offsetA  = DIM*jnrA;
276             j_coord_offsetB  = DIM*jnrB;
277             j_coord_offsetC  = DIM*jnrC;
278             j_coord_offsetD  = DIM*jnrD;
279
280             /* load j atom coordinates */
281             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
282                                               x+j_coord_offsetC,x+j_coord_offsetD,
283                                               &jx0,&jy0,&jz0);
284
285             /* Calculate displacement vector */
286             dx00             = _mm_sub_ps(ix0,jx0);
287             dy00             = _mm_sub_ps(iy0,jy0);
288             dz00             = _mm_sub_ps(iz0,jz0);
289
290             /* Calculate squared distance and things based on it */
291             rsq00            = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
292
293             rinv00           = gmx_mm_invsqrt_ps(rsq00);
294
295             rinvsq00         = _mm_mul_ps(rinv00,rinv00);
296
297             /* Load parameters for j particles */
298             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
299                                                               charge+jnrC+0,charge+jnrD+0);
300             vdwjidx0A        = 2*vdwtype[jnrA+0];
301             vdwjidx0B        = 2*vdwtype[jnrB+0];
302             vdwjidx0C        = 2*vdwtype[jnrC+0];
303             vdwjidx0D        = 2*vdwtype[jnrD+0];
304
305             /**************************
306              * CALCULATE INTERACTIONS *
307              **************************/
308
309             if (gmx_mm_any_lt(rsq00,rcutoff2))
310             {
311
312             /* Compute parameters for interactions between i and j atoms */
313             qq00             = _mm_mul_ps(iq0,jq0);
314             gmx_mm_load_4pair_swizzle_ps(vdwparam+vdwioffset0+vdwjidx0A,
315                                          vdwparam+vdwioffset0+vdwjidx0B,
316                                          vdwparam+vdwioffset0+vdwjidx0C,
317                                          vdwparam+vdwioffset0+vdwjidx0D,
318                                          &c6_00,&c12_00);
319
320             /* REACTION-FIELD ELECTROSTATICS */
321             velec            = _mm_mul_ps(qq00,_mm_sub_ps(_mm_add_ps(rinv00,_mm_mul_ps(krf,rsq00)),crf));
322             felec            = _mm_mul_ps(qq00,_mm_sub_ps(_mm_mul_ps(rinv00,rinvsq00),krf2));
323
324             /* LENNARD-JONES DISPERSION/REPULSION */
325
326             rinvsix          = _mm_mul_ps(_mm_mul_ps(rinvsq00,rinvsq00),rinvsq00);
327             vvdw6            = _mm_mul_ps(c6_00,rinvsix);
328             vvdw12           = _mm_mul_ps(c12_00,_mm_mul_ps(rinvsix,rinvsix));
329             vvdw             = _mm_sub_ps(_mm_mul_ps( _mm_sub_ps(vvdw12 , _mm_mul_ps(c12_00,_mm_mul_ps(sh_vdw_invrcut6,sh_vdw_invrcut6))), one_twelfth) ,
330                                           _mm_mul_ps( _mm_sub_ps(vvdw6,_mm_mul_ps(c6_00,sh_vdw_invrcut6)),one_sixth));
331             fvdw             = _mm_mul_ps(_mm_sub_ps(vvdw12,vvdw6),rinvsq00);
332
333             cutoff_mask      = _mm_cmplt_ps(rsq00,rcutoff2);
334
335             /* Update potential sum for this i atom from the interaction with this j atom. */
336             velec            = _mm_and_ps(velec,cutoff_mask);
337             velec            = _mm_andnot_ps(dummy_mask,velec);
338             velecsum         = _mm_add_ps(velecsum,velec);
339             vvdw             = _mm_and_ps(vvdw,cutoff_mask);
340             vvdw             = _mm_andnot_ps(dummy_mask,vvdw);
341             vvdwsum          = _mm_add_ps(vvdwsum,vvdw);
342
343             fscal            = _mm_add_ps(felec,fvdw);
344
345             fscal            = _mm_and_ps(fscal,cutoff_mask);
346
347             fscal            = _mm_andnot_ps(dummy_mask,fscal);
348
349             /* Calculate temporary vectorial force */
350             tx               = _mm_mul_ps(fscal,dx00);
351             ty               = _mm_mul_ps(fscal,dy00);
352             tz               = _mm_mul_ps(fscal,dz00);
353
354             /* Update vectorial force */
355             fix0             = _mm_add_ps(fix0,tx);
356             fiy0             = _mm_add_ps(fiy0,ty);
357             fiz0             = _mm_add_ps(fiz0,tz);
358
359             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
360             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
361             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
362             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
363             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
364
365             }
366
367             /* Inner loop uses 54 flops */
368         }
369
370         /* End of innermost loop */
371
372         gmx_mm_update_iforce_1atom_swizzle_ps(fix0,fiy0,fiz0,
373                                               f+i_coord_offset,fshift+i_shift_offset);
374
375         ggid                        = gid[iidx];
376         /* Update potential energies */
377         gmx_mm_update_1pot_ps(velecsum,kernel_data->energygrp_elec+ggid);
378         gmx_mm_update_1pot_ps(vvdwsum,kernel_data->energygrp_vdw+ggid);
379
380         /* Increment number of inner iterations */
381         inneriter                  += j_index_end - j_index_start;
382
383         /* Outer loop uses 9 flops */
384     }
385
386     /* Increment number of outer iterations */
387     outeriter        += nri;
388
389     /* Update outer/inner flops */
390
391     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_VF,outeriter*9 + inneriter*54);
392 }
393 /*
394  * Gromacs nonbonded kernel:   nb_kernel_ElecRFCut_VdwLJSh_GeomP1P1_F_sse4_1_single
395  * Electrostatics interaction: ReactionField
396  * VdW interaction:            LennardJones
397  * Geometry:                   Particle-Particle
398  * Calculate force/pot:        Force
399  */
400 void
401 nb_kernel_ElecRFCut_VdwLJSh_GeomP1P1_F_sse4_1_single
402                     (t_nblist * gmx_restrict                nlist,
403                      rvec * gmx_restrict                    xx,
404                      rvec * gmx_restrict                    ff,
405                      t_forcerec * gmx_restrict              fr,
406                      t_mdatoms * gmx_restrict               mdatoms,
407                      nb_kernel_data_t * gmx_restrict        kernel_data,
408                      t_nrnb * gmx_restrict                  nrnb)
409 {
410     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or 
411      * just 0 for non-waters.
412      * Suffixes A,B,C,D refer to j loop unrolling done with SSE, e.g. for the four different
413      * jnr indices corresponding to data put in the four positions in the SIMD register.
414      */
415     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
416     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
417     int              jnrA,jnrB,jnrC,jnrD;
418     int              jnrlistA,jnrlistB,jnrlistC,jnrlistD;
419     int              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
420     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
421     real             rcutoff_scalar;
422     real             *shiftvec,*fshift,*x,*f;
423     real             *fjptrA,*fjptrB,*fjptrC,*fjptrD;
424     real             scratch[4*DIM];
425     __m128           tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
426     int              vdwioffset0;
427     __m128           ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
428     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
429     __m128           jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
430     __m128           dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
431     __m128           velec,felec,velecsum,facel,crf,krf,krf2;
432     real             *charge;
433     int              nvdwtype;
434     __m128           rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
435     int              *vdwtype;
436     real             *vdwparam;
437     __m128           one_sixth   = _mm_set1_ps(1.0/6.0);
438     __m128           one_twelfth = _mm_set1_ps(1.0/12.0);
439     __m128           dummy_mask,cutoff_mask;
440     __m128           signbit = _mm_castsi128_ps( _mm_set1_epi32(0x80000000) );
441     __m128           one     = _mm_set1_ps(1.0);
442     __m128           two     = _mm_set1_ps(2.0);
443     x                = xx[0];
444     f                = ff[0];
445
446     nri              = nlist->nri;
447     iinr             = nlist->iinr;
448     jindex           = nlist->jindex;
449     jjnr             = nlist->jjnr;
450     shiftidx         = nlist->shift;
451     gid              = nlist->gid;
452     shiftvec         = fr->shift_vec[0];
453     fshift           = fr->fshift[0];
454     facel            = _mm_set1_ps(fr->epsfac);
455     charge           = mdatoms->chargeA;
456     krf              = _mm_set1_ps(fr->ic->k_rf);
457     krf2             = _mm_set1_ps(fr->ic->k_rf*2.0);
458     crf              = _mm_set1_ps(fr->ic->c_rf);
459     nvdwtype         = fr->ntype;
460     vdwparam         = fr->nbfp;
461     vdwtype          = mdatoms->typeA;
462
463     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
464     rcutoff_scalar   = fr->rcoulomb;
465     rcutoff          = _mm_set1_ps(rcutoff_scalar);
466     rcutoff2         = _mm_mul_ps(rcutoff,rcutoff);
467
468     sh_vdw_invrcut6  = _mm_set1_ps(fr->ic->sh_invrc6);
469     rvdw             = _mm_set1_ps(fr->rvdw);
470
471     /* Avoid stupid compiler warnings */
472     jnrA = jnrB = jnrC = jnrD = 0;
473     j_coord_offsetA = 0;
474     j_coord_offsetB = 0;
475     j_coord_offsetC = 0;
476     j_coord_offsetD = 0;
477
478     outeriter        = 0;
479     inneriter        = 0;
480
481     for(iidx=0;iidx<4*DIM;iidx++)
482     {
483         scratch[iidx] = 0.0;
484     }
485
486     /* Start outer loop over neighborlists */
487     for(iidx=0; iidx<nri; iidx++)
488     {
489         /* Load shift vector for this list */
490         i_shift_offset   = DIM*shiftidx[iidx];
491
492         /* Load limits for loop over neighbors */
493         j_index_start    = jindex[iidx];
494         j_index_end      = jindex[iidx+1];
495
496         /* Get outer coordinate index */
497         inr              = iinr[iidx];
498         i_coord_offset   = DIM*inr;
499
500         /* Load i particle coords and add shift vector */
501         gmx_mm_load_shift_and_1rvec_broadcast_ps(shiftvec+i_shift_offset,x+i_coord_offset,&ix0,&iy0,&iz0);
502
503         fix0             = _mm_setzero_ps();
504         fiy0             = _mm_setzero_ps();
505         fiz0             = _mm_setzero_ps();
506
507         /* Load parameters for i particles */
508         iq0              = _mm_mul_ps(facel,_mm_load1_ps(charge+inr+0));
509         vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
510
511         /* Start inner kernel loop */
512         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
513         {
514
515             /* Get j neighbor index, and coordinate index */
516             jnrA             = jjnr[jidx];
517             jnrB             = jjnr[jidx+1];
518             jnrC             = jjnr[jidx+2];
519             jnrD             = jjnr[jidx+3];
520             j_coord_offsetA  = DIM*jnrA;
521             j_coord_offsetB  = DIM*jnrB;
522             j_coord_offsetC  = DIM*jnrC;
523             j_coord_offsetD  = DIM*jnrD;
524
525             /* load j atom coordinates */
526             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
527                                               x+j_coord_offsetC,x+j_coord_offsetD,
528                                               &jx0,&jy0,&jz0);
529
530             /* Calculate displacement vector */
531             dx00             = _mm_sub_ps(ix0,jx0);
532             dy00             = _mm_sub_ps(iy0,jy0);
533             dz00             = _mm_sub_ps(iz0,jz0);
534
535             /* Calculate squared distance and things based on it */
536             rsq00            = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
537
538             rinv00           = gmx_mm_invsqrt_ps(rsq00);
539
540             rinvsq00         = _mm_mul_ps(rinv00,rinv00);
541
542             /* Load parameters for j particles */
543             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
544                                                               charge+jnrC+0,charge+jnrD+0);
545             vdwjidx0A        = 2*vdwtype[jnrA+0];
546             vdwjidx0B        = 2*vdwtype[jnrB+0];
547             vdwjidx0C        = 2*vdwtype[jnrC+0];
548             vdwjidx0D        = 2*vdwtype[jnrD+0];
549
550             /**************************
551              * CALCULATE INTERACTIONS *
552              **************************/
553
554             if (gmx_mm_any_lt(rsq00,rcutoff2))
555             {
556
557             /* Compute parameters for interactions between i and j atoms */
558             qq00             = _mm_mul_ps(iq0,jq0);
559             gmx_mm_load_4pair_swizzle_ps(vdwparam+vdwioffset0+vdwjidx0A,
560                                          vdwparam+vdwioffset0+vdwjidx0B,
561                                          vdwparam+vdwioffset0+vdwjidx0C,
562                                          vdwparam+vdwioffset0+vdwjidx0D,
563                                          &c6_00,&c12_00);
564
565             /* REACTION-FIELD ELECTROSTATICS */
566             felec            = _mm_mul_ps(qq00,_mm_sub_ps(_mm_mul_ps(rinv00,rinvsq00),krf2));
567
568             /* LENNARD-JONES DISPERSION/REPULSION */
569
570             rinvsix          = _mm_mul_ps(_mm_mul_ps(rinvsq00,rinvsq00),rinvsq00);
571             fvdw             = _mm_mul_ps(_mm_sub_ps(_mm_mul_ps(c12_00,rinvsix),c6_00),_mm_mul_ps(rinvsix,rinvsq00));
572
573             cutoff_mask      = _mm_cmplt_ps(rsq00,rcutoff2);
574
575             fscal            = _mm_add_ps(felec,fvdw);
576
577             fscal            = _mm_and_ps(fscal,cutoff_mask);
578
579             /* Calculate temporary vectorial force */
580             tx               = _mm_mul_ps(fscal,dx00);
581             ty               = _mm_mul_ps(fscal,dy00);
582             tz               = _mm_mul_ps(fscal,dz00);
583
584             /* Update vectorial force */
585             fix0             = _mm_add_ps(fix0,tx);
586             fiy0             = _mm_add_ps(fiy0,ty);
587             fiz0             = _mm_add_ps(fiz0,tz);
588
589             fjptrA             = f+j_coord_offsetA;
590             fjptrB             = f+j_coord_offsetB;
591             fjptrC             = f+j_coord_offsetC;
592             fjptrD             = f+j_coord_offsetD;
593             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
594
595             }
596
597             /* Inner loop uses 37 flops */
598         }
599
600         if(jidx<j_index_end)
601         {
602
603             /* Get j neighbor index, and coordinate index */
604             jnrlistA         = jjnr[jidx];
605             jnrlistB         = jjnr[jidx+1];
606             jnrlistC         = jjnr[jidx+2];
607             jnrlistD         = jjnr[jidx+3];
608             /* Sign of each element will be negative for non-real atoms.
609              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
610              * so use it as val = _mm_andnot_ps(mask,val) to clear dummy entries.
611              */
612             dummy_mask = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
613             jnrA       = (jnrlistA>=0) ? jnrlistA : 0;
614             jnrB       = (jnrlistB>=0) ? jnrlistB : 0;
615             jnrC       = (jnrlistC>=0) ? jnrlistC : 0;
616             jnrD       = (jnrlistD>=0) ? jnrlistD : 0;
617             j_coord_offsetA  = DIM*jnrA;
618             j_coord_offsetB  = DIM*jnrB;
619             j_coord_offsetC  = DIM*jnrC;
620             j_coord_offsetD  = DIM*jnrD;
621
622             /* load j atom coordinates */
623             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
624                                               x+j_coord_offsetC,x+j_coord_offsetD,
625                                               &jx0,&jy0,&jz0);
626
627             /* Calculate displacement vector */
628             dx00             = _mm_sub_ps(ix0,jx0);
629             dy00             = _mm_sub_ps(iy0,jy0);
630             dz00             = _mm_sub_ps(iz0,jz0);
631
632             /* Calculate squared distance and things based on it */
633             rsq00            = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
634
635             rinv00           = gmx_mm_invsqrt_ps(rsq00);
636
637             rinvsq00         = _mm_mul_ps(rinv00,rinv00);
638
639             /* Load parameters for j particles */
640             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
641                                                               charge+jnrC+0,charge+jnrD+0);
642             vdwjidx0A        = 2*vdwtype[jnrA+0];
643             vdwjidx0B        = 2*vdwtype[jnrB+0];
644             vdwjidx0C        = 2*vdwtype[jnrC+0];
645             vdwjidx0D        = 2*vdwtype[jnrD+0];
646
647             /**************************
648              * CALCULATE INTERACTIONS *
649              **************************/
650
651             if (gmx_mm_any_lt(rsq00,rcutoff2))
652             {
653
654             /* Compute parameters for interactions between i and j atoms */
655             qq00             = _mm_mul_ps(iq0,jq0);
656             gmx_mm_load_4pair_swizzle_ps(vdwparam+vdwioffset0+vdwjidx0A,
657                                          vdwparam+vdwioffset0+vdwjidx0B,
658                                          vdwparam+vdwioffset0+vdwjidx0C,
659                                          vdwparam+vdwioffset0+vdwjidx0D,
660                                          &c6_00,&c12_00);
661
662             /* REACTION-FIELD ELECTROSTATICS */
663             felec            = _mm_mul_ps(qq00,_mm_sub_ps(_mm_mul_ps(rinv00,rinvsq00),krf2));
664
665             /* LENNARD-JONES DISPERSION/REPULSION */
666
667             rinvsix          = _mm_mul_ps(_mm_mul_ps(rinvsq00,rinvsq00),rinvsq00);
668             fvdw             = _mm_mul_ps(_mm_sub_ps(_mm_mul_ps(c12_00,rinvsix),c6_00),_mm_mul_ps(rinvsix,rinvsq00));
669
670             cutoff_mask      = _mm_cmplt_ps(rsq00,rcutoff2);
671
672             fscal            = _mm_add_ps(felec,fvdw);
673
674             fscal            = _mm_and_ps(fscal,cutoff_mask);
675
676             fscal            = _mm_andnot_ps(dummy_mask,fscal);
677
678             /* Calculate temporary vectorial force */
679             tx               = _mm_mul_ps(fscal,dx00);
680             ty               = _mm_mul_ps(fscal,dy00);
681             tz               = _mm_mul_ps(fscal,dz00);
682
683             /* Update vectorial force */
684             fix0             = _mm_add_ps(fix0,tx);
685             fiy0             = _mm_add_ps(fiy0,ty);
686             fiz0             = _mm_add_ps(fiz0,tz);
687
688             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
689             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
690             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
691             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
692             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,tx,ty,tz);
693
694             }
695
696             /* Inner loop uses 37 flops */
697         }
698
699         /* End of innermost loop */
700
701         gmx_mm_update_iforce_1atom_swizzle_ps(fix0,fiy0,fiz0,
702                                               f+i_coord_offset,fshift+i_shift_offset);
703
704         /* Increment number of inner iterations */
705         inneriter                  += j_index_end - j_index_start;
706
707         /* Outer loop uses 7 flops */
708     }
709
710     /* Increment number of outer iterations */
711     outeriter        += nri;
712
713     /* Update outer/inner flops */
714
715     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_F,outeriter*7 + inneriter*37);
716 }