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