Fujitsu Sparc64 acceleration and general fixes for non-x86 builds
[alexxy/gromacs.git] / src / gmxlib / nonbonded / nb_kernel_sparc64_hpc_ace_double / nb_kernel_ElecRFCut_VdwLJSw_GeomW4P1_sparc64_hpc_ace_double.c
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2012, by the GROMACS development team, led by
5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
6  * others, as listed in the AUTHORS file in the top-level source
7  * directory and at http://www.gromacs.org.
8  *
9  * GROMACS is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public License
11  * as published by the Free Software Foundation; either version 2.1
12  * of the License, or (at your option) any later version.
13  *
14  * GROMACS is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17  * Lesser General Public License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with GROMACS; if not, see
21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
23  *
24  * If you want to redistribute modifications to GROMACS, please
25  * consider that scientific software is very special. Version
26  * control is crucial - bugs must be traceable. We will be happy to
27  * consider code for inclusion in the official distribution, but
28  * derived work must not be called official GROMACS. Details are found
29  * in the README & COPYING files - if they are missing, get the
30  * official version at http://www.gromacs.org.
31  *
32  * To help us fund GROMACS development, we humbly ask that you cite
33  * the research papers on the package. Check out http://www.gromacs.org.
34  */
35 /*
36  * Note: this file was generated by the GROMACS sparc64_hpc_ace_double kernel generator.
37  */
38 #ifdef HAVE_CONFIG_H
39 #include <config.h>
40 #endif
41
42 #include <math.h>
43
44 #include "../nb_kernel.h"
45 #include "types/simple.h"
46 #include "vec.h"
47 #include "nrnb.h"
48
49 #include "kernelutil_sparc64_hpc_ace_double.h"
50
51 /*
52  * Gromacs nonbonded kernel:   nb_kernel_ElecRFCut_VdwLJSw_GeomW4P1_VF_sparc64_hpc_ace_double
53  * Electrostatics interaction: ReactionField
54  * VdW interaction:            LennardJones
55  * Geometry:                   Water4-Particle
56  * Calculate force/pot:        PotentialAndForce
57  */
58 void
59 nb_kernel_ElecRFCut_VdwLJSw_GeomW4P1_VF_sparc64_hpc_ace_double
60                     (t_nblist * gmx_restrict                nlist,
61                      rvec * gmx_restrict                    xx,
62                      rvec * gmx_restrict                    ff,
63                      t_forcerec * gmx_restrict              fr,
64                      t_mdatoms * gmx_restrict               mdatoms,
65                      nb_kernel_data_t * gmx_restrict        kernel_data,
66                      t_nrnb * gmx_restrict                  nrnb)
67 {
68     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
69      * just 0 for non-waters.
70      * Suffixes A,B refer to j loop unrolling done with double precision SIMD, e.g. for the two different
71      * jnr indices corresponding to data put in the four positions in the SIMD register.
72      */
73     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
74     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
75     int              jnrA,jnrB;
76     int              j_coord_offsetA,j_coord_offsetB;
77     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
78     real             rcutoff_scalar;
79     real             *shiftvec,*fshift,*x,*f;
80     _fjsp_v2r8       tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
81     int              vdwioffset0;
82     _fjsp_v2r8       ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
83     int              vdwioffset1;
84     _fjsp_v2r8       ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
85     int              vdwioffset2;
86     _fjsp_v2r8       ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
87     int              vdwioffset3;
88     _fjsp_v2r8       ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
89     int              vdwjidx0A,vdwjidx0B;
90     _fjsp_v2r8       jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
91     _fjsp_v2r8       dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
92     _fjsp_v2r8       dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
93     _fjsp_v2r8       dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
94     _fjsp_v2r8       dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
95     _fjsp_v2r8       velec,felec,velecsum,facel,crf,krf,krf2;
96     real             *charge;
97     int              nvdwtype;
98     _fjsp_v2r8       rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
99     int              *vdwtype;
100     real             *vdwparam;
101     _fjsp_v2r8       one_sixth   = gmx_fjsp_set1_v2r8(1.0/6.0);
102     _fjsp_v2r8       one_twelfth = gmx_fjsp_set1_v2r8(1.0/12.0);
103     _fjsp_v2r8       rswitch,swV3,swV4,swV5,swF2,swF3,swF4,d,d2,sw,dsw;
104     real             rswitch_scalar,d_scalar;
105     _fjsp_v2r8       itab_tmp;
106     _fjsp_v2r8       dummy_mask,cutoff_mask;
107     _fjsp_v2r8       one     = gmx_fjsp_set1_v2r8(1.0);
108     _fjsp_v2r8       two     = gmx_fjsp_set1_v2r8(2.0);
109     union { _fjsp_v2r8 simd; long long int i[2]; } vfconv,gbconv,ewconv;
110
111     x                = xx[0];
112     f                = ff[0];
113
114     nri              = nlist->nri;
115     iinr             = nlist->iinr;
116     jindex           = nlist->jindex;
117     jjnr             = nlist->jjnr;
118     shiftidx         = nlist->shift;
119     gid              = nlist->gid;
120     shiftvec         = fr->shift_vec[0];
121     fshift           = fr->fshift[0];
122     facel            = gmx_fjsp_set1_v2r8(fr->epsfac);
123     charge           = mdatoms->chargeA;
124     krf              = gmx_fjsp_set1_v2r8(fr->ic->k_rf);
125     krf2             = gmx_fjsp_set1_v2r8(fr->ic->k_rf*2.0);
126     crf              = gmx_fjsp_set1_v2r8(fr->ic->c_rf);
127     nvdwtype         = fr->ntype;
128     vdwparam         = fr->nbfp;
129     vdwtype          = mdatoms->typeA;
130
131     /* Setup water-specific parameters */
132     inr              = nlist->iinr[0];
133     iq1              = _fjsp_mul_v2r8(facel,gmx_fjsp_set1_v2r8(charge[inr+1]));
134     iq2              = _fjsp_mul_v2r8(facel,gmx_fjsp_set1_v2r8(charge[inr+2]));
135     iq3              = _fjsp_mul_v2r8(facel,gmx_fjsp_set1_v2r8(charge[inr+3]));
136     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
137
138     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
139     rcutoff_scalar   = fr->rcoulomb;
140     rcutoff          = gmx_fjsp_set1_v2r8(rcutoff_scalar);
141     rcutoff2         = _fjsp_mul_v2r8(rcutoff,rcutoff);
142
143     rswitch_scalar   = fr->rvdw_switch;
144     rswitch          = gmx_fjsp_set1_v2r8(rswitch_scalar);
145     /* Setup switch parameters */
146     d_scalar         = rcutoff_scalar-rswitch_scalar;
147     d                = gmx_fjsp_set1_v2r8(d_scalar);
148     swV3             = gmx_fjsp_set1_v2r8(-10.0/(d_scalar*d_scalar*d_scalar));
149     swV4             = gmx_fjsp_set1_v2r8( 15.0/(d_scalar*d_scalar*d_scalar*d_scalar));
150     swV5             = gmx_fjsp_set1_v2r8( -6.0/(d_scalar*d_scalar*d_scalar*d_scalar*d_scalar));
151     swF2             = gmx_fjsp_set1_v2r8(-30.0/(d_scalar*d_scalar*d_scalar));
152     swF3             = gmx_fjsp_set1_v2r8( 60.0/(d_scalar*d_scalar*d_scalar*d_scalar));
153     swF4             = gmx_fjsp_set1_v2r8(-30.0/(d_scalar*d_scalar*d_scalar*d_scalar*d_scalar));
154
155     /* Avoid stupid compiler warnings */
156     jnrA = jnrB = 0;
157     j_coord_offsetA = 0;
158     j_coord_offsetB = 0;
159
160     outeriter        = 0;
161     inneriter        = 0;
162
163     /* Start outer loop over neighborlists */
164     for(iidx=0; iidx<nri; iidx++)
165     {
166         /* Load shift vector for this list */
167         i_shift_offset   = DIM*shiftidx[iidx];
168
169         /* Load limits for loop over neighbors */
170         j_index_start    = jindex[iidx];
171         j_index_end      = jindex[iidx+1];
172
173         /* Get outer coordinate index */
174         inr              = iinr[iidx];
175         i_coord_offset   = DIM*inr;
176
177         /* Load i particle coords and add shift vector */
178         gmx_fjsp_load_shift_and_4rvec_broadcast_v2r8(shiftvec+i_shift_offset,x+i_coord_offset,
179                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
180
181         fix0             = _fjsp_setzero_v2r8();
182         fiy0             = _fjsp_setzero_v2r8();
183         fiz0             = _fjsp_setzero_v2r8();
184         fix1             = _fjsp_setzero_v2r8();
185         fiy1             = _fjsp_setzero_v2r8();
186         fiz1             = _fjsp_setzero_v2r8();
187         fix2             = _fjsp_setzero_v2r8();
188         fiy2             = _fjsp_setzero_v2r8();
189         fiz2             = _fjsp_setzero_v2r8();
190         fix3             = _fjsp_setzero_v2r8();
191         fiy3             = _fjsp_setzero_v2r8();
192         fiz3             = _fjsp_setzero_v2r8();
193
194         /* Reset potential sums */
195         velecsum         = _fjsp_setzero_v2r8();
196         vvdwsum          = _fjsp_setzero_v2r8();
197
198         /* Start inner kernel loop */
199         for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
200         {
201
202             /* Get j neighbor index, and coordinate index */
203             jnrA             = jjnr[jidx];
204             jnrB             = jjnr[jidx+1];
205             j_coord_offsetA  = DIM*jnrA;
206             j_coord_offsetB  = DIM*jnrB;
207
208             /* load j atom coordinates */
209             gmx_fjsp_load_1rvec_2ptr_swizzle_v2r8(x+j_coord_offsetA,x+j_coord_offsetB,
210                                               &jx0,&jy0,&jz0);
211
212             /* Calculate displacement vector */
213             dx00             = _fjsp_sub_v2r8(ix0,jx0);
214             dy00             = _fjsp_sub_v2r8(iy0,jy0);
215             dz00             = _fjsp_sub_v2r8(iz0,jz0);
216             dx10             = _fjsp_sub_v2r8(ix1,jx0);
217             dy10             = _fjsp_sub_v2r8(iy1,jy0);
218             dz10             = _fjsp_sub_v2r8(iz1,jz0);
219             dx20             = _fjsp_sub_v2r8(ix2,jx0);
220             dy20             = _fjsp_sub_v2r8(iy2,jy0);
221             dz20             = _fjsp_sub_v2r8(iz2,jz0);
222             dx30             = _fjsp_sub_v2r8(ix3,jx0);
223             dy30             = _fjsp_sub_v2r8(iy3,jy0);
224             dz30             = _fjsp_sub_v2r8(iz3,jz0);
225
226             /* Calculate squared distance and things based on it */
227             rsq00            = gmx_fjsp_calc_rsq_v2r8(dx00,dy00,dz00);
228             rsq10            = gmx_fjsp_calc_rsq_v2r8(dx10,dy10,dz10);
229             rsq20            = gmx_fjsp_calc_rsq_v2r8(dx20,dy20,dz20);
230             rsq30            = gmx_fjsp_calc_rsq_v2r8(dx30,dy30,dz30);
231
232             rinv00           = gmx_fjsp_invsqrt_v2r8(rsq00);
233             rinv10           = gmx_fjsp_invsqrt_v2r8(rsq10);
234             rinv20           = gmx_fjsp_invsqrt_v2r8(rsq20);
235             rinv30           = gmx_fjsp_invsqrt_v2r8(rsq30);
236
237             rinvsq00         = _fjsp_mul_v2r8(rinv00,rinv00);
238             rinvsq10         = _fjsp_mul_v2r8(rinv10,rinv10);
239             rinvsq20         = _fjsp_mul_v2r8(rinv20,rinv20);
240             rinvsq30         = _fjsp_mul_v2r8(rinv30,rinv30);
241
242             /* Load parameters for j particles */
243             jq0              = gmx_fjsp_load_2real_swizzle_v2r8(charge+jnrA+0,charge+jnrB+0);
244             vdwjidx0A        = 2*vdwtype[jnrA+0];
245             vdwjidx0B        = 2*vdwtype[jnrB+0];
246
247             fjx0             = _fjsp_setzero_v2r8();
248             fjy0             = _fjsp_setzero_v2r8();
249             fjz0             = _fjsp_setzero_v2r8();
250
251             /**************************
252              * CALCULATE INTERACTIONS *
253              **************************/
254
255             if (gmx_fjsp_any_lt_v2r8(rsq00,rcutoff2))
256             {
257
258             r00              = _fjsp_mul_v2r8(rsq00,rinv00);
259
260             /* Compute parameters for interactions between i and j atoms */
261             gmx_fjsp_load_2pair_swizzle_v2r8(vdwparam+vdwioffset0+vdwjidx0A,
262                                          vdwparam+vdwioffset0+vdwjidx0B,&c6_00,&c12_00);
263
264             /* LENNARD-JONES DISPERSION/REPULSION */
265
266             rinvsix          = _fjsp_mul_v2r8(_fjsp_mul_v2r8(rinvsq00,rinvsq00),rinvsq00);
267             vvdw6            = _fjsp_mul_v2r8(c6_00,rinvsix);
268             vvdw12           = _fjsp_mul_v2r8(c12_00,_fjsp_mul_v2r8(rinvsix,rinvsix));
269             vvdw             = _fjsp_msub_v2r8( vvdw12,one_twelfth, _fjsp_mul_v2r8(vvdw6,one_sixth) );
270             fvdw             = _fjsp_mul_v2r8(_fjsp_sub_v2r8(vvdw12,vvdw6),rinvsq00);
271
272             d                = _fjsp_sub_v2r8(r00,rswitch);
273             d                = _fjsp_max_v2r8(d,_fjsp_setzero_v2r8());
274             d2               = _fjsp_mul_v2r8(d,d);
275             sw               = _fjsp_add_v2r8(one,_fjsp_mul_v2r8(d2,_fjsp_mul_v2r8(d,_fjsp_madd_v2r8(d,_fjsp_madd_v2r8(d,swV5,swV4),swV3))));
276
277             dsw              = _fjsp_mul_v2r8(d2,_fjsp_madd_v2r8(d,_fjsp_madd_v2r8(d,swF4,swF3),swF2));
278
279             /* Evaluate switch function */
280             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
281             fvdw             = _fjsp_msub_v2r8( fvdw,sw , _fjsp_mul_v2r8(rinv00,_fjsp_mul_v2r8(vvdw,dsw)) );
282             vvdw             = _fjsp_mul_v2r8(vvdw,sw);
283             cutoff_mask      = _fjsp_cmplt_v2r8(rsq00,rcutoff2);
284
285             /* Update potential sum for this i atom from the interaction with this j atom. */
286             vvdw             = _fjsp_and_v2r8(vvdw,cutoff_mask);
287             vvdwsum          = _fjsp_add_v2r8(vvdwsum,vvdw);
288
289             fscal            = fvdw;
290
291             fscal            = _fjsp_and_v2r8(fscal,cutoff_mask);
292
293             /* Update vectorial force */
294             fix0             = _fjsp_madd_v2r8(dx00,fscal,fix0);
295             fiy0             = _fjsp_madd_v2r8(dy00,fscal,fiy0);
296             fiz0             = _fjsp_madd_v2r8(dz00,fscal,fiz0);
297             
298             fjx0             = _fjsp_madd_v2r8(dx00,fscal,fjx0);
299             fjy0             = _fjsp_madd_v2r8(dy00,fscal,fjy0);
300             fjz0             = _fjsp_madd_v2r8(dz00,fscal,fjz0);
301
302             }
303
304             /**************************
305              * CALCULATE INTERACTIONS *
306              **************************/
307
308             if (gmx_fjsp_any_lt_v2r8(rsq10,rcutoff2))
309             {
310
311             /* Compute parameters for interactions between i and j atoms */
312             qq10             = _fjsp_mul_v2r8(iq1,jq0);
313
314             /* REACTION-FIELD ELECTROSTATICS */
315             velec            = _fjsp_mul_v2r8(qq10,_fjsp_sub_v2r8(_fjsp_madd_v2r8(krf,rsq10,rinv10),crf));
316             felec            = _fjsp_mul_v2r8(qq10,_fjsp_msub_v2r8(rinv10,rinvsq10,krf2));
317
318             cutoff_mask      = _fjsp_cmplt_v2r8(rsq10,rcutoff2);
319
320             /* Update potential sum for this i atom from the interaction with this j atom. */
321             velec            = _fjsp_and_v2r8(velec,cutoff_mask);
322             velecsum         = _fjsp_add_v2r8(velecsum,velec);
323
324             fscal            = felec;
325
326             fscal            = _fjsp_and_v2r8(fscal,cutoff_mask);
327
328             /* Update vectorial force */
329             fix1             = _fjsp_madd_v2r8(dx10,fscal,fix1);
330             fiy1             = _fjsp_madd_v2r8(dy10,fscal,fiy1);
331             fiz1             = _fjsp_madd_v2r8(dz10,fscal,fiz1);
332             
333             fjx0             = _fjsp_madd_v2r8(dx10,fscal,fjx0);
334             fjy0             = _fjsp_madd_v2r8(dy10,fscal,fjy0);
335             fjz0             = _fjsp_madd_v2r8(dz10,fscal,fjz0);
336
337             }
338
339             /**************************
340              * CALCULATE INTERACTIONS *
341              **************************/
342
343             if (gmx_fjsp_any_lt_v2r8(rsq20,rcutoff2))
344             {
345
346             /* Compute parameters for interactions between i and j atoms */
347             qq20             = _fjsp_mul_v2r8(iq2,jq0);
348
349             /* REACTION-FIELD ELECTROSTATICS */
350             velec            = _fjsp_mul_v2r8(qq20,_fjsp_sub_v2r8(_fjsp_madd_v2r8(krf,rsq20,rinv20),crf));
351             felec            = _fjsp_mul_v2r8(qq20,_fjsp_msub_v2r8(rinv20,rinvsq20,krf2));
352
353             cutoff_mask      = _fjsp_cmplt_v2r8(rsq20,rcutoff2);
354
355             /* Update potential sum for this i atom from the interaction with this j atom. */
356             velec            = _fjsp_and_v2r8(velec,cutoff_mask);
357             velecsum         = _fjsp_add_v2r8(velecsum,velec);
358
359             fscal            = felec;
360
361             fscal            = _fjsp_and_v2r8(fscal,cutoff_mask);
362
363             /* Update vectorial force */
364             fix2             = _fjsp_madd_v2r8(dx20,fscal,fix2);
365             fiy2             = _fjsp_madd_v2r8(dy20,fscal,fiy2);
366             fiz2             = _fjsp_madd_v2r8(dz20,fscal,fiz2);
367             
368             fjx0             = _fjsp_madd_v2r8(dx20,fscal,fjx0);
369             fjy0             = _fjsp_madd_v2r8(dy20,fscal,fjy0);
370             fjz0             = _fjsp_madd_v2r8(dz20,fscal,fjz0);
371
372             }
373
374             /**************************
375              * CALCULATE INTERACTIONS *
376              **************************/
377
378             if (gmx_fjsp_any_lt_v2r8(rsq30,rcutoff2))
379             {
380
381             /* Compute parameters for interactions between i and j atoms */
382             qq30             = _fjsp_mul_v2r8(iq3,jq0);
383
384             /* REACTION-FIELD ELECTROSTATICS */
385             velec            = _fjsp_mul_v2r8(qq30,_fjsp_sub_v2r8(_fjsp_madd_v2r8(krf,rsq30,rinv30),crf));
386             felec            = _fjsp_mul_v2r8(qq30,_fjsp_msub_v2r8(rinv30,rinvsq30,krf2));
387
388             cutoff_mask      = _fjsp_cmplt_v2r8(rsq30,rcutoff2);
389
390             /* Update potential sum for this i atom from the interaction with this j atom. */
391             velec            = _fjsp_and_v2r8(velec,cutoff_mask);
392             velecsum         = _fjsp_add_v2r8(velecsum,velec);
393
394             fscal            = felec;
395
396             fscal            = _fjsp_and_v2r8(fscal,cutoff_mask);
397
398             /* Update vectorial force */
399             fix3             = _fjsp_madd_v2r8(dx30,fscal,fix3);
400             fiy3             = _fjsp_madd_v2r8(dy30,fscal,fiy3);
401             fiz3             = _fjsp_madd_v2r8(dz30,fscal,fiz3);
402             
403             fjx0             = _fjsp_madd_v2r8(dx30,fscal,fjx0);
404             fjy0             = _fjsp_madd_v2r8(dy30,fscal,fjy0);
405             fjz0             = _fjsp_madd_v2r8(dz30,fscal,fjz0);
406
407             }
408
409             gmx_fjsp_decrement_1rvec_2ptr_swizzle_v2r8(f+j_coord_offsetA,f+j_coord_offsetB,fjx0,fjy0,fjz0);
410
411             /* Inner loop uses 182 flops */
412         }
413
414         if(jidx<j_index_end)
415         {
416
417             jnrA             = jjnr[jidx];
418             j_coord_offsetA  = DIM*jnrA;
419
420             /* load j atom coordinates */
421             gmx_fjsp_load_1rvec_1ptr_swizzle_v2r8(x+j_coord_offsetA,
422                                               &jx0,&jy0,&jz0);
423
424             /* Calculate displacement vector */
425             dx00             = _fjsp_sub_v2r8(ix0,jx0);
426             dy00             = _fjsp_sub_v2r8(iy0,jy0);
427             dz00             = _fjsp_sub_v2r8(iz0,jz0);
428             dx10             = _fjsp_sub_v2r8(ix1,jx0);
429             dy10             = _fjsp_sub_v2r8(iy1,jy0);
430             dz10             = _fjsp_sub_v2r8(iz1,jz0);
431             dx20             = _fjsp_sub_v2r8(ix2,jx0);
432             dy20             = _fjsp_sub_v2r8(iy2,jy0);
433             dz20             = _fjsp_sub_v2r8(iz2,jz0);
434             dx30             = _fjsp_sub_v2r8(ix3,jx0);
435             dy30             = _fjsp_sub_v2r8(iy3,jy0);
436             dz30             = _fjsp_sub_v2r8(iz3,jz0);
437
438             /* Calculate squared distance and things based on it */
439             rsq00            = gmx_fjsp_calc_rsq_v2r8(dx00,dy00,dz00);
440             rsq10            = gmx_fjsp_calc_rsq_v2r8(dx10,dy10,dz10);
441             rsq20            = gmx_fjsp_calc_rsq_v2r8(dx20,dy20,dz20);
442             rsq30            = gmx_fjsp_calc_rsq_v2r8(dx30,dy30,dz30);
443
444             rinv00           = gmx_fjsp_invsqrt_v2r8(rsq00);
445             rinv10           = gmx_fjsp_invsqrt_v2r8(rsq10);
446             rinv20           = gmx_fjsp_invsqrt_v2r8(rsq20);
447             rinv30           = gmx_fjsp_invsqrt_v2r8(rsq30);
448
449             rinvsq00         = _fjsp_mul_v2r8(rinv00,rinv00);
450             rinvsq10         = _fjsp_mul_v2r8(rinv10,rinv10);
451             rinvsq20         = _fjsp_mul_v2r8(rinv20,rinv20);
452             rinvsq30         = _fjsp_mul_v2r8(rinv30,rinv30);
453
454             /* Load parameters for j particles */
455             jq0              = _fjsp_loadl_v2r8(_fjsp_setzero_v2r8(),charge+jnrA+0);
456             vdwjidx0A        = 2*vdwtype[jnrA+0];
457
458             fjx0             = _fjsp_setzero_v2r8();
459             fjy0             = _fjsp_setzero_v2r8();
460             fjz0             = _fjsp_setzero_v2r8();
461
462             /**************************
463              * CALCULATE INTERACTIONS *
464              **************************/
465
466             if (gmx_fjsp_any_lt_v2r8(rsq00,rcutoff2))
467             {
468
469             r00              = _fjsp_mul_v2r8(rsq00,rinv00);
470
471             /* Compute parameters for interactions between i and j atoms */
472             gmx_fjsp_load_1pair_swizzle_v2r8(vdwparam+vdwioffset0+vdwjidx0A,&c6_00,&c12_00);
473
474             /* LENNARD-JONES DISPERSION/REPULSION */
475
476             rinvsix          = _fjsp_mul_v2r8(_fjsp_mul_v2r8(rinvsq00,rinvsq00),rinvsq00);
477             vvdw6            = _fjsp_mul_v2r8(c6_00,rinvsix);
478             vvdw12           = _fjsp_mul_v2r8(c12_00,_fjsp_mul_v2r8(rinvsix,rinvsix));
479             vvdw             = _fjsp_msub_v2r8( vvdw12,one_twelfth, _fjsp_mul_v2r8(vvdw6,one_sixth) );
480             fvdw             = _fjsp_mul_v2r8(_fjsp_sub_v2r8(vvdw12,vvdw6),rinvsq00);
481
482             d                = _fjsp_sub_v2r8(r00,rswitch);
483             d                = _fjsp_max_v2r8(d,_fjsp_setzero_v2r8());
484             d2               = _fjsp_mul_v2r8(d,d);
485             sw               = _fjsp_add_v2r8(one,_fjsp_mul_v2r8(d2,_fjsp_mul_v2r8(d,_fjsp_madd_v2r8(d,_fjsp_madd_v2r8(d,swV5,swV4),swV3))));
486
487             dsw              = _fjsp_mul_v2r8(d2,_fjsp_madd_v2r8(d,_fjsp_madd_v2r8(d,swF4,swF3),swF2));
488
489             /* Evaluate switch function */
490             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
491             fvdw             = _fjsp_msub_v2r8( fvdw,sw , _fjsp_mul_v2r8(rinv00,_fjsp_mul_v2r8(vvdw,dsw)) );
492             vvdw             = _fjsp_mul_v2r8(vvdw,sw);
493             cutoff_mask      = _fjsp_cmplt_v2r8(rsq00,rcutoff2);
494
495             /* Update potential sum for this i atom from the interaction with this j atom. */
496             vvdw             = _fjsp_and_v2r8(vvdw,cutoff_mask);
497             vvdw             = _fjsp_unpacklo_v2r8(vvdw,_fjsp_setzero_v2r8());
498             vvdwsum          = _fjsp_add_v2r8(vvdwsum,vvdw);
499
500             fscal            = fvdw;
501
502             fscal            = _fjsp_and_v2r8(fscal,cutoff_mask);
503
504             fscal            = _fjsp_unpacklo_v2r8(fscal,_fjsp_setzero_v2r8());
505
506             /* Update vectorial force */
507             fix0             = _fjsp_madd_v2r8(dx00,fscal,fix0);
508             fiy0             = _fjsp_madd_v2r8(dy00,fscal,fiy0);
509             fiz0             = _fjsp_madd_v2r8(dz00,fscal,fiz0);
510             
511             fjx0             = _fjsp_madd_v2r8(dx00,fscal,fjx0);
512             fjy0             = _fjsp_madd_v2r8(dy00,fscal,fjy0);
513             fjz0             = _fjsp_madd_v2r8(dz00,fscal,fjz0);
514
515             }
516
517             /**************************
518              * CALCULATE INTERACTIONS *
519              **************************/
520
521             if (gmx_fjsp_any_lt_v2r8(rsq10,rcutoff2))
522             {
523
524             /* Compute parameters for interactions between i and j atoms */
525             qq10             = _fjsp_mul_v2r8(iq1,jq0);
526
527             /* REACTION-FIELD ELECTROSTATICS */
528             velec            = _fjsp_mul_v2r8(qq10,_fjsp_sub_v2r8(_fjsp_madd_v2r8(krf,rsq10,rinv10),crf));
529             felec            = _fjsp_mul_v2r8(qq10,_fjsp_msub_v2r8(rinv10,rinvsq10,krf2));
530
531             cutoff_mask      = _fjsp_cmplt_v2r8(rsq10,rcutoff2);
532
533             /* Update potential sum for this i atom from the interaction with this j atom. */
534             velec            = _fjsp_and_v2r8(velec,cutoff_mask);
535             velec            = _fjsp_unpacklo_v2r8(velec,_fjsp_setzero_v2r8());
536             velecsum         = _fjsp_add_v2r8(velecsum,velec);
537
538             fscal            = felec;
539
540             fscal            = _fjsp_and_v2r8(fscal,cutoff_mask);
541
542             fscal            = _fjsp_unpacklo_v2r8(fscal,_fjsp_setzero_v2r8());
543
544             /* Update vectorial force */
545             fix1             = _fjsp_madd_v2r8(dx10,fscal,fix1);
546             fiy1             = _fjsp_madd_v2r8(dy10,fscal,fiy1);
547             fiz1             = _fjsp_madd_v2r8(dz10,fscal,fiz1);
548             
549             fjx0             = _fjsp_madd_v2r8(dx10,fscal,fjx0);
550             fjy0             = _fjsp_madd_v2r8(dy10,fscal,fjy0);
551             fjz0             = _fjsp_madd_v2r8(dz10,fscal,fjz0);
552
553             }
554
555             /**************************
556              * CALCULATE INTERACTIONS *
557              **************************/
558
559             if (gmx_fjsp_any_lt_v2r8(rsq20,rcutoff2))
560             {
561
562             /* Compute parameters for interactions between i and j atoms */
563             qq20             = _fjsp_mul_v2r8(iq2,jq0);
564
565             /* REACTION-FIELD ELECTROSTATICS */
566             velec            = _fjsp_mul_v2r8(qq20,_fjsp_sub_v2r8(_fjsp_madd_v2r8(krf,rsq20,rinv20),crf));
567             felec            = _fjsp_mul_v2r8(qq20,_fjsp_msub_v2r8(rinv20,rinvsq20,krf2));
568
569             cutoff_mask      = _fjsp_cmplt_v2r8(rsq20,rcutoff2);
570
571             /* Update potential sum for this i atom from the interaction with this j atom. */
572             velec            = _fjsp_and_v2r8(velec,cutoff_mask);
573             velec            = _fjsp_unpacklo_v2r8(velec,_fjsp_setzero_v2r8());
574             velecsum         = _fjsp_add_v2r8(velecsum,velec);
575
576             fscal            = felec;
577
578             fscal            = _fjsp_and_v2r8(fscal,cutoff_mask);
579
580             fscal            = _fjsp_unpacklo_v2r8(fscal,_fjsp_setzero_v2r8());
581
582             /* Update vectorial force */
583             fix2             = _fjsp_madd_v2r8(dx20,fscal,fix2);
584             fiy2             = _fjsp_madd_v2r8(dy20,fscal,fiy2);
585             fiz2             = _fjsp_madd_v2r8(dz20,fscal,fiz2);
586             
587             fjx0             = _fjsp_madd_v2r8(dx20,fscal,fjx0);
588             fjy0             = _fjsp_madd_v2r8(dy20,fscal,fjy0);
589             fjz0             = _fjsp_madd_v2r8(dz20,fscal,fjz0);
590
591             }
592
593             /**************************
594              * CALCULATE INTERACTIONS *
595              **************************/
596
597             if (gmx_fjsp_any_lt_v2r8(rsq30,rcutoff2))
598             {
599
600             /* Compute parameters for interactions between i and j atoms */
601             qq30             = _fjsp_mul_v2r8(iq3,jq0);
602
603             /* REACTION-FIELD ELECTROSTATICS */
604             velec            = _fjsp_mul_v2r8(qq30,_fjsp_sub_v2r8(_fjsp_madd_v2r8(krf,rsq30,rinv30),crf));
605             felec            = _fjsp_mul_v2r8(qq30,_fjsp_msub_v2r8(rinv30,rinvsq30,krf2));
606
607             cutoff_mask      = _fjsp_cmplt_v2r8(rsq30,rcutoff2);
608
609             /* Update potential sum for this i atom from the interaction with this j atom. */
610             velec            = _fjsp_and_v2r8(velec,cutoff_mask);
611             velec            = _fjsp_unpacklo_v2r8(velec,_fjsp_setzero_v2r8());
612             velecsum         = _fjsp_add_v2r8(velecsum,velec);
613
614             fscal            = felec;
615
616             fscal            = _fjsp_and_v2r8(fscal,cutoff_mask);
617
618             fscal            = _fjsp_unpacklo_v2r8(fscal,_fjsp_setzero_v2r8());
619
620             /* Update vectorial force */
621             fix3             = _fjsp_madd_v2r8(dx30,fscal,fix3);
622             fiy3             = _fjsp_madd_v2r8(dy30,fscal,fiy3);
623             fiz3             = _fjsp_madd_v2r8(dz30,fscal,fiz3);
624             
625             fjx0             = _fjsp_madd_v2r8(dx30,fscal,fjx0);
626             fjy0             = _fjsp_madd_v2r8(dy30,fscal,fjy0);
627             fjz0             = _fjsp_madd_v2r8(dz30,fscal,fjz0);
628
629             }
630
631             gmx_fjsp_decrement_1rvec_1ptr_swizzle_v2r8(f+j_coord_offsetA,fjx0,fjy0,fjz0);
632
633             /* Inner loop uses 182 flops */
634         }
635
636         /* End of innermost loop */
637
638         gmx_fjsp_update_iforce_4atom_swizzle_v2r8(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
639                                               f+i_coord_offset,fshift+i_shift_offset);
640
641         ggid                        = gid[iidx];
642         /* Update potential energies */
643         gmx_fjsp_update_1pot_v2r8(velecsum,kernel_data->energygrp_elec+ggid);
644         gmx_fjsp_update_1pot_v2r8(vvdwsum,kernel_data->energygrp_vdw+ggid);
645
646         /* Increment number of inner iterations */
647         inneriter                  += j_index_end - j_index_start;
648
649         /* Outer loop uses 26 flops */
650     }
651
652     /* Increment number of outer iterations */
653     outeriter        += nri;
654
655     /* Update outer/inner flops */
656
657     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4_VF,outeriter*26 + inneriter*182);
658 }
659 /*
660  * Gromacs nonbonded kernel:   nb_kernel_ElecRFCut_VdwLJSw_GeomW4P1_F_sparc64_hpc_ace_double
661  * Electrostatics interaction: ReactionField
662  * VdW interaction:            LennardJones
663  * Geometry:                   Water4-Particle
664  * Calculate force/pot:        Force
665  */
666 void
667 nb_kernel_ElecRFCut_VdwLJSw_GeomW4P1_F_sparc64_hpc_ace_double
668                     (t_nblist * gmx_restrict                nlist,
669                      rvec * gmx_restrict                    xx,
670                      rvec * gmx_restrict                    ff,
671                      t_forcerec * gmx_restrict              fr,
672                      t_mdatoms * gmx_restrict               mdatoms,
673                      nb_kernel_data_t * gmx_restrict        kernel_data,
674                      t_nrnb * gmx_restrict                  nrnb)
675 {
676     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
677      * just 0 for non-waters.
678      * Suffixes A,B refer to j loop unrolling done with double precision SIMD, e.g. for the two different
679      * jnr indices corresponding to data put in the four positions in the SIMD register.
680      */
681     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
682     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
683     int              jnrA,jnrB;
684     int              j_coord_offsetA,j_coord_offsetB;
685     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
686     real             rcutoff_scalar;
687     real             *shiftvec,*fshift,*x,*f;
688     _fjsp_v2r8       tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
689     int              vdwioffset0;
690     _fjsp_v2r8       ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
691     int              vdwioffset1;
692     _fjsp_v2r8       ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
693     int              vdwioffset2;
694     _fjsp_v2r8       ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
695     int              vdwioffset3;
696     _fjsp_v2r8       ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
697     int              vdwjidx0A,vdwjidx0B;
698     _fjsp_v2r8       jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
699     _fjsp_v2r8       dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
700     _fjsp_v2r8       dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
701     _fjsp_v2r8       dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
702     _fjsp_v2r8       dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
703     _fjsp_v2r8       velec,felec,velecsum,facel,crf,krf,krf2;
704     real             *charge;
705     int              nvdwtype;
706     _fjsp_v2r8       rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
707     int              *vdwtype;
708     real             *vdwparam;
709     _fjsp_v2r8       one_sixth   = gmx_fjsp_set1_v2r8(1.0/6.0);
710     _fjsp_v2r8       one_twelfth = gmx_fjsp_set1_v2r8(1.0/12.0);
711     _fjsp_v2r8       rswitch,swV3,swV4,swV5,swF2,swF3,swF4,d,d2,sw,dsw;
712     real             rswitch_scalar,d_scalar;
713     _fjsp_v2r8       itab_tmp;
714     _fjsp_v2r8       dummy_mask,cutoff_mask;
715     _fjsp_v2r8       one     = gmx_fjsp_set1_v2r8(1.0);
716     _fjsp_v2r8       two     = gmx_fjsp_set1_v2r8(2.0);
717     union { _fjsp_v2r8 simd; long long int i[2]; } vfconv,gbconv,ewconv;
718
719     x                = xx[0];
720     f                = ff[0];
721
722     nri              = nlist->nri;
723     iinr             = nlist->iinr;
724     jindex           = nlist->jindex;
725     jjnr             = nlist->jjnr;
726     shiftidx         = nlist->shift;
727     gid              = nlist->gid;
728     shiftvec         = fr->shift_vec[0];
729     fshift           = fr->fshift[0];
730     facel            = gmx_fjsp_set1_v2r8(fr->epsfac);
731     charge           = mdatoms->chargeA;
732     krf              = gmx_fjsp_set1_v2r8(fr->ic->k_rf);
733     krf2             = gmx_fjsp_set1_v2r8(fr->ic->k_rf*2.0);
734     crf              = gmx_fjsp_set1_v2r8(fr->ic->c_rf);
735     nvdwtype         = fr->ntype;
736     vdwparam         = fr->nbfp;
737     vdwtype          = mdatoms->typeA;
738
739     /* Setup water-specific parameters */
740     inr              = nlist->iinr[0];
741     iq1              = _fjsp_mul_v2r8(facel,gmx_fjsp_set1_v2r8(charge[inr+1]));
742     iq2              = _fjsp_mul_v2r8(facel,gmx_fjsp_set1_v2r8(charge[inr+2]));
743     iq3              = _fjsp_mul_v2r8(facel,gmx_fjsp_set1_v2r8(charge[inr+3]));
744     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
745
746     /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
747     rcutoff_scalar   = fr->rcoulomb;
748     rcutoff          = gmx_fjsp_set1_v2r8(rcutoff_scalar);
749     rcutoff2         = _fjsp_mul_v2r8(rcutoff,rcutoff);
750
751     rswitch_scalar   = fr->rvdw_switch;
752     rswitch          = gmx_fjsp_set1_v2r8(rswitch_scalar);
753     /* Setup switch parameters */
754     d_scalar         = rcutoff_scalar-rswitch_scalar;
755     d                = gmx_fjsp_set1_v2r8(d_scalar);
756     swV3             = gmx_fjsp_set1_v2r8(-10.0/(d_scalar*d_scalar*d_scalar));
757     swV4             = gmx_fjsp_set1_v2r8( 15.0/(d_scalar*d_scalar*d_scalar*d_scalar));
758     swV5             = gmx_fjsp_set1_v2r8( -6.0/(d_scalar*d_scalar*d_scalar*d_scalar*d_scalar));
759     swF2             = gmx_fjsp_set1_v2r8(-30.0/(d_scalar*d_scalar*d_scalar));
760     swF3             = gmx_fjsp_set1_v2r8( 60.0/(d_scalar*d_scalar*d_scalar*d_scalar));
761     swF4             = gmx_fjsp_set1_v2r8(-30.0/(d_scalar*d_scalar*d_scalar*d_scalar*d_scalar));
762
763     /* Avoid stupid compiler warnings */
764     jnrA = jnrB = 0;
765     j_coord_offsetA = 0;
766     j_coord_offsetB = 0;
767
768     outeriter        = 0;
769     inneriter        = 0;
770
771     /* Start outer loop over neighborlists */
772     for(iidx=0; iidx<nri; iidx++)
773     {
774         /* Load shift vector for this list */
775         i_shift_offset   = DIM*shiftidx[iidx];
776
777         /* Load limits for loop over neighbors */
778         j_index_start    = jindex[iidx];
779         j_index_end      = jindex[iidx+1];
780
781         /* Get outer coordinate index */
782         inr              = iinr[iidx];
783         i_coord_offset   = DIM*inr;
784
785         /* Load i particle coords and add shift vector */
786         gmx_fjsp_load_shift_and_4rvec_broadcast_v2r8(shiftvec+i_shift_offset,x+i_coord_offset,
787                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
788
789         fix0             = _fjsp_setzero_v2r8();
790         fiy0             = _fjsp_setzero_v2r8();
791         fiz0             = _fjsp_setzero_v2r8();
792         fix1             = _fjsp_setzero_v2r8();
793         fiy1             = _fjsp_setzero_v2r8();
794         fiz1             = _fjsp_setzero_v2r8();
795         fix2             = _fjsp_setzero_v2r8();
796         fiy2             = _fjsp_setzero_v2r8();
797         fiz2             = _fjsp_setzero_v2r8();
798         fix3             = _fjsp_setzero_v2r8();
799         fiy3             = _fjsp_setzero_v2r8();
800         fiz3             = _fjsp_setzero_v2r8();
801
802         /* Start inner kernel loop */
803         for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
804         {
805
806             /* Get j neighbor index, and coordinate index */
807             jnrA             = jjnr[jidx];
808             jnrB             = jjnr[jidx+1];
809             j_coord_offsetA  = DIM*jnrA;
810             j_coord_offsetB  = DIM*jnrB;
811
812             /* load j atom coordinates */
813             gmx_fjsp_load_1rvec_2ptr_swizzle_v2r8(x+j_coord_offsetA,x+j_coord_offsetB,
814                                               &jx0,&jy0,&jz0);
815
816             /* Calculate displacement vector */
817             dx00             = _fjsp_sub_v2r8(ix0,jx0);
818             dy00             = _fjsp_sub_v2r8(iy0,jy0);
819             dz00             = _fjsp_sub_v2r8(iz0,jz0);
820             dx10             = _fjsp_sub_v2r8(ix1,jx0);
821             dy10             = _fjsp_sub_v2r8(iy1,jy0);
822             dz10             = _fjsp_sub_v2r8(iz1,jz0);
823             dx20             = _fjsp_sub_v2r8(ix2,jx0);
824             dy20             = _fjsp_sub_v2r8(iy2,jy0);
825             dz20             = _fjsp_sub_v2r8(iz2,jz0);
826             dx30             = _fjsp_sub_v2r8(ix3,jx0);
827             dy30             = _fjsp_sub_v2r8(iy3,jy0);
828             dz30             = _fjsp_sub_v2r8(iz3,jz0);
829
830             /* Calculate squared distance and things based on it */
831             rsq00            = gmx_fjsp_calc_rsq_v2r8(dx00,dy00,dz00);
832             rsq10            = gmx_fjsp_calc_rsq_v2r8(dx10,dy10,dz10);
833             rsq20            = gmx_fjsp_calc_rsq_v2r8(dx20,dy20,dz20);
834             rsq30            = gmx_fjsp_calc_rsq_v2r8(dx30,dy30,dz30);
835
836             rinv00           = gmx_fjsp_invsqrt_v2r8(rsq00);
837             rinv10           = gmx_fjsp_invsqrt_v2r8(rsq10);
838             rinv20           = gmx_fjsp_invsqrt_v2r8(rsq20);
839             rinv30           = gmx_fjsp_invsqrt_v2r8(rsq30);
840
841             rinvsq00         = _fjsp_mul_v2r8(rinv00,rinv00);
842             rinvsq10         = _fjsp_mul_v2r8(rinv10,rinv10);
843             rinvsq20         = _fjsp_mul_v2r8(rinv20,rinv20);
844             rinvsq30         = _fjsp_mul_v2r8(rinv30,rinv30);
845
846             /* Load parameters for j particles */
847             jq0              = gmx_fjsp_load_2real_swizzle_v2r8(charge+jnrA+0,charge+jnrB+0);
848             vdwjidx0A        = 2*vdwtype[jnrA+0];
849             vdwjidx0B        = 2*vdwtype[jnrB+0];
850
851             fjx0             = _fjsp_setzero_v2r8();
852             fjy0             = _fjsp_setzero_v2r8();
853             fjz0             = _fjsp_setzero_v2r8();
854
855             /**************************
856              * CALCULATE INTERACTIONS *
857              **************************/
858
859             if (gmx_fjsp_any_lt_v2r8(rsq00,rcutoff2))
860             {
861
862             r00              = _fjsp_mul_v2r8(rsq00,rinv00);
863
864             /* Compute parameters for interactions between i and j atoms */
865             gmx_fjsp_load_2pair_swizzle_v2r8(vdwparam+vdwioffset0+vdwjidx0A,
866                                          vdwparam+vdwioffset0+vdwjidx0B,&c6_00,&c12_00);
867
868             /* LENNARD-JONES DISPERSION/REPULSION */
869
870             rinvsix          = _fjsp_mul_v2r8(_fjsp_mul_v2r8(rinvsq00,rinvsq00),rinvsq00);
871             vvdw6            = _fjsp_mul_v2r8(c6_00,rinvsix);
872             vvdw12           = _fjsp_mul_v2r8(c12_00,_fjsp_mul_v2r8(rinvsix,rinvsix));
873             vvdw             = _fjsp_msub_v2r8( vvdw12,one_twelfth, _fjsp_mul_v2r8(vvdw6,one_sixth) );
874             fvdw             = _fjsp_mul_v2r8(_fjsp_sub_v2r8(vvdw12,vvdw6),rinvsq00);
875
876             d                = _fjsp_sub_v2r8(r00,rswitch);
877             d                = _fjsp_max_v2r8(d,_fjsp_setzero_v2r8());
878             d2               = _fjsp_mul_v2r8(d,d);
879             sw               = _fjsp_add_v2r8(one,_fjsp_mul_v2r8(d2,_fjsp_mul_v2r8(d,_fjsp_madd_v2r8(d,_fjsp_madd_v2r8(d,swV5,swV4),swV3))));
880
881             dsw              = _fjsp_mul_v2r8(d2,_fjsp_madd_v2r8(d,_fjsp_madd_v2r8(d,swF4,swF3),swF2));
882
883             /* Evaluate switch function */
884             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
885             fvdw             = _fjsp_msub_v2r8( fvdw,sw , _fjsp_mul_v2r8(rinv00,_fjsp_mul_v2r8(vvdw,dsw)) );
886             cutoff_mask      = _fjsp_cmplt_v2r8(rsq00,rcutoff2);
887
888             fscal            = fvdw;
889
890             fscal            = _fjsp_and_v2r8(fscal,cutoff_mask);
891
892             /* Update vectorial force */
893             fix0             = _fjsp_madd_v2r8(dx00,fscal,fix0);
894             fiy0             = _fjsp_madd_v2r8(dy00,fscal,fiy0);
895             fiz0             = _fjsp_madd_v2r8(dz00,fscal,fiz0);
896             
897             fjx0             = _fjsp_madd_v2r8(dx00,fscal,fjx0);
898             fjy0             = _fjsp_madd_v2r8(dy00,fscal,fjy0);
899             fjz0             = _fjsp_madd_v2r8(dz00,fscal,fjz0);
900
901             }
902
903             /**************************
904              * CALCULATE INTERACTIONS *
905              **************************/
906
907             if (gmx_fjsp_any_lt_v2r8(rsq10,rcutoff2))
908             {
909
910             /* Compute parameters for interactions between i and j atoms */
911             qq10             = _fjsp_mul_v2r8(iq1,jq0);
912
913             /* REACTION-FIELD ELECTROSTATICS */
914             felec            = _fjsp_mul_v2r8(qq10,_fjsp_msub_v2r8(rinv10,rinvsq10,krf2));
915
916             cutoff_mask      = _fjsp_cmplt_v2r8(rsq10,rcutoff2);
917
918             fscal            = felec;
919
920             fscal            = _fjsp_and_v2r8(fscal,cutoff_mask);
921
922             /* Update vectorial force */
923             fix1             = _fjsp_madd_v2r8(dx10,fscal,fix1);
924             fiy1             = _fjsp_madd_v2r8(dy10,fscal,fiy1);
925             fiz1             = _fjsp_madd_v2r8(dz10,fscal,fiz1);
926             
927             fjx0             = _fjsp_madd_v2r8(dx10,fscal,fjx0);
928             fjy0             = _fjsp_madd_v2r8(dy10,fscal,fjy0);
929             fjz0             = _fjsp_madd_v2r8(dz10,fscal,fjz0);
930
931             }
932
933             /**************************
934              * CALCULATE INTERACTIONS *
935              **************************/
936
937             if (gmx_fjsp_any_lt_v2r8(rsq20,rcutoff2))
938             {
939
940             /* Compute parameters for interactions between i and j atoms */
941             qq20             = _fjsp_mul_v2r8(iq2,jq0);
942
943             /* REACTION-FIELD ELECTROSTATICS */
944             felec            = _fjsp_mul_v2r8(qq20,_fjsp_msub_v2r8(rinv20,rinvsq20,krf2));
945
946             cutoff_mask      = _fjsp_cmplt_v2r8(rsq20,rcutoff2);
947
948             fscal            = felec;
949
950             fscal            = _fjsp_and_v2r8(fscal,cutoff_mask);
951
952             /* Update vectorial force */
953             fix2             = _fjsp_madd_v2r8(dx20,fscal,fix2);
954             fiy2             = _fjsp_madd_v2r8(dy20,fscal,fiy2);
955             fiz2             = _fjsp_madd_v2r8(dz20,fscal,fiz2);
956             
957             fjx0             = _fjsp_madd_v2r8(dx20,fscal,fjx0);
958             fjy0             = _fjsp_madd_v2r8(dy20,fscal,fjy0);
959             fjz0             = _fjsp_madd_v2r8(dz20,fscal,fjz0);
960
961             }
962
963             /**************************
964              * CALCULATE INTERACTIONS *
965              **************************/
966
967             if (gmx_fjsp_any_lt_v2r8(rsq30,rcutoff2))
968             {
969
970             /* Compute parameters for interactions between i and j atoms */
971             qq30             = _fjsp_mul_v2r8(iq3,jq0);
972
973             /* REACTION-FIELD ELECTROSTATICS */
974             felec            = _fjsp_mul_v2r8(qq30,_fjsp_msub_v2r8(rinv30,rinvsq30,krf2));
975
976             cutoff_mask      = _fjsp_cmplt_v2r8(rsq30,rcutoff2);
977
978             fscal            = felec;
979
980             fscal            = _fjsp_and_v2r8(fscal,cutoff_mask);
981
982             /* Update vectorial force */
983             fix3             = _fjsp_madd_v2r8(dx30,fscal,fix3);
984             fiy3             = _fjsp_madd_v2r8(dy30,fscal,fiy3);
985             fiz3             = _fjsp_madd_v2r8(dz30,fscal,fiz3);
986             
987             fjx0             = _fjsp_madd_v2r8(dx30,fscal,fjx0);
988             fjy0             = _fjsp_madd_v2r8(dy30,fscal,fjy0);
989             fjz0             = _fjsp_madd_v2r8(dz30,fscal,fjz0);
990
991             }
992
993             gmx_fjsp_decrement_1rvec_2ptr_swizzle_v2r8(f+j_coord_offsetA,f+j_coord_offsetB,fjx0,fjy0,fjz0);
994
995             /* Inner loop uses 161 flops */
996         }
997
998         if(jidx<j_index_end)
999         {
1000
1001             jnrA             = jjnr[jidx];
1002             j_coord_offsetA  = DIM*jnrA;
1003
1004             /* load j atom coordinates */
1005             gmx_fjsp_load_1rvec_1ptr_swizzle_v2r8(x+j_coord_offsetA,
1006                                               &jx0,&jy0,&jz0);
1007
1008             /* Calculate displacement vector */
1009             dx00             = _fjsp_sub_v2r8(ix0,jx0);
1010             dy00             = _fjsp_sub_v2r8(iy0,jy0);
1011             dz00             = _fjsp_sub_v2r8(iz0,jz0);
1012             dx10             = _fjsp_sub_v2r8(ix1,jx0);
1013             dy10             = _fjsp_sub_v2r8(iy1,jy0);
1014             dz10             = _fjsp_sub_v2r8(iz1,jz0);
1015             dx20             = _fjsp_sub_v2r8(ix2,jx0);
1016             dy20             = _fjsp_sub_v2r8(iy2,jy0);
1017             dz20             = _fjsp_sub_v2r8(iz2,jz0);
1018             dx30             = _fjsp_sub_v2r8(ix3,jx0);
1019             dy30             = _fjsp_sub_v2r8(iy3,jy0);
1020             dz30             = _fjsp_sub_v2r8(iz3,jz0);
1021
1022             /* Calculate squared distance and things based on it */
1023             rsq00            = gmx_fjsp_calc_rsq_v2r8(dx00,dy00,dz00);
1024             rsq10            = gmx_fjsp_calc_rsq_v2r8(dx10,dy10,dz10);
1025             rsq20            = gmx_fjsp_calc_rsq_v2r8(dx20,dy20,dz20);
1026             rsq30            = gmx_fjsp_calc_rsq_v2r8(dx30,dy30,dz30);
1027
1028             rinv00           = gmx_fjsp_invsqrt_v2r8(rsq00);
1029             rinv10           = gmx_fjsp_invsqrt_v2r8(rsq10);
1030             rinv20           = gmx_fjsp_invsqrt_v2r8(rsq20);
1031             rinv30           = gmx_fjsp_invsqrt_v2r8(rsq30);
1032
1033             rinvsq00         = _fjsp_mul_v2r8(rinv00,rinv00);
1034             rinvsq10         = _fjsp_mul_v2r8(rinv10,rinv10);
1035             rinvsq20         = _fjsp_mul_v2r8(rinv20,rinv20);
1036             rinvsq30         = _fjsp_mul_v2r8(rinv30,rinv30);
1037
1038             /* Load parameters for j particles */
1039             jq0              = _fjsp_loadl_v2r8(_fjsp_setzero_v2r8(),charge+jnrA+0);
1040             vdwjidx0A        = 2*vdwtype[jnrA+0];
1041
1042             fjx0             = _fjsp_setzero_v2r8();
1043             fjy0             = _fjsp_setzero_v2r8();
1044             fjz0             = _fjsp_setzero_v2r8();
1045
1046             /**************************
1047              * CALCULATE INTERACTIONS *
1048              **************************/
1049
1050             if (gmx_fjsp_any_lt_v2r8(rsq00,rcutoff2))
1051             {
1052
1053             r00              = _fjsp_mul_v2r8(rsq00,rinv00);
1054
1055             /* Compute parameters for interactions between i and j atoms */
1056             gmx_fjsp_load_1pair_swizzle_v2r8(vdwparam+vdwioffset0+vdwjidx0A,&c6_00,&c12_00);
1057
1058             /* LENNARD-JONES DISPERSION/REPULSION */
1059
1060             rinvsix          = _fjsp_mul_v2r8(_fjsp_mul_v2r8(rinvsq00,rinvsq00),rinvsq00);
1061             vvdw6            = _fjsp_mul_v2r8(c6_00,rinvsix);
1062             vvdw12           = _fjsp_mul_v2r8(c12_00,_fjsp_mul_v2r8(rinvsix,rinvsix));
1063             vvdw             = _fjsp_msub_v2r8( vvdw12,one_twelfth, _fjsp_mul_v2r8(vvdw6,one_sixth) );
1064             fvdw             = _fjsp_mul_v2r8(_fjsp_sub_v2r8(vvdw12,vvdw6),rinvsq00);
1065
1066             d                = _fjsp_sub_v2r8(r00,rswitch);
1067             d                = _fjsp_max_v2r8(d,_fjsp_setzero_v2r8());
1068             d2               = _fjsp_mul_v2r8(d,d);
1069             sw               = _fjsp_add_v2r8(one,_fjsp_mul_v2r8(d2,_fjsp_mul_v2r8(d,_fjsp_madd_v2r8(d,_fjsp_madd_v2r8(d,swV5,swV4),swV3))));
1070
1071             dsw              = _fjsp_mul_v2r8(d2,_fjsp_madd_v2r8(d,_fjsp_madd_v2r8(d,swF4,swF3),swF2));
1072
1073             /* Evaluate switch function */
1074             /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1075             fvdw             = _fjsp_msub_v2r8( fvdw,sw , _fjsp_mul_v2r8(rinv00,_fjsp_mul_v2r8(vvdw,dsw)) );
1076             cutoff_mask      = _fjsp_cmplt_v2r8(rsq00,rcutoff2);
1077
1078             fscal            = fvdw;
1079
1080             fscal            = _fjsp_and_v2r8(fscal,cutoff_mask);
1081
1082             fscal            = _fjsp_unpacklo_v2r8(fscal,_fjsp_setzero_v2r8());
1083
1084             /* Update vectorial force */
1085             fix0             = _fjsp_madd_v2r8(dx00,fscal,fix0);
1086             fiy0             = _fjsp_madd_v2r8(dy00,fscal,fiy0);
1087             fiz0             = _fjsp_madd_v2r8(dz00,fscal,fiz0);
1088             
1089             fjx0             = _fjsp_madd_v2r8(dx00,fscal,fjx0);
1090             fjy0             = _fjsp_madd_v2r8(dy00,fscal,fjy0);
1091             fjz0             = _fjsp_madd_v2r8(dz00,fscal,fjz0);
1092
1093             }
1094
1095             /**************************
1096              * CALCULATE INTERACTIONS *
1097              **************************/
1098
1099             if (gmx_fjsp_any_lt_v2r8(rsq10,rcutoff2))
1100             {
1101
1102             /* Compute parameters for interactions between i and j atoms */
1103             qq10             = _fjsp_mul_v2r8(iq1,jq0);
1104
1105             /* REACTION-FIELD ELECTROSTATICS */
1106             felec            = _fjsp_mul_v2r8(qq10,_fjsp_msub_v2r8(rinv10,rinvsq10,krf2));
1107
1108             cutoff_mask      = _fjsp_cmplt_v2r8(rsq10,rcutoff2);
1109
1110             fscal            = felec;
1111
1112             fscal            = _fjsp_and_v2r8(fscal,cutoff_mask);
1113
1114             fscal            = _fjsp_unpacklo_v2r8(fscal,_fjsp_setzero_v2r8());
1115
1116             /* Update vectorial force */
1117             fix1             = _fjsp_madd_v2r8(dx10,fscal,fix1);
1118             fiy1             = _fjsp_madd_v2r8(dy10,fscal,fiy1);
1119             fiz1             = _fjsp_madd_v2r8(dz10,fscal,fiz1);
1120             
1121             fjx0             = _fjsp_madd_v2r8(dx10,fscal,fjx0);
1122             fjy0             = _fjsp_madd_v2r8(dy10,fscal,fjy0);
1123             fjz0             = _fjsp_madd_v2r8(dz10,fscal,fjz0);
1124
1125             }
1126
1127             /**************************
1128              * CALCULATE INTERACTIONS *
1129              **************************/
1130
1131             if (gmx_fjsp_any_lt_v2r8(rsq20,rcutoff2))
1132             {
1133
1134             /* Compute parameters for interactions between i and j atoms */
1135             qq20             = _fjsp_mul_v2r8(iq2,jq0);
1136
1137             /* REACTION-FIELD ELECTROSTATICS */
1138             felec            = _fjsp_mul_v2r8(qq20,_fjsp_msub_v2r8(rinv20,rinvsq20,krf2));
1139
1140             cutoff_mask      = _fjsp_cmplt_v2r8(rsq20,rcutoff2);
1141
1142             fscal            = felec;
1143
1144             fscal            = _fjsp_and_v2r8(fscal,cutoff_mask);
1145
1146             fscal            = _fjsp_unpacklo_v2r8(fscal,_fjsp_setzero_v2r8());
1147
1148             /* Update vectorial force */
1149             fix2             = _fjsp_madd_v2r8(dx20,fscal,fix2);
1150             fiy2             = _fjsp_madd_v2r8(dy20,fscal,fiy2);
1151             fiz2             = _fjsp_madd_v2r8(dz20,fscal,fiz2);
1152             
1153             fjx0             = _fjsp_madd_v2r8(dx20,fscal,fjx0);
1154             fjy0             = _fjsp_madd_v2r8(dy20,fscal,fjy0);
1155             fjz0             = _fjsp_madd_v2r8(dz20,fscal,fjz0);
1156
1157             }
1158
1159             /**************************
1160              * CALCULATE INTERACTIONS *
1161              **************************/
1162
1163             if (gmx_fjsp_any_lt_v2r8(rsq30,rcutoff2))
1164             {
1165
1166             /* Compute parameters for interactions between i and j atoms */
1167             qq30             = _fjsp_mul_v2r8(iq3,jq0);
1168
1169             /* REACTION-FIELD ELECTROSTATICS */
1170             felec            = _fjsp_mul_v2r8(qq30,_fjsp_msub_v2r8(rinv30,rinvsq30,krf2));
1171
1172             cutoff_mask      = _fjsp_cmplt_v2r8(rsq30,rcutoff2);
1173
1174             fscal            = felec;
1175
1176             fscal            = _fjsp_and_v2r8(fscal,cutoff_mask);
1177
1178             fscal            = _fjsp_unpacklo_v2r8(fscal,_fjsp_setzero_v2r8());
1179
1180             /* Update vectorial force */
1181             fix3             = _fjsp_madd_v2r8(dx30,fscal,fix3);
1182             fiy3             = _fjsp_madd_v2r8(dy30,fscal,fiy3);
1183             fiz3             = _fjsp_madd_v2r8(dz30,fscal,fiz3);
1184             
1185             fjx0             = _fjsp_madd_v2r8(dx30,fscal,fjx0);
1186             fjy0             = _fjsp_madd_v2r8(dy30,fscal,fjy0);
1187             fjz0             = _fjsp_madd_v2r8(dz30,fscal,fjz0);
1188
1189             }
1190
1191             gmx_fjsp_decrement_1rvec_1ptr_swizzle_v2r8(f+j_coord_offsetA,fjx0,fjy0,fjz0);
1192
1193             /* Inner loop uses 161 flops */
1194         }
1195
1196         /* End of innermost loop */
1197
1198         gmx_fjsp_update_iforce_4atom_swizzle_v2r8(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
1199                                               f+i_coord_offset,fshift+i_shift_offset);
1200
1201         /* Increment number of inner iterations */
1202         inneriter                  += j_index_end - j_index_start;
1203
1204         /* Outer loop uses 24 flops */
1205     }
1206
1207     /* Increment number of outer iterations */
1208     outeriter        += nri;
1209
1210     /* Update outer/inner flops */
1211
1212     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4_F,outeriter*24 + inneriter*161);
1213 }