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