32bb365c287a2e1307aeb8b4bb44eb2e885ceeec
[alexxy/gromacs.git] / src / gmxlib / nonbonded / nb_kernel_avx_128_fma_double / nb_kernel_ElecRF_VdwLJ_GeomW3P1_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_GeomW3P1_VF_avx_128_fma_double
38  * Electrostatics interaction: ReactionField
39  * VdW interaction:            LennardJones
40  * Geometry:                   Water3-Particle
41  * Calculate force/pot:        PotentialAndForce
42  */
43 void
44 nb_kernel_ElecRF_VdwLJ_GeomW3P1_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              vdwjidx0A,vdwjidx0B;
73     __m128d          jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
74     __m128d          dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
75     __m128d          dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
76     __m128d          dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
77     __m128d          velec,felec,velecsum,facel,crf,krf,krf2;
78     real             *charge;
79     int              nvdwtype;
80     __m128d          rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
81     int              *vdwtype;
82     real             *vdwparam;
83     __m128d          one_sixth   = _mm_set1_pd(1.0/6.0);
84     __m128d          one_twelfth = _mm_set1_pd(1.0/12.0);
85     __m128d          dummy_mask,cutoff_mask;
86     __m128d          signbit   = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
87     __m128d          one     = _mm_set1_pd(1.0);
88     __m128d          two     = _mm_set1_pd(2.0);
89     x                = xx[0];
90     f                = ff[0];
91
92     nri              = nlist->nri;
93     iinr             = nlist->iinr;
94     jindex           = nlist->jindex;
95     jjnr             = nlist->jjnr;
96     shiftidx         = nlist->shift;
97     gid              = nlist->gid;
98     shiftvec         = fr->shift_vec[0];
99     fshift           = fr->fshift[0];
100     facel            = _mm_set1_pd(fr->epsfac);
101     charge           = mdatoms->chargeA;
102     krf              = _mm_set1_pd(fr->ic->k_rf);
103     krf2             = _mm_set1_pd(fr->ic->k_rf*2.0);
104     crf              = _mm_set1_pd(fr->ic->c_rf);
105     nvdwtype         = fr->ntype;
106     vdwparam         = fr->nbfp;
107     vdwtype          = mdatoms->typeA;
108
109     /* Setup water-specific parameters */
110     inr              = nlist->iinr[0];
111     iq0              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+0]));
112     iq1              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+1]));
113     iq2              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+2]));
114     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
115
116     /* Avoid stupid compiler warnings */
117     jnrA = jnrB = 0;
118     j_coord_offsetA = 0;
119     j_coord_offsetB = 0;
120
121     outeriter        = 0;
122     inneriter        = 0;
123
124     /* Start outer loop over neighborlists */
125     for(iidx=0; iidx<nri; iidx++)
126     {
127         /* Load shift vector for this list */
128         i_shift_offset   = DIM*shiftidx[iidx];
129
130         /* Load limits for loop over neighbors */
131         j_index_start    = jindex[iidx];
132         j_index_end      = jindex[iidx+1];
133
134         /* Get outer coordinate index */
135         inr              = iinr[iidx];
136         i_coord_offset   = DIM*inr;
137
138         /* Load i particle coords and add shift vector */
139         gmx_mm_load_shift_and_3rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,
140                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2);
141
142         fix0             = _mm_setzero_pd();
143         fiy0             = _mm_setzero_pd();
144         fiz0             = _mm_setzero_pd();
145         fix1             = _mm_setzero_pd();
146         fiy1             = _mm_setzero_pd();
147         fiz1             = _mm_setzero_pd();
148         fix2             = _mm_setzero_pd();
149         fiy2             = _mm_setzero_pd();
150         fiz2             = _mm_setzero_pd();
151
152         /* Reset potential sums */
153         velecsum         = _mm_setzero_pd();
154         vvdwsum          = _mm_setzero_pd();
155
156         /* Start inner kernel loop */
157         for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
158         {
159
160             /* Get j neighbor index, and coordinate index */
161             jnrA             = jjnr[jidx];
162             jnrB             = jjnr[jidx+1];
163             j_coord_offsetA  = DIM*jnrA;
164             j_coord_offsetB  = DIM*jnrB;
165
166             /* load j atom coordinates */
167             gmx_mm_load_1rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
168                                               &jx0,&jy0,&jz0);
169
170             /* Calculate displacement vector */
171             dx00             = _mm_sub_pd(ix0,jx0);
172             dy00             = _mm_sub_pd(iy0,jy0);
173             dz00             = _mm_sub_pd(iz0,jz0);
174             dx10             = _mm_sub_pd(ix1,jx0);
175             dy10             = _mm_sub_pd(iy1,jy0);
176             dz10             = _mm_sub_pd(iz1,jz0);
177             dx20             = _mm_sub_pd(ix2,jx0);
178             dy20             = _mm_sub_pd(iy2,jy0);
179             dz20             = _mm_sub_pd(iz2,jz0);
180
181             /* Calculate squared distance and things based on it */
182             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
183             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
184             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
185
186             rinv00           = gmx_mm_invsqrt_pd(rsq00);
187             rinv10           = gmx_mm_invsqrt_pd(rsq10);
188             rinv20           = gmx_mm_invsqrt_pd(rsq20);
189
190             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
191             rinvsq10         = _mm_mul_pd(rinv10,rinv10);
192             rinvsq20         = _mm_mul_pd(rinv20,rinv20);
193
194             /* Load parameters for j particles */
195             jq0              = gmx_mm_load_2real_swizzle_pd(charge+jnrA+0,charge+jnrB+0);
196             vdwjidx0A        = 2*vdwtype[jnrA+0];
197             vdwjidx0B        = 2*vdwtype[jnrB+0];
198
199             fjx0             = _mm_setzero_pd();
200             fjy0             = _mm_setzero_pd();
201             fjz0             = _mm_setzero_pd();
202
203             /**************************
204              * CALCULATE INTERACTIONS *
205              **************************/
206
207             /* Compute parameters for interactions between i and j atoms */
208             qq00             = _mm_mul_pd(iq0,jq0);
209             gmx_mm_load_2pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,
210                                          vdwparam+vdwioffset0+vdwjidx0B,&c6_00,&c12_00);
211
212             /* REACTION-FIELD ELECTROSTATICS */
213             velec            = _mm_mul_pd(qq00,_mm_sub_pd(_mm_macc_pd(krf,rsq00,rinv00),crf));
214             felec            = _mm_mul_pd(qq00,_mm_msub_pd(rinv00,rinvsq00,krf2));
215
216             /* LENNARD-JONES DISPERSION/REPULSION */
217
218             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
219             vvdw6            = _mm_mul_pd(c6_00,rinvsix);
220             vvdw12           = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
221             vvdw             = _mm_msub_pd( vvdw12,one_twelfth, _mm_mul_pd(vvdw6,one_sixth) );
222             fvdw             = _mm_mul_pd(_mm_sub_pd(vvdw12,vvdw6),rinvsq00);
223
224             /* Update potential sum for this i atom from the interaction with this j atom. */
225             velecsum         = _mm_add_pd(velecsum,velec);
226             vvdwsum          = _mm_add_pd(vvdwsum,vvdw);
227
228             fscal            = _mm_add_pd(felec,fvdw);
229
230             /* Update vectorial force */
231             fix0             = _mm_macc_pd(dx00,fscal,fix0);
232             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
233             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
234             
235             fjx0             = _mm_macc_pd(dx00,fscal,fjx0);
236             fjy0             = _mm_macc_pd(dy00,fscal,fjy0);
237             fjz0             = _mm_macc_pd(dz00,fscal,fjz0);
238
239             /**************************
240              * CALCULATE INTERACTIONS *
241              **************************/
242
243             /* Compute parameters for interactions between i and j atoms */
244             qq10             = _mm_mul_pd(iq1,jq0);
245
246             /* REACTION-FIELD ELECTROSTATICS */
247             velec            = _mm_mul_pd(qq10,_mm_sub_pd(_mm_macc_pd(krf,rsq10,rinv10),crf));
248             felec            = _mm_mul_pd(qq10,_mm_msub_pd(rinv10,rinvsq10,krf2));
249
250             /* Update potential sum for this i atom from the interaction with this j atom. */
251             velecsum         = _mm_add_pd(velecsum,velec);
252
253             fscal            = felec;
254
255             /* Update vectorial force */
256             fix1             = _mm_macc_pd(dx10,fscal,fix1);
257             fiy1             = _mm_macc_pd(dy10,fscal,fiy1);
258             fiz1             = _mm_macc_pd(dz10,fscal,fiz1);
259             
260             fjx0             = _mm_macc_pd(dx10,fscal,fjx0);
261             fjy0             = _mm_macc_pd(dy10,fscal,fjy0);
262             fjz0             = _mm_macc_pd(dz10,fscal,fjz0);
263
264             /**************************
265              * CALCULATE INTERACTIONS *
266              **************************/
267
268             /* Compute parameters for interactions between i and j atoms */
269             qq20             = _mm_mul_pd(iq2,jq0);
270
271             /* REACTION-FIELD ELECTROSTATICS */
272             velec            = _mm_mul_pd(qq20,_mm_sub_pd(_mm_macc_pd(krf,rsq20,rinv20),crf));
273             felec            = _mm_mul_pd(qq20,_mm_msub_pd(rinv20,rinvsq20,krf2));
274
275             /* Update potential sum for this i atom from the interaction with this j atom. */
276             velecsum         = _mm_add_pd(velecsum,velec);
277
278             fscal            = felec;
279
280             /* Update vectorial force */
281             fix2             = _mm_macc_pd(dx20,fscal,fix2);
282             fiy2             = _mm_macc_pd(dy20,fscal,fiy2);
283             fiz2             = _mm_macc_pd(dz20,fscal,fiz2);
284             
285             fjx0             = _mm_macc_pd(dx20,fscal,fjx0);
286             fjy0             = _mm_macc_pd(dy20,fscal,fjy0);
287             fjz0             = _mm_macc_pd(dz20,fscal,fjz0);
288
289             gmx_mm_decrement_1rvec_2ptr_swizzle_pd(f+j_coord_offsetA,f+j_coord_offsetB,fjx0,fjy0,fjz0);
290
291             /* Inner loop uses 120 flops */
292         }
293
294         if(jidx<j_index_end)
295         {
296
297             jnrA             = jjnr[jidx];
298             j_coord_offsetA  = DIM*jnrA;
299
300             /* load j atom coordinates */
301             gmx_mm_load_1rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
302                                               &jx0,&jy0,&jz0);
303
304             /* Calculate displacement vector */
305             dx00             = _mm_sub_pd(ix0,jx0);
306             dy00             = _mm_sub_pd(iy0,jy0);
307             dz00             = _mm_sub_pd(iz0,jz0);
308             dx10             = _mm_sub_pd(ix1,jx0);
309             dy10             = _mm_sub_pd(iy1,jy0);
310             dz10             = _mm_sub_pd(iz1,jz0);
311             dx20             = _mm_sub_pd(ix2,jx0);
312             dy20             = _mm_sub_pd(iy2,jy0);
313             dz20             = _mm_sub_pd(iz2,jz0);
314
315             /* Calculate squared distance and things based on it */
316             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
317             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
318             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
319
320             rinv00           = gmx_mm_invsqrt_pd(rsq00);
321             rinv10           = gmx_mm_invsqrt_pd(rsq10);
322             rinv20           = gmx_mm_invsqrt_pd(rsq20);
323
324             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
325             rinvsq10         = _mm_mul_pd(rinv10,rinv10);
326             rinvsq20         = _mm_mul_pd(rinv20,rinv20);
327
328             /* Load parameters for j particles */
329             jq0              = _mm_load_sd(charge+jnrA+0);
330             vdwjidx0A        = 2*vdwtype[jnrA+0];
331
332             fjx0             = _mm_setzero_pd();
333             fjy0             = _mm_setzero_pd();
334             fjz0             = _mm_setzero_pd();
335
336             /**************************
337              * CALCULATE INTERACTIONS *
338              **************************/
339
340             /* Compute parameters for interactions between i and j atoms */
341             qq00             = _mm_mul_pd(iq0,jq0);
342             gmx_mm_load_1pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,&c6_00,&c12_00);
343
344             /* REACTION-FIELD ELECTROSTATICS */
345             velec            = _mm_mul_pd(qq00,_mm_sub_pd(_mm_macc_pd(krf,rsq00,rinv00),crf));
346             felec            = _mm_mul_pd(qq00,_mm_msub_pd(rinv00,rinvsq00,krf2));
347
348             /* LENNARD-JONES DISPERSION/REPULSION */
349
350             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
351             vvdw6            = _mm_mul_pd(c6_00,rinvsix);
352             vvdw12           = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
353             vvdw             = _mm_msub_pd( vvdw12,one_twelfth, _mm_mul_pd(vvdw6,one_sixth) );
354             fvdw             = _mm_mul_pd(_mm_sub_pd(vvdw12,vvdw6),rinvsq00);
355
356             /* Update potential sum for this i atom from the interaction with this j atom. */
357             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
358             velecsum         = _mm_add_pd(velecsum,velec);
359             vvdw             = _mm_unpacklo_pd(vvdw,_mm_setzero_pd());
360             vvdwsum          = _mm_add_pd(vvdwsum,vvdw);
361
362             fscal            = _mm_add_pd(felec,fvdw);
363
364             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
365
366             /* Update vectorial force */
367             fix0             = _mm_macc_pd(dx00,fscal,fix0);
368             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
369             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
370             
371             fjx0             = _mm_macc_pd(dx00,fscal,fjx0);
372             fjy0             = _mm_macc_pd(dy00,fscal,fjy0);
373             fjz0             = _mm_macc_pd(dz00,fscal,fjz0);
374
375             /**************************
376              * CALCULATE INTERACTIONS *
377              **************************/
378
379             /* Compute parameters for interactions between i and j atoms */
380             qq10             = _mm_mul_pd(iq1,jq0);
381
382             /* REACTION-FIELD ELECTROSTATICS */
383             velec            = _mm_mul_pd(qq10,_mm_sub_pd(_mm_macc_pd(krf,rsq10,rinv10),crf));
384             felec            = _mm_mul_pd(qq10,_mm_msub_pd(rinv10,rinvsq10,krf2));
385
386             /* Update potential sum for this i atom from the interaction with this j atom. */
387             velec            = _mm_unpacklo_pd(velec,_mm_setzero_pd());
388             velecsum         = _mm_add_pd(velecsum,velec);
389
390             fscal            = felec;
391
392             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
393
394             /* Update vectorial force */
395             fix1             = _mm_macc_pd(dx10,fscal,fix1);
396             fiy1             = _mm_macc_pd(dy10,fscal,fiy1);
397             fiz1             = _mm_macc_pd(dz10,fscal,fiz1);
398             
399             fjx0             = _mm_macc_pd(dx10,fscal,fjx0);
400             fjy0             = _mm_macc_pd(dy10,fscal,fjy0);
401             fjz0             = _mm_macc_pd(dz10,fscal,fjz0);
402
403             /**************************
404              * CALCULATE INTERACTIONS *
405              **************************/
406
407             /* Compute parameters for interactions between i and j atoms */
408             qq20             = _mm_mul_pd(iq2,jq0);
409
410             /* REACTION-FIELD ELECTROSTATICS */
411             velec            = _mm_mul_pd(qq20,_mm_sub_pd(_mm_macc_pd(krf,rsq20,rinv20),crf));
412             felec            = _mm_mul_pd(qq20,_mm_msub_pd(rinv20,rinvsq20,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             fix2             = _mm_macc_pd(dx20,fscal,fix2);
424             fiy2             = _mm_macc_pd(dy20,fscal,fiy2);
425             fiz2             = _mm_macc_pd(dz20,fscal,fiz2);
426             
427             fjx0             = _mm_macc_pd(dx20,fscal,fjx0);
428             fjy0             = _mm_macc_pd(dy20,fscal,fjy0);
429             fjz0             = _mm_macc_pd(dz20,fscal,fjz0);
430
431             gmx_mm_decrement_1rvec_1ptr_swizzle_pd(f+j_coord_offsetA,fjx0,fjy0,fjz0);
432
433             /* Inner loop uses 120 flops */
434         }
435
436         /* End of innermost loop */
437
438         gmx_mm_update_iforce_3atom_swizzle_pd(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,
439                                               f+i_coord_offset,fshift+i_shift_offset);
440
441         ggid                        = gid[iidx];
442         /* Update potential energies */
443         gmx_mm_update_1pot_pd(velecsum,kernel_data->energygrp_elec+ggid);
444         gmx_mm_update_1pot_pd(vvdwsum,kernel_data->energygrp_vdw+ggid);
445
446         /* Increment number of inner iterations */
447         inneriter                  += j_index_end - j_index_start;
448
449         /* Outer loop uses 20 flops */
450     }
451
452     /* Increment number of outer iterations */
453     outeriter        += nri;
454
455     /* Update outer/inner flops */
456
457     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W3_VF,outeriter*20 + inneriter*120);
458 }
459 /*
460  * Gromacs nonbonded kernel:   nb_kernel_ElecRF_VdwLJ_GeomW3P1_F_avx_128_fma_double
461  * Electrostatics interaction: ReactionField
462  * VdW interaction:            LennardJones
463  * Geometry:                   Water3-Particle
464  * Calculate force/pot:        Force
465  */
466 void
467 nb_kernel_ElecRF_VdwLJ_GeomW3P1_F_avx_128_fma_double
468                     (t_nblist * gmx_restrict                nlist,
469                      rvec * gmx_restrict                    xx,
470                      rvec * gmx_restrict                    ff,
471                      t_forcerec * gmx_restrict              fr,
472                      t_mdatoms * gmx_restrict               mdatoms,
473                      nb_kernel_data_t * gmx_restrict        kernel_data,
474                      t_nrnb * gmx_restrict                  nrnb)
475 {
476     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
477      * just 0 for non-waters.
478      * Suffixes A,B refer to j loop unrolling done with SSE double precision, e.g. for the two different
479      * jnr indices corresponding to data put in the four positions in the SIMD register.
480      */
481     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
482     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
483     int              jnrA,jnrB;
484     int              j_coord_offsetA,j_coord_offsetB;
485     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
486     real             rcutoff_scalar;
487     real             *shiftvec,*fshift,*x,*f;
488     __m128d          tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
489     int              vdwioffset0;
490     __m128d          ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
491     int              vdwioffset1;
492     __m128d          ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
493     int              vdwioffset2;
494     __m128d          ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
495     int              vdwjidx0A,vdwjidx0B;
496     __m128d          jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
497     __m128d          dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
498     __m128d          dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
499     __m128d          dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
500     __m128d          velec,felec,velecsum,facel,crf,krf,krf2;
501     real             *charge;
502     int              nvdwtype;
503     __m128d          rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
504     int              *vdwtype;
505     real             *vdwparam;
506     __m128d          one_sixth   = _mm_set1_pd(1.0/6.0);
507     __m128d          one_twelfth = _mm_set1_pd(1.0/12.0);
508     __m128d          dummy_mask,cutoff_mask;
509     __m128d          signbit   = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
510     __m128d          one     = _mm_set1_pd(1.0);
511     __m128d          two     = _mm_set1_pd(2.0);
512     x                = xx[0];
513     f                = ff[0];
514
515     nri              = nlist->nri;
516     iinr             = nlist->iinr;
517     jindex           = nlist->jindex;
518     jjnr             = nlist->jjnr;
519     shiftidx         = nlist->shift;
520     gid              = nlist->gid;
521     shiftvec         = fr->shift_vec[0];
522     fshift           = fr->fshift[0];
523     facel            = _mm_set1_pd(fr->epsfac);
524     charge           = mdatoms->chargeA;
525     krf              = _mm_set1_pd(fr->ic->k_rf);
526     krf2             = _mm_set1_pd(fr->ic->k_rf*2.0);
527     crf              = _mm_set1_pd(fr->ic->c_rf);
528     nvdwtype         = fr->ntype;
529     vdwparam         = fr->nbfp;
530     vdwtype          = mdatoms->typeA;
531
532     /* Setup water-specific parameters */
533     inr              = nlist->iinr[0];
534     iq0              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+0]));
535     iq1              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+1]));
536     iq2              = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+2]));
537     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
538
539     /* Avoid stupid compiler warnings */
540     jnrA = jnrB = 0;
541     j_coord_offsetA = 0;
542     j_coord_offsetB = 0;
543
544     outeriter        = 0;
545     inneriter        = 0;
546
547     /* Start outer loop over neighborlists */
548     for(iidx=0; iidx<nri; iidx++)
549     {
550         /* Load shift vector for this list */
551         i_shift_offset   = DIM*shiftidx[iidx];
552
553         /* Load limits for loop over neighbors */
554         j_index_start    = jindex[iidx];
555         j_index_end      = jindex[iidx+1];
556
557         /* Get outer coordinate index */
558         inr              = iinr[iidx];
559         i_coord_offset   = DIM*inr;
560
561         /* Load i particle coords and add shift vector */
562         gmx_mm_load_shift_and_3rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,
563                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2);
564
565         fix0             = _mm_setzero_pd();
566         fiy0             = _mm_setzero_pd();
567         fiz0             = _mm_setzero_pd();
568         fix1             = _mm_setzero_pd();
569         fiy1             = _mm_setzero_pd();
570         fiz1             = _mm_setzero_pd();
571         fix2             = _mm_setzero_pd();
572         fiy2             = _mm_setzero_pd();
573         fiz2             = _mm_setzero_pd();
574
575         /* Start inner kernel loop */
576         for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
577         {
578
579             /* Get j neighbor index, and coordinate index */
580             jnrA             = jjnr[jidx];
581             jnrB             = jjnr[jidx+1];
582             j_coord_offsetA  = DIM*jnrA;
583             j_coord_offsetB  = DIM*jnrB;
584
585             /* load j atom coordinates */
586             gmx_mm_load_1rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
587                                               &jx0,&jy0,&jz0);
588
589             /* Calculate displacement vector */
590             dx00             = _mm_sub_pd(ix0,jx0);
591             dy00             = _mm_sub_pd(iy0,jy0);
592             dz00             = _mm_sub_pd(iz0,jz0);
593             dx10             = _mm_sub_pd(ix1,jx0);
594             dy10             = _mm_sub_pd(iy1,jy0);
595             dz10             = _mm_sub_pd(iz1,jz0);
596             dx20             = _mm_sub_pd(ix2,jx0);
597             dy20             = _mm_sub_pd(iy2,jy0);
598             dz20             = _mm_sub_pd(iz2,jz0);
599
600             /* Calculate squared distance and things based on it */
601             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
602             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
603             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
604
605             rinv00           = gmx_mm_invsqrt_pd(rsq00);
606             rinv10           = gmx_mm_invsqrt_pd(rsq10);
607             rinv20           = gmx_mm_invsqrt_pd(rsq20);
608
609             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
610             rinvsq10         = _mm_mul_pd(rinv10,rinv10);
611             rinvsq20         = _mm_mul_pd(rinv20,rinv20);
612
613             /* Load parameters for j particles */
614             jq0              = gmx_mm_load_2real_swizzle_pd(charge+jnrA+0,charge+jnrB+0);
615             vdwjidx0A        = 2*vdwtype[jnrA+0];
616             vdwjidx0B        = 2*vdwtype[jnrB+0];
617
618             fjx0             = _mm_setzero_pd();
619             fjy0             = _mm_setzero_pd();
620             fjz0             = _mm_setzero_pd();
621
622             /**************************
623              * CALCULATE INTERACTIONS *
624              **************************/
625
626             /* Compute parameters for interactions between i and j atoms */
627             qq00             = _mm_mul_pd(iq0,jq0);
628             gmx_mm_load_2pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,
629                                          vdwparam+vdwioffset0+vdwjidx0B,&c6_00,&c12_00);
630
631             /* REACTION-FIELD ELECTROSTATICS */
632             felec            = _mm_mul_pd(qq00,_mm_msub_pd(rinv00,rinvsq00,krf2));
633
634             /* LENNARD-JONES DISPERSION/REPULSION */
635
636             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
637             fvdw             = _mm_mul_pd(_mm_msub_pd(c12_00,rinvsix,c6_00),_mm_mul_pd(rinvsix,rinvsq00));
638
639             fscal            = _mm_add_pd(felec,fvdw);
640
641             /* Update vectorial force */
642             fix0             = _mm_macc_pd(dx00,fscal,fix0);
643             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
644             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
645             
646             fjx0             = _mm_macc_pd(dx00,fscal,fjx0);
647             fjy0             = _mm_macc_pd(dy00,fscal,fjy0);
648             fjz0             = _mm_macc_pd(dz00,fscal,fjz0);
649
650             /**************************
651              * CALCULATE INTERACTIONS *
652              **************************/
653
654             /* Compute parameters for interactions between i and j atoms */
655             qq10             = _mm_mul_pd(iq1,jq0);
656
657             /* REACTION-FIELD ELECTROSTATICS */
658             felec            = _mm_mul_pd(qq10,_mm_msub_pd(rinv10,rinvsq10,krf2));
659
660             fscal            = felec;
661
662             /* Update vectorial force */
663             fix1             = _mm_macc_pd(dx10,fscal,fix1);
664             fiy1             = _mm_macc_pd(dy10,fscal,fiy1);
665             fiz1             = _mm_macc_pd(dz10,fscal,fiz1);
666             
667             fjx0             = _mm_macc_pd(dx10,fscal,fjx0);
668             fjy0             = _mm_macc_pd(dy10,fscal,fjy0);
669             fjz0             = _mm_macc_pd(dz10,fscal,fjz0);
670
671             /**************************
672              * CALCULATE INTERACTIONS *
673              **************************/
674
675             /* Compute parameters for interactions between i and j atoms */
676             qq20             = _mm_mul_pd(iq2,jq0);
677
678             /* REACTION-FIELD ELECTROSTATICS */
679             felec            = _mm_mul_pd(qq20,_mm_msub_pd(rinv20,rinvsq20,krf2));
680
681             fscal            = felec;
682
683             /* Update vectorial force */
684             fix2             = _mm_macc_pd(dx20,fscal,fix2);
685             fiy2             = _mm_macc_pd(dy20,fscal,fiy2);
686             fiz2             = _mm_macc_pd(dz20,fscal,fiz2);
687             
688             fjx0             = _mm_macc_pd(dx20,fscal,fjx0);
689             fjy0             = _mm_macc_pd(dy20,fscal,fjy0);
690             fjz0             = _mm_macc_pd(dz20,fscal,fjz0);
691
692             gmx_mm_decrement_1rvec_2ptr_swizzle_pd(f+j_coord_offsetA,f+j_coord_offsetB,fjx0,fjy0,fjz0);
693
694             /* Inner loop uses 100 flops */
695         }
696
697         if(jidx<j_index_end)
698         {
699
700             jnrA             = jjnr[jidx];
701             j_coord_offsetA  = DIM*jnrA;
702
703             /* load j atom coordinates */
704             gmx_mm_load_1rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
705                                               &jx0,&jy0,&jz0);
706
707             /* Calculate displacement vector */
708             dx00             = _mm_sub_pd(ix0,jx0);
709             dy00             = _mm_sub_pd(iy0,jy0);
710             dz00             = _mm_sub_pd(iz0,jz0);
711             dx10             = _mm_sub_pd(ix1,jx0);
712             dy10             = _mm_sub_pd(iy1,jy0);
713             dz10             = _mm_sub_pd(iz1,jz0);
714             dx20             = _mm_sub_pd(ix2,jx0);
715             dy20             = _mm_sub_pd(iy2,jy0);
716             dz20             = _mm_sub_pd(iz2,jz0);
717
718             /* Calculate squared distance and things based on it */
719             rsq00            = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
720             rsq10            = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
721             rsq20            = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
722
723             rinv00           = gmx_mm_invsqrt_pd(rsq00);
724             rinv10           = gmx_mm_invsqrt_pd(rsq10);
725             rinv20           = gmx_mm_invsqrt_pd(rsq20);
726
727             rinvsq00         = _mm_mul_pd(rinv00,rinv00);
728             rinvsq10         = _mm_mul_pd(rinv10,rinv10);
729             rinvsq20         = _mm_mul_pd(rinv20,rinv20);
730
731             /* Load parameters for j particles */
732             jq0              = _mm_load_sd(charge+jnrA+0);
733             vdwjidx0A        = 2*vdwtype[jnrA+0];
734
735             fjx0             = _mm_setzero_pd();
736             fjy0             = _mm_setzero_pd();
737             fjz0             = _mm_setzero_pd();
738
739             /**************************
740              * CALCULATE INTERACTIONS *
741              **************************/
742
743             /* Compute parameters for interactions between i and j atoms */
744             qq00             = _mm_mul_pd(iq0,jq0);
745             gmx_mm_load_1pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,&c6_00,&c12_00);
746
747             /* REACTION-FIELD ELECTROSTATICS */
748             felec            = _mm_mul_pd(qq00,_mm_msub_pd(rinv00,rinvsq00,krf2));
749
750             /* LENNARD-JONES DISPERSION/REPULSION */
751
752             rinvsix          = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
753             fvdw             = _mm_mul_pd(_mm_msub_pd(c12_00,rinvsix,c6_00),_mm_mul_pd(rinvsix,rinvsq00));
754
755             fscal            = _mm_add_pd(felec,fvdw);
756
757             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
758
759             /* Update vectorial force */
760             fix0             = _mm_macc_pd(dx00,fscal,fix0);
761             fiy0             = _mm_macc_pd(dy00,fscal,fiy0);
762             fiz0             = _mm_macc_pd(dz00,fscal,fiz0);
763             
764             fjx0             = _mm_macc_pd(dx00,fscal,fjx0);
765             fjy0             = _mm_macc_pd(dy00,fscal,fjy0);
766             fjz0             = _mm_macc_pd(dz00,fscal,fjz0);
767
768             /**************************
769              * CALCULATE INTERACTIONS *
770              **************************/
771
772             /* Compute parameters for interactions between i and j atoms */
773             qq10             = _mm_mul_pd(iq1,jq0);
774
775             /* REACTION-FIELD ELECTROSTATICS */
776             felec            = _mm_mul_pd(qq10,_mm_msub_pd(rinv10,rinvsq10,krf2));
777
778             fscal            = felec;
779
780             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
781
782             /* Update vectorial force */
783             fix1             = _mm_macc_pd(dx10,fscal,fix1);
784             fiy1             = _mm_macc_pd(dy10,fscal,fiy1);
785             fiz1             = _mm_macc_pd(dz10,fscal,fiz1);
786             
787             fjx0             = _mm_macc_pd(dx10,fscal,fjx0);
788             fjy0             = _mm_macc_pd(dy10,fscal,fjy0);
789             fjz0             = _mm_macc_pd(dz10,fscal,fjz0);
790
791             /**************************
792              * CALCULATE INTERACTIONS *
793              **************************/
794
795             /* Compute parameters for interactions between i and j atoms */
796             qq20             = _mm_mul_pd(iq2,jq0);
797
798             /* REACTION-FIELD ELECTROSTATICS */
799             felec            = _mm_mul_pd(qq20,_mm_msub_pd(rinv20,rinvsq20,krf2));
800
801             fscal            = felec;
802
803             fscal            = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
804
805             /* Update vectorial force */
806             fix2             = _mm_macc_pd(dx20,fscal,fix2);
807             fiy2             = _mm_macc_pd(dy20,fscal,fiy2);
808             fiz2             = _mm_macc_pd(dz20,fscal,fiz2);
809             
810             fjx0             = _mm_macc_pd(dx20,fscal,fjx0);
811             fjy0             = _mm_macc_pd(dy20,fscal,fjy0);
812             fjz0             = _mm_macc_pd(dz20,fscal,fjz0);
813
814             gmx_mm_decrement_1rvec_1ptr_swizzle_pd(f+j_coord_offsetA,fjx0,fjy0,fjz0);
815
816             /* Inner loop uses 100 flops */
817         }
818
819         /* End of innermost loop */
820
821         gmx_mm_update_iforce_3atom_swizzle_pd(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,
822                                               f+i_coord_offset,fshift+i_shift_offset);
823
824         /* Increment number of inner iterations */
825         inneriter                  += j_index_end - j_index_start;
826
827         /* Outer loop uses 18 flops */
828     }
829
830     /* Increment number of outer iterations */
831     outeriter        += nri;
832
833     /* Update outer/inner flops */
834
835     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W3_F,outeriter*18 + inneriter*100);
836 }