2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013,2014,2015,2017,2018, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
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.
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.
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.
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.
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.
36 * Note: this file was generated by the GROMACS avx_128_fma_double kernel generator.
44 #include "../nb_kernel.h"
45 #include "gromacs/gmxlib/nrnb.h"
47 #include "kernelutil_x86_avx_128_fma_double.h"
50 * Gromacs nonbonded kernel: nb_kernel_ElecEwSh_VdwLJEwSh_GeomW4P1_VF_avx_128_fma_double
51 * Electrostatics interaction: Ewald
52 * VdW interaction: LJEwald
53 * Geometry: Water4-Particle
54 * Calculate force/pot: PotentialAndForce
57 nb_kernel_ElecEwSh_VdwLJEwSh_GeomW4P1_VF_avx_128_fma_double
58 (t_nblist * gmx_restrict nlist,
59 rvec * gmx_restrict xx,
60 rvec * gmx_restrict ff,
61 struct t_forcerec * gmx_restrict fr,
62 t_mdatoms * gmx_restrict mdatoms,
63 nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
64 t_nrnb * gmx_restrict nrnb)
66 /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
67 * just 0 for non-waters.
68 * Suffixes A,B refer to j loop unrolling done with SSE double precision, e.g. for the two different
69 * jnr indices corresponding to data put in the four positions in the SIMD register.
71 int i_shift_offset,i_coord_offset,outeriter,inneriter;
72 int j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
74 int j_coord_offsetA,j_coord_offsetB;
75 int *iinr,*jindex,*jjnr,*shiftidx,*gid;
77 real *shiftvec,*fshift,*x,*f;
78 __m128d tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
80 __m128d ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
82 __m128d ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
84 __m128d ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
86 __m128d ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
87 int vdwjidx0A,vdwjidx0B;
88 __m128d jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
89 __m128d dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
90 __m128d dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
91 __m128d dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
92 __m128d dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
93 __m128d velec,felec,velecsum,facel,crf,krf,krf2;
96 __m128d rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
99 __m128d one_sixth = _mm_set1_pd(1.0/6.0);
100 __m128d one_twelfth = _mm_set1_pd(1.0/12.0);
106 __m128d ewclj,ewclj2,ewclj6,ewcljrsq,poly,exponent,f6A,f6B,sh_lj_ewald;
107 __m128d one_half = _mm_set1_pd(0.5);
108 __m128d minus_one = _mm_set1_pd(-1.0);
110 __m128d ewtabscale,eweps,twoeweps,sh_ewald,ewrt,ewtabhalfspace,ewtabF,ewtabFn,ewtabD,ewtabV;
112 __m128d dummy_mask,cutoff_mask;
113 __m128d signbit = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
114 __m128d one = _mm_set1_pd(1.0);
115 __m128d two = _mm_set1_pd(2.0);
121 jindex = nlist->jindex;
123 shiftidx = nlist->shift;
125 shiftvec = fr->shift_vec[0];
126 fshift = fr->fshift[0];
127 facel = _mm_set1_pd(fr->ic->epsfac);
128 charge = mdatoms->chargeA;
129 nvdwtype = fr->ntype;
131 vdwtype = mdatoms->typeA;
132 vdwgridparam = fr->ljpme_c6grid;
133 sh_lj_ewald = _mm_set1_pd(fr->ic->sh_lj_ewald);
134 ewclj = _mm_set1_pd(fr->ic->ewaldcoeff_lj);
135 ewclj2 = _mm_mul_pd(minus_one,_mm_mul_pd(ewclj,ewclj));
137 sh_ewald = _mm_set1_pd(fr->ic->sh_ewald);
138 ewtab = fr->ic->tabq_coul_FDV0;
139 ewtabscale = _mm_set1_pd(fr->ic->tabq_scale);
140 ewtabhalfspace = _mm_set1_pd(0.5/fr->ic->tabq_scale);
142 /* Setup water-specific parameters */
143 inr = nlist->iinr[0];
144 iq1 = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+1]));
145 iq2 = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+2]));
146 iq3 = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+3]));
147 vdwioffset0 = 2*nvdwtype*vdwtype[inr+0];
149 /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
150 rcutoff_scalar = fr->ic->rcoulomb;
151 rcutoff = _mm_set1_pd(rcutoff_scalar);
152 rcutoff2 = _mm_mul_pd(rcutoff,rcutoff);
154 sh_vdw_invrcut6 = _mm_set1_pd(fr->ic->sh_invrc6);
155 rvdw = _mm_set1_pd(fr->ic->rvdw);
157 /* Avoid stupid compiler warnings */
165 /* Start outer loop over neighborlists */
166 for(iidx=0; iidx<nri; iidx++)
168 /* Load shift vector for this list */
169 i_shift_offset = DIM*shiftidx[iidx];
171 /* Load limits for loop over neighbors */
172 j_index_start = jindex[iidx];
173 j_index_end = jindex[iidx+1];
175 /* Get outer coordinate index */
177 i_coord_offset = DIM*inr;
179 /* Load i particle coords and add shift vector */
180 gmx_mm_load_shift_and_4rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,
181 &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
183 fix0 = _mm_setzero_pd();
184 fiy0 = _mm_setzero_pd();
185 fiz0 = _mm_setzero_pd();
186 fix1 = _mm_setzero_pd();
187 fiy1 = _mm_setzero_pd();
188 fiz1 = _mm_setzero_pd();
189 fix2 = _mm_setzero_pd();
190 fiy2 = _mm_setzero_pd();
191 fiz2 = _mm_setzero_pd();
192 fix3 = _mm_setzero_pd();
193 fiy3 = _mm_setzero_pd();
194 fiz3 = _mm_setzero_pd();
196 /* Reset potential sums */
197 velecsum = _mm_setzero_pd();
198 vvdwsum = _mm_setzero_pd();
200 /* Start inner kernel loop */
201 for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
204 /* Get j neighbor index, and coordinate index */
207 j_coord_offsetA = DIM*jnrA;
208 j_coord_offsetB = DIM*jnrB;
210 /* load j atom coordinates */
211 gmx_mm_load_1rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
214 /* Calculate displacement vector */
215 dx00 = _mm_sub_pd(ix0,jx0);
216 dy00 = _mm_sub_pd(iy0,jy0);
217 dz00 = _mm_sub_pd(iz0,jz0);
218 dx10 = _mm_sub_pd(ix1,jx0);
219 dy10 = _mm_sub_pd(iy1,jy0);
220 dz10 = _mm_sub_pd(iz1,jz0);
221 dx20 = _mm_sub_pd(ix2,jx0);
222 dy20 = _mm_sub_pd(iy2,jy0);
223 dz20 = _mm_sub_pd(iz2,jz0);
224 dx30 = _mm_sub_pd(ix3,jx0);
225 dy30 = _mm_sub_pd(iy3,jy0);
226 dz30 = _mm_sub_pd(iz3,jz0);
228 /* Calculate squared distance and things based on it */
229 rsq00 = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
230 rsq10 = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
231 rsq20 = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
232 rsq30 = gmx_mm_calc_rsq_pd(dx30,dy30,dz30);
234 rinv00 = avx128fma_invsqrt_d(rsq00);
235 rinv10 = avx128fma_invsqrt_d(rsq10);
236 rinv20 = avx128fma_invsqrt_d(rsq20);
237 rinv30 = avx128fma_invsqrt_d(rsq30);
239 rinvsq00 = _mm_mul_pd(rinv00,rinv00);
240 rinvsq10 = _mm_mul_pd(rinv10,rinv10);
241 rinvsq20 = _mm_mul_pd(rinv20,rinv20);
242 rinvsq30 = _mm_mul_pd(rinv30,rinv30);
244 /* Load parameters for j particles */
245 jq0 = gmx_mm_load_2real_swizzle_pd(charge+jnrA+0,charge+jnrB+0);
246 vdwjidx0A = 2*vdwtype[jnrA+0];
247 vdwjidx0B = 2*vdwtype[jnrB+0];
249 fjx0 = _mm_setzero_pd();
250 fjy0 = _mm_setzero_pd();
251 fjz0 = _mm_setzero_pd();
253 /**************************
254 * CALCULATE INTERACTIONS *
255 **************************/
257 if (gmx_mm_any_lt(rsq00,rcutoff2))
260 r00 = _mm_mul_pd(rsq00,rinv00);
262 /* Compute parameters for interactions between i and j atoms */
263 gmx_mm_load_2pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,
264 vdwparam+vdwioffset0+vdwjidx0B,&c6_00,&c12_00);
265 c6grid_00 = gmx_mm_load_2real_swizzle_pd(vdwgridparam+vdwioffset0+vdwjidx0A,
266 vdwgridparam+vdwioffset0+vdwjidx0B);
268 /* Analytical LJ-PME */
269 rinvsix = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
270 ewcljrsq = _mm_mul_pd(ewclj2,rsq00);
271 ewclj6 = _mm_mul_pd(ewclj2,_mm_mul_pd(ewclj2,ewclj2));
272 exponent = avx128fma_exp_d(ewcljrsq);
273 /* poly = exp(-(beta*r)^2) * (1 + (beta*r)^2 + (beta*r)^4 /2) */
274 poly = _mm_mul_pd(exponent,_mm_macc_pd(_mm_mul_pd(ewcljrsq,ewcljrsq),one_half,_mm_sub_pd(one,ewcljrsq)));
275 /* vvdw6 = [C6 - C6grid * (1-poly)]/r6 */
276 vvdw6 = _mm_mul_pd(_mm_macc_pd(-c6grid_00,_mm_sub_pd(one,poly),c6_00),rinvsix);
277 vvdw12 = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
278 vvdw = _mm_msub_pd(_mm_nmacc_pd(c12_00,_mm_mul_pd(sh_vdw_invrcut6,sh_vdw_invrcut6),vvdw12),one_twelfth,
279 _mm_mul_pd(_mm_sub_pd(vvdw6,_mm_macc_pd(c6grid_00,sh_lj_ewald,_mm_mul_pd(c6_00,sh_vdw_invrcut6))),one_sixth));
280 /* fvdw = vvdw12/r - (vvdw6/r + (C6grid * exponent * beta^6)/r) */
281 fvdw = _mm_mul_pd(_mm_add_pd(vvdw12,_mm_msub_pd(_mm_mul_pd(c6grid_00,one_sixth),_mm_mul_pd(exponent,ewclj6),vvdw6)),rinvsq00);
283 cutoff_mask = _mm_cmplt_pd(rsq00,rcutoff2);
285 /* Update potential sum for this i atom from the interaction with this j atom. */
286 vvdw = _mm_and_pd(vvdw,cutoff_mask);
287 vvdwsum = _mm_add_pd(vvdwsum,vvdw);
291 fscal = _mm_and_pd(fscal,cutoff_mask);
293 /* Update vectorial force */
294 fix0 = _mm_macc_pd(dx00,fscal,fix0);
295 fiy0 = _mm_macc_pd(dy00,fscal,fiy0);
296 fiz0 = _mm_macc_pd(dz00,fscal,fiz0);
298 fjx0 = _mm_macc_pd(dx00,fscal,fjx0);
299 fjy0 = _mm_macc_pd(dy00,fscal,fjy0);
300 fjz0 = _mm_macc_pd(dz00,fscal,fjz0);
304 /**************************
305 * CALCULATE INTERACTIONS *
306 **************************/
308 if (gmx_mm_any_lt(rsq10,rcutoff2))
311 r10 = _mm_mul_pd(rsq10,rinv10);
313 /* Compute parameters for interactions between i and j atoms */
314 qq10 = _mm_mul_pd(iq1,jq0);
316 /* EWALD ELECTROSTATICS */
318 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
319 ewrt = _mm_mul_pd(r10,ewtabscale);
320 ewitab = _mm_cvttpd_epi32(ewrt);
322 eweps = _mm_frcz_pd(ewrt);
324 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
326 twoeweps = _mm_add_pd(eweps,eweps);
327 ewitab = _mm_slli_epi32(ewitab,2);
328 ewtabF = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
329 ewtabD = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
330 GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
331 ewtabV = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
332 ewtabFn = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
333 GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
334 felec = _mm_macc_pd(eweps,ewtabD,ewtabF);
335 velec = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
336 velec = _mm_mul_pd(qq10,_mm_sub_pd(_mm_sub_pd(rinv10,sh_ewald),velec));
337 felec = _mm_mul_pd(_mm_mul_pd(qq10,rinv10),_mm_sub_pd(rinvsq10,felec));
339 cutoff_mask = _mm_cmplt_pd(rsq10,rcutoff2);
341 /* Update potential sum for this i atom from the interaction with this j atom. */
342 velec = _mm_and_pd(velec,cutoff_mask);
343 velecsum = _mm_add_pd(velecsum,velec);
347 fscal = _mm_and_pd(fscal,cutoff_mask);
349 /* Update vectorial force */
350 fix1 = _mm_macc_pd(dx10,fscal,fix1);
351 fiy1 = _mm_macc_pd(dy10,fscal,fiy1);
352 fiz1 = _mm_macc_pd(dz10,fscal,fiz1);
354 fjx0 = _mm_macc_pd(dx10,fscal,fjx0);
355 fjy0 = _mm_macc_pd(dy10,fscal,fjy0);
356 fjz0 = _mm_macc_pd(dz10,fscal,fjz0);
360 /**************************
361 * CALCULATE INTERACTIONS *
362 **************************/
364 if (gmx_mm_any_lt(rsq20,rcutoff2))
367 r20 = _mm_mul_pd(rsq20,rinv20);
369 /* Compute parameters for interactions between i and j atoms */
370 qq20 = _mm_mul_pd(iq2,jq0);
372 /* EWALD ELECTROSTATICS */
374 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
375 ewrt = _mm_mul_pd(r20,ewtabscale);
376 ewitab = _mm_cvttpd_epi32(ewrt);
378 eweps = _mm_frcz_pd(ewrt);
380 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
382 twoeweps = _mm_add_pd(eweps,eweps);
383 ewitab = _mm_slli_epi32(ewitab,2);
384 ewtabF = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
385 ewtabD = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
386 GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
387 ewtabV = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
388 ewtabFn = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
389 GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
390 felec = _mm_macc_pd(eweps,ewtabD,ewtabF);
391 velec = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
392 velec = _mm_mul_pd(qq20,_mm_sub_pd(_mm_sub_pd(rinv20,sh_ewald),velec));
393 felec = _mm_mul_pd(_mm_mul_pd(qq20,rinv20),_mm_sub_pd(rinvsq20,felec));
395 cutoff_mask = _mm_cmplt_pd(rsq20,rcutoff2);
397 /* Update potential sum for this i atom from the interaction with this j atom. */
398 velec = _mm_and_pd(velec,cutoff_mask);
399 velecsum = _mm_add_pd(velecsum,velec);
403 fscal = _mm_and_pd(fscal,cutoff_mask);
405 /* Update vectorial force */
406 fix2 = _mm_macc_pd(dx20,fscal,fix2);
407 fiy2 = _mm_macc_pd(dy20,fscal,fiy2);
408 fiz2 = _mm_macc_pd(dz20,fscal,fiz2);
410 fjx0 = _mm_macc_pd(dx20,fscal,fjx0);
411 fjy0 = _mm_macc_pd(dy20,fscal,fjy0);
412 fjz0 = _mm_macc_pd(dz20,fscal,fjz0);
416 /**************************
417 * CALCULATE INTERACTIONS *
418 **************************/
420 if (gmx_mm_any_lt(rsq30,rcutoff2))
423 r30 = _mm_mul_pd(rsq30,rinv30);
425 /* Compute parameters for interactions between i and j atoms */
426 qq30 = _mm_mul_pd(iq3,jq0);
428 /* EWALD ELECTROSTATICS */
430 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
431 ewrt = _mm_mul_pd(r30,ewtabscale);
432 ewitab = _mm_cvttpd_epi32(ewrt);
434 eweps = _mm_frcz_pd(ewrt);
436 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
438 twoeweps = _mm_add_pd(eweps,eweps);
439 ewitab = _mm_slli_epi32(ewitab,2);
440 ewtabF = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
441 ewtabD = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,1) );
442 GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
443 ewtabV = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
444 ewtabFn = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,1) +2);
445 GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
446 felec = _mm_macc_pd(eweps,ewtabD,ewtabF);
447 velec = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
448 velec = _mm_mul_pd(qq30,_mm_sub_pd(_mm_sub_pd(rinv30,sh_ewald),velec));
449 felec = _mm_mul_pd(_mm_mul_pd(qq30,rinv30),_mm_sub_pd(rinvsq30,felec));
451 cutoff_mask = _mm_cmplt_pd(rsq30,rcutoff2);
453 /* Update potential sum for this i atom from the interaction with this j atom. */
454 velec = _mm_and_pd(velec,cutoff_mask);
455 velecsum = _mm_add_pd(velecsum,velec);
459 fscal = _mm_and_pd(fscal,cutoff_mask);
461 /* Update vectorial force */
462 fix3 = _mm_macc_pd(dx30,fscal,fix3);
463 fiy3 = _mm_macc_pd(dy30,fscal,fiy3);
464 fiz3 = _mm_macc_pd(dz30,fscal,fiz3);
466 fjx0 = _mm_macc_pd(dx30,fscal,fjx0);
467 fjy0 = _mm_macc_pd(dy30,fscal,fjy0);
468 fjz0 = _mm_macc_pd(dz30,fscal,fjz0);
472 gmx_mm_decrement_1rvec_2ptr_swizzle_pd(f+j_coord_offsetA,f+j_coord_offsetB,fjx0,fjy0,fjz0);
474 /* Inner loop uses 208 flops */
481 j_coord_offsetA = DIM*jnrA;
483 /* load j atom coordinates */
484 gmx_mm_load_1rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
487 /* Calculate displacement vector */
488 dx00 = _mm_sub_pd(ix0,jx0);
489 dy00 = _mm_sub_pd(iy0,jy0);
490 dz00 = _mm_sub_pd(iz0,jz0);
491 dx10 = _mm_sub_pd(ix1,jx0);
492 dy10 = _mm_sub_pd(iy1,jy0);
493 dz10 = _mm_sub_pd(iz1,jz0);
494 dx20 = _mm_sub_pd(ix2,jx0);
495 dy20 = _mm_sub_pd(iy2,jy0);
496 dz20 = _mm_sub_pd(iz2,jz0);
497 dx30 = _mm_sub_pd(ix3,jx0);
498 dy30 = _mm_sub_pd(iy3,jy0);
499 dz30 = _mm_sub_pd(iz3,jz0);
501 /* Calculate squared distance and things based on it */
502 rsq00 = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
503 rsq10 = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
504 rsq20 = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
505 rsq30 = gmx_mm_calc_rsq_pd(dx30,dy30,dz30);
507 rinv00 = avx128fma_invsqrt_d(rsq00);
508 rinv10 = avx128fma_invsqrt_d(rsq10);
509 rinv20 = avx128fma_invsqrt_d(rsq20);
510 rinv30 = avx128fma_invsqrt_d(rsq30);
512 rinvsq00 = _mm_mul_pd(rinv00,rinv00);
513 rinvsq10 = _mm_mul_pd(rinv10,rinv10);
514 rinvsq20 = _mm_mul_pd(rinv20,rinv20);
515 rinvsq30 = _mm_mul_pd(rinv30,rinv30);
517 /* Load parameters for j particles */
518 jq0 = _mm_load_sd(charge+jnrA+0);
519 vdwjidx0A = 2*vdwtype[jnrA+0];
521 fjx0 = _mm_setzero_pd();
522 fjy0 = _mm_setzero_pd();
523 fjz0 = _mm_setzero_pd();
525 /**************************
526 * CALCULATE INTERACTIONS *
527 **************************/
529 if (gmx_mm_any_lt(rsq00,rcutoff2))
532 r00 = _mm_mul_pd(rsq00,rinv00);
534 /* Compute parameters for interactions between i and j atoms */
535 gmx_mm_load_1pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,&c6_00,&c12_00);
536 c6grid_00 = gmx_mm_load_1real_pd(vdwgridparam+vdwioffset0+vdwjidx0A);
538 /* Analytical LJ-PME */
539 rinvsix = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
540 ewcljrsq = _mm_mul_pd(ewclj2,rsq00);
541 ewclj6 = _mm_mul_pd(ewclj2,_mm_mul_pd(ewclj2,ewclj2));
542 exponent = avx128fma_exp_d(ewcljrsq);
543 /* poly = exp(-(beta*r)^2) * (1 + (beta*r)^2 + (beta*r)^4 /2) */
544 poly = _mm_mul_pd(exponent,_mm_macc_pd(_mm_mul_pd(ewcljrsq,ewcljrsq),one_half,_mm_sub_pd(one,ewcljrsq)));
545 /* vvdw6 = [C6 - C6grid * (1-poly)]/r6 */
546 vvdw6 = _mm_mul_pd(_mm_macc_pd(-c6grid_00,_mm_sub_pd(one,poly),c6_00),rinvsix);
547 vvdw12 = _mm_mul_pd(c12_00,_mm_mul_pd(rinvsix,rinvsix));
548 vvdw = _mm_msub_pd(_mm_nmacc_pd(c12_00,_mm_mul_pd(sh_vdw_invrcut6,sh_vdw_invrcut6),vvdw12),one_twelfth,
549 _mm_mul_pd(_mm_sub_pd(vvdw6,_mm_macc_pd(c6grid_00,sh_lj_ewald,_mm_mul_pd(c6_00,sh_vdw_invrcut6))),one_sixth));
550 /* fvdw = vvdw12/r - (vvdw6/r + (C6grid * exponent * beta^6)/r) */
551 fvdw = _mm_mul_pd(_mm_add_pd(vvdw12,_mm_msub_pd(_mm_mul_pd(c6grid_00,one_sixth),_mm_mul_pd(exponent,ewclj6),vvdw6)),rinvsq00);
553 cutoff_mask = _mm_cmplt_pd(rsq00,rcutoff2);
555 /* Update potential sum for this i atom from the interaction with this j atom. */
556 vvdw = _mm_and_pd(vvdw,cutoff_mask);
557 vvdw = _mm_unpacklo_pd(vvdw,_mm_setzero_pd());
558 vvdwsum = _mm_add_pd(vvdwsum,vvdw);
562 fscal = _mm_and_pd(fscal,cutoff_mask);
564 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
566 /* Update vectorial force */
567 fix0 = _mm_macc_pd(dx00,fscal,fix0);
568 fiy0 = _mm_macc_pd(dy00,fscal,fiy0);
569 fiz0 = _mm_macc_pd(dz00,fscal,fiz0);
571 fjx0 = _mm_macc_pd(dx00,fscal,fjx0);
572 fjy0 = _mm_macc_pd(dy00,fscal,fjy0);
573 fjz0 = _mm_macc_pd(dz00,fscal,fjz0);
577 /**************************
578 * CALCULATE INTERACTIONS *
579 **************************/
581 if (gmx_mm_any_lt(rsq10,rcutoff2))
584 r10 = _mm_mul_pd(rsq10,rinv10);
586 /* Compute parameters for interactions between i and j atoms */
587 qq10 = _mm_mul_pd(iq1,jq0);
589 /* EWALD ELECTROSTATICS */
591 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
592 ewrt = _mm_mul_pd(r10,ewtabscale);
593 ewitab = _mm_cvttpd_epi32(ewrt);
595 eweps = _mm_frcz_pd(ewrt);
597 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
599 twoeweps = _mm_add_pd(eweps,eweps);
600 ewitab = _mm_slli_epi32(ewitab,2);
601 ewtabF = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
602 ewtabD = _mm_setzero_pd();
603 GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
604 ewtabV = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
605 ewtabFn = _mm_setzero_pd();
606 GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
607 felec = _mm_macc_pd(eweps,ewtabD,ewtabF);
608 velec = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
609 velec = _mm_mul_pd(qq10,_mm_sub_pd(_mm_sub_pd(rinv10,sh_ewald),velec));
610 felec = _mm_mul_pd(_mm_mul_pd(qq10,rinv10),_mm_sub_pd(rinvsq10,felec));
612 cutoff_mask = _mm_cmplt_pd(rsq10,rcutoff2);
614 /* Update potential sum for this i atom from the interaction with this j atom. */
615 velec = _mm_and_pd(velec,cutoff_mask);
616 velec = _mm_unpacklo_pd(velec,_mm_setzero_pd());
617 velecsum = _mm_add_pd(velecsum,velec);
621 fscal = _mm_and_pd(fscal,cutoff_mask);
623 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
625 /* Update vectorial force */
626 fix1 = _mm_macc_pd(dx10,fscal,fix1);
627 fiy1 = _mm_macc_pd(dy10,fscal,fiy1);
628 fiz1 = _mm_macc_pd(dz10,fscal,fiz1);
630 fjx0 = _mm_macc_pd(dx10,fscal,fjx0);
631 fjy0 = _mm_macc_pd(dy10,fscal,fjy0);
632 fjz0 = _mm_macc_pd(dz10,fscal,fjz0);
636 /**************************
637 * CALCULATE INTERACTIONS *
638 **************************/
640 if (gmx_mm_any_lt(rsq20,rcutoff2))
643 r20 = _mm_mul_pd(rsq20,rinv20);
645 /* Compute parameters for interactions between i and j atoms */
646 qq20 = _mm_mul_pd(iq2,jq0);
648 /* EWALD ELECTROSTATICS */
650 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
651 ewrt = _mm_mul_pd(r20,ewtabscale);
652 ewitab = _mm_cvttpd_epi32(ewrt);
654 eweps = _mm_frcz_pd(ewrt);
656 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
658 twoeweps = _mm_add_pd(eweps,eweps);
659 ewitab = _mm_slli_epi32(ewitab,2);
660 ewtabF = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
661 ewtabD = _mm_setzero_pd();
662 GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
663 ewtabV = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
664 ewtabFn = _mm_setzero_pd();
665 GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
666 felec = _mm_macc_pd(eweps,ewtabD,ewtabF);
667 velec = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
668 velec = _mm_mul_pd(qq20,_mm_sub_pd(_mm_sub_pd(rinv20,sh_ewald),velec));
669 felec = _mm_mul_pd(_mm_mul_pd(qq20,rinv20),_mm_sub_pd(rinvsq20,felec));
671 cutoff_mask = _mm_cmplt_pd(rsq20,rcutoff2);
673 /* Update potential sum for this i atom from the interaction with this j atom. */
674 velec = _mm_and_pd(velec,cutoff_mask);
675 velec = _mm_unpacklo_pd(velec,_mm_setzero_pd());
676 velecsum = _mm_add_pd(velecsum,velec);
680 fscal = _mm_and_pd(fscal,cutoff_mask);
682 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
684 /* Update vectorial force */
685 fix2 = _mm_macc_pd(dx20,fscal,fix2);
686 fiy2 = _mm_macc_pd(dy20,fscal,fiy2);
687 fiz2 = _mm_macc_pd(dz20,fscal,fiz2);
689 fjx0 = _mm_macc_pd(dx20,fscal,fjx0);
690 fjy0 = _mm_macc_pd(dy20,fscal,fjy0);
691 fjz0 = _mm_macc_pd(dz20,fscal,fjz0);
695 /**************************
696 * CALCULATE INTERACTIONS *
697 **************************/
699 if (gmx_mm_any_lt(rsq30,rcutoff2))
702 r30 = _mm_mul_pd(rsq30,rinv30);
704 /* Compute parameters for interactions between i and j atoms */
705 qq30 = _mm_mul_pd(iq3,jq0);
707 /* EWALD ELECTROSTATICS */
709 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
710 ewrt = _mm_mul_pd(r30,ewtabscale);
711 ewitab = _mm_cvttpd_epi32(ewrt);
713 eweps = _mm_frcz_pd(ewrt);
715 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
717 twoeweps = _mm_add_pd(eweps,eweps);
718 ewitab = _mm_slli_epi32(ewitab,2);
719 ewtabF = _mm_load_pd( ewtab + _mm_extract_epi32(ewitab,0) );
720 ewtabD = _mm_setzero_pd();
721 GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
722 ewtabV = _mm_load_sd( ewtab + _mm_extract_epi32(ewitab,0) +2);
723 ewtabFn = _mm_setzero_pd();
724 GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
725 felec = _mm_macc_pd(eweps,ewtabD,ewtabF);
726 velec = _mm_nmacc_pd(_mm_mul_pd(ewtabhalfspace,eweps) ,_mm_add_pd(ewtabF,felec), ewtabV);
727 velec = _mm_mul_pd(qq30,_mm_sub_pd(_mm_sub_pd(rinv30,sh_ewald),velec));
728 felec = _mm_mul_pd(_mm_mul_pd(qq30,rinv30),_mm_sub_pd(rinvsq30,felec));
730 cutoff_mask = _mm_cmplt_pd(rsq30,rcutoff2);
732 /* Update potential sum for this i atom from the interaction with this j atom. */
733 velec = _mm_and_pd(velec,cutoff_mask);
734 velec = _mm_unpacklo_pd(velec,_mm_setzero_pd());
735 velecsum = _mm_add_pd(velecsum,velec);
739 fscal = _mm_and_pd(fscal,cutoff_mask);
741 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
743 /* Update vectorial force */
744 fix3 = _mm_macc_pd(dx30,fscal,fix3);
745 fiy3 = _mm_macc_pd(dy30,fscal,fiy3);
746 fiz3 = _mm_macc_pd(dz30,fscal,fiz3);
748 fjx0 = _mm_macc_pd(dx30,fscal,fjx0);
749 fjy0 = _mm_macc_pd(dy30,fscal,fjy0);
750 fjz0 = _mm_macc_pd(dz30,fscal,fjz0);
754 gmx_mm_decrement_1rvec_1ptr_swizzle_pd(f+j_coord_offsetA,fjx0,fjy0,fjz0);
756 /* Inner loop uses 208 flops */
759 /* End of innermost loop */
761 gmx_mm_update_iforce_4atom_swizzle_pd(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
762 f+i_coord_offset,fshift+i_shift_offset);
765 /* Update potential energies */
766 gmx_mm_update_1pot_pd(velecsum,kernel_data->energygrp_elec+ggid);
767 gmx_mm_update_1pot_pd(vvdwsum,kernel_data->energygrp_vdw+ggid);
769 /* Increment number of inner iterations */
770 inneriter += j_index_end - j_index_start;
772 /* Outer loop uses 26 flops */
775 /* Increment number of outer iterations */
778 /* Update outer/inner flops */
780 inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4_VF,outeriter*26 + inneriter*208);
783 * Gromacs nonbonded kernel: nb_kernel_ElecEwSh_VdwLJEwSh_GeomW4P1_F_avx_128_fma_double
784 * Electrostatics interaction: Ewald
785 * VdW interaction: LJEwald
786 * Geometry: Water4-Particle
787 * Calculate force/pot: Force
790 nb_kernel_ElecEwSh_VdwLJEwSh_GeomW4P1_F_avx_128_fma_double
791 (t_nblist * gmx_restrict nlist,
792 rvec * gmx_restrict xx,
793 rvec * gmx_restrict ff,
794 struct t_forcerec * gmx_restrict fr,
795 t_mdatoms * gmx_restrict mdatoms,
796 nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
797 t_nrnb * gmx_restrict nrnb)
799 /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
800 * just 0 for non-waters.
801 * Suffixes A,B refer to j loop unrolling done with SSE double precision, e.g. for the two different
802 * jnr indices corresponding to data put in the four positions in the SIMD register.
804 int i_shift_offset,i_coord_offset,outeriter,inneriter;
805 int j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
807 int j_coord_offsetA,j_coord_offsetB;
808 int *iinr,*jindex,*jjnr,*shiftidx,*gid;
810 real *shiftvec,*fshift,*x,*f;
811 __m128d tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
813 __m128d ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
815 __m128d ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
817 __m128d ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
819 __m128d ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
820 int vdwjidx0A,vdwjidx0B;
821 __m128d jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
822 __m128d dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
823 __m128d dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
824 __m128d dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
825 __m128d dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
826 __m128d velec,felec,velecsum,facel,crf,krf,krf2;
829 __m128d rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
832 __m128d one_sixth = _mm_set1_pd(1.0/6.0);
833 __m128d one_twelfth = _mm_set1_pd(1.0/12.0);
839 __m128d ewclj,ewclj2,ewclj6,ewcljrsq,poly,exponent,f6A,f6B,sh_lj_ewald;
840 __m128d one_half = _mm_set1_pd(0.5);
841 __m128d minus_one = _mm_set1_pd(-1.0);
843 __m128d ewtabscale,eweps,twoeweps,sh_ewald,ewrt,ewtabhalfspace,ewtabF,ewtabFn,ewtabD,ewtabV;
845 __m128d dummy_mask,cutoff_mask;
846 __m128d signbit = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
847 __m128d one = _mm_set1_pd(1.0);
848 __m128d two = _mm_set1_pd(2.0);
854 jindex = nlist->jindex;
856 shiftidx = nlist->shift;
858 shiftvec = fr->shift_vec[0];
859 fshift = fr->fshift[0];
860 facel = _mm_set1_pd(fr->ic->epsfac);
861 charge = mdatoms->chargeA;
862 nvdwtype = fr->ntype;
864 vdwtype = mdatoms->typeA;
865 vdwgridparam = fr->ljpme_c6grid;
866 sh_lj_ewald = _mm_set1_pd(fr->ic->sh_lj_ewald);
867 ewclj = _mm_set1_pd(fr->ic->ewaldcoeff_lj);
868 ewclj2 = _mm_mul_pd(minus_one,_mm_mul_pd(ewclj,ewclj));
870 sh_ewald = _mm_set1_pd(fr->ic->sh_ewald);
871 ewtab = fr->ic->tabq_coul_F;
872 ewtabscale = _mm_set1_pd(fr->ic->tabq_scale);
873 ewtabhalfspace = _mm_set1_pd(0.5/fr->ic->tabq_scale);
875 /* Setup water-specific parameters */
876 inr = nlist->iinr[0];
877 iq1 = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+1]));
878 iq2 = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+2]));
879 iq3 = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+3]));
880 vdwioffset0 = 2*nvdwtype*vdwtype[inr+0];
882 /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
883 rcutoff_scalar = fr->ic->rcoulomb;
884 rcutoff = _mm_set1_pd(rcutoff_scalar);
885 rcutoff2 = _mm_mul_pd(rcutoff,rcutoff);
887 sh_vdw_invrcut6 = _mm_set1_pd(fr->ic->sh_invrc6);
888 rvdw = _mm_set1_pd(fr->ic->rvdw);
890 /* Avoid stupid compiler warnings */
898 /* Start outer loop over neighborlists */
899 for(iidx=0; iidx<nri; iidx++)
901 /* Load shift vector for this list */
902 i_shift_offset = DIM*shiftidx[iidx];
904 /* Load limits for loop over neighbors */
905 j_index_start = jindex[iidx];
906 j_index_end = jindex[iidx+1];
908 /* Get outer coordinate index */
910 i_coord_offset = DIM*inr;
912 /* Load i particle coords and add shift vector */
913 gmx_mm_load_shift_and_4rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,
914 &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
916 fix0 = _mm_setzero_pd();
917 fiy0 = _mm_setzero_pd();
918 fiz0 = _mm_setzero_pd();
919 fix1 = _mm_setzero_pd();
920 fiy1 = _mm_setzero_pd();
921 fiz1 = _mm_setzero_pd();
922 fix2 = _mm_setzero_pd();
923 fiy2 = _mm_setzero_pd();
924 fiz2 = _mm_setzero_pd();
925 fix3 = _mm_setzero_pd();
926 fiy3 = _mm_setzero_pd();
927 fiz3 = _mm_setzero_pd();
929 /* Start inner kernel loop */
930 for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
933 /* Get j neighbor index, and coordinate index */
936 j_coord_offsetA = DIM*jnrA;
937 j_coord_offsetB = DIM*jnrB;
939 /* load j atom coordinates */
940 gmx_mm_load_1rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
943 /* Calculate displacement vector */
944 dx00 = _mm_sub_pd(ix0,jx0);
945 dy00 = _mm_sub_pd(iy0,jy0);
946 dz00 = _mm_sub_pd(iz0,jz0);
947 dx10 = _mm_sub_pd(ix1,jx0);
948 dy10 = _mm_sub_pd(iy1,jy0);
949 dz10 = _mm_sub_pd(iz1,jz0);
950 dx20 = _mm_sub_pd(ix2,jx0);
951 dy20 = _mm_sub_pd(iy2,jy0);
952 dz20 = _mm_sub_pd(iz2,jz0);
953 dx30 = _mm_sub_pd(ix3,jx0);
954 dy30 = _mm_sub_pd(iy3,jy0);
955 dz30 = _mm_sub_pd(iz3,jz0);
957 /* Calculate squared distance and things based on it */
958 rsq00 = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
959 rsq10 = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
960 rsq20 = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
961 rsq30 = gmx_mm_calc_rsq_pd(dx30,dy30,dz30);
963 rinv00 = avx128fma_invsqrt_d(rsq00);
964 rinv10 = avx128fma_invsqrt_d(rsq10);
965 rinv20 = avx128fma_invsqrt_d(rsq20);
966 rinv30 = avx128fma_invsqrt_d(rsq30);
968 rinvsq00 = _mm_mul_pd(rinv00,rinv00);
969 rinvsq10 = _mm_mul_pd(rinv10,rinv10);
970 rinvsq20 = _mm_mul_pd(rinv20,rinv20);
971 rinvsq30 = _mm_mul_pd(rinv30,rinv30);
973 /* Load parameters for j particles */
974 jq0 = gmx_mm_load_2real_swizzle_pd(charge+jnrA+0,charge+jnrB+0);
975 vdwjidx0A = 2*vdwtype[jnrA+0];
976 vdwjidx0B = 2*vdwtype[jnrB+0];
978 fjx0 = _mm_setzero_pd();
979 fjy0 = _mm_setzero_pd();
980 fjz0 = _mm_setzero_pd();
982 /**************************
983 * CALCULATE INTERACTIONS *
984 **************************/
986 if (gmx_mm_any_lt(rsq00,rcutoff2))
989 r00 = _mm_mul_pd(rsq00,rinv00);
991 /* Compute parameters for interactions between i and j atoms */
992 gmx_mm_load_2pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,
993 vdwparam+vdwioffset0+vdwjidx0B,&c6_00,&c12_00);
994 c6grid_00 = gmx_mm_load_2real_swizzle_pd(vdwgridparam+vdwioffset0+vdwjidx0A,
995 vdwgridparam+vdwioffset0+vdwjidx0B);
997 /* Analytical LJ-PME */
998 rinvsix = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
999 ewcljrsq = _mm_mul_pd(ewclj2,rsq00);
1000 ewclj6 = _mm_mul_pd(ewclj2,_mm_mul_pd(ewclj2,ewclj2));
1001 exponent = avx128fma_exp_d(ewcljrsq);
1002 /* poly = exp(-(beta*r)^2) * (1 + (beta*r)^2 + (beta*r)^4 /2) */
1003 poly = _mm_mul_pd(exponent,_mm_macc_pd(_mm_mul_pd(ewcljrsq,ewcljrsq),one_half,_mm_sub_pd(one,ewcljrsq)));
1004 /* f6A = 6 * C6grid * (1 - poly) */
1005 f6A = _mm_mul_pd(c6grid_00,_mm_sub_pd(one,poly));
1006 /* f6B = C6grid * exponent * beta^6 */
1007 f6B = _mm_mul_pd(_mm_mul_pd(c6grid_00,one_sixth),_mm_mul_pd(exponent,ewclj6));
1008 /* fvdw = 12*C12/r13 - ((6*C6 - f6A)/r6 + f6B)/r */
1009 fvdw = _mm_mul_pd(_mm_macc_pd(_mm_msub_pd(c12_00,rinvsix,_mm_sub_pd(c6_00,f6A)),rinvsix,f6B),rinvsq00);
1011 cutoff_mask = _mm_cmplt_pd(rsq00,rcutoff2);
1015 fscal = _mm_and_pd(fscal,cutoff_mask);
1017 /* Update vectorial force */
1018 fix0 = _mm_macc_pd(dx00,fscal,fix0);
1019 fiy0 = _mm_macc_pd(dy00,fscal,fiy0);
1020 fiz0 = _mm_macc_pd(dz00,fscal,fiz0);
1022 fjx0 = _mm_macc_pd(dx00,fscal,fjx0);
1023 fjy0 = _mm_macc_pd(dy00,fscal,fjy0);
1024 fjz0 = _mm_macc_pd(dz00,fscal,fjz0);
1028 /**************************
1029 * CALCULATE INTERACTIONS *
1030 **************************/
1032 if (gmx_mm_any_lt(rsq10,rcutoff2))
1035 r10 = _mm_mul_pd(rsq10,rinv10);
1037 /* Compute parameters for interactions between i and j atoms */
1038 qq10 = _mm_mul_pd(iq1,jq0);
1040 /* EWALD ELECTROSTATICS */
1042 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1043 ewrt = _mm_mul_pd(r10,ewtabscale);
1044 ewitab = _mm_cvttpd_epi32(ewrt);
1046 eweps = _mm_frcz_pd(ewrt);
1048 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1050 twoeweps = _mm_add_pd(eweps,eweps);
1051 gmx_mm_load_2pair_swizzle_pd(ewtab+_mm_extract_epi32(ewitab,0),ewtab+_mm_extract_epi32(ewitab,1),
1053 felec = _mm_macc_pd(eweps,ewtabFn,_mm_mul_pd( _mm_sub_pd(one,eweps),ewtabF));
1054 felec = _mm_mul_pd(_mm_mul_pd(qq10,rinv10),_mm_sub_pd(rinvsq10,felec));
1056 cutoff_mask = _mm_cmplt_pd(rsq10,rcutoff2);
1060 fscal = _mm_and_pd(fscal,cutoff_mask);
1062 /* Update vectorial force */
1063 fix1 = _mm_macc_pd(dx10,fscal,fix1);
1064 fiy1 = _mm_macc_pd(dy10,fscal,fiy1);
1065 fiz1 = _mm_macc_pd(dz10,fscal,fiz1);
1067 fjx0 = _mm_macc_pd(dx10,fscal,fjx0);
1068 fjy0 = _mm_macc_pd(dy10,fscal,fjy0);
1069 fjz0 = _mm_macc_pd(dz10,fscal,fjz0);
1073 /**************************
1074 * CALCULATE INTERACTIONS *
1075 **************************/
1077 if (gmx_mm_any_lt(rsq20,rcutoff2))
1080 r20 = _mm_mul_pd(rsq20,rinv20);
1082 /* Compute parameters for interactions between i and j atoms */
1083 qq20 = _mm_mul_pd(iq2,jq0);
1085 /* EWALD ELECTROSTATICS */
1087 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1088 ewrt = _mm_mul_pd(r20,ewtabscale);
1089 ewitab = _mm_cvttpd_epi32(ewrt);
1091 eweps = _mm_frcz_pd(ewrt);
1093 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1095 twoeweps = _mm_add_pd(eweps,eweps);
1096 gmx_mm_load_2pair_swizzle_pd(ewtab+_mm_extract_epi32(ewitab,0),ewtab+_mm_extract_epi32(ewitab,1),
1098 felec = _mm_macc_pd(eweps,ewtabFn,_mm_mul_pd( _mm_sub_pd(one,eweps),ewtabF));
1099 felec = _mm_mul_pd(_mm_mul_pd(qq20,rinv20),_mm_sub_pd(rinvsq20,felec));
1101 cutoff_mask = _mm_cmplt_pd(rsq20,rcutoff2);
1105 fscal = _mm_and_pd(fscal,cutoff_mask);
1107 /* Update vectorial force */
1108 fix2 = _mm_macc_pd(dx20,fscal,fix2);
1109 fiy2 = _mm_macc_pd(dy20,fscal,fiy2);
1110 fiz2 = _mm_macc_pd(dz20,fscal,fiz2);
1112 fjx0 = _mm_macc_pd(dx20,fscal,fjx0);
1113 fjy0 = _mm_macc_pd(dy20,fscal,fjy0);
1114 fjz0 = _mm_macc_pd(dz20,fscal,fjz0);
1118 /**************************
1119 * CALCULATE INTERACTIONS *
1120 **************************/
1122 if (gmx_mm_any_lt(rsq30,rcutoff2))
1125 r30 = _mm_mul_pd(rsq30,rinv30);
1127 /* Compute parameters for interactions between i and j atoms */
1128 qq30 = _mm_mul_pd(iq3,jq0);
1130 /* EWALD ELECTROSTATICS */
1132 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1133 ewrt = _mm_mul_pd(r30,ewtabscale);
1134 ewitab = _mm_cvttpd_epi32(ewrt);
1136 eweps = _mm_frcz_pd(ewrt);
1138 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1140 twoeweps = _mm_add_pd(eweps,eweps);
1141 gmx_mm_load_2pair_swizzle_pd(ewtab+_mm_extract_epi32(ewitab,0),ewtab+_mm_extract_epi32(ewitab,1),
1143 felec = _mm_macc_pd(eweps,ewtabFn,_mm_mul_pd( _mm_sub_pd(one,eweps),ewtabF));
1144 felec = _mm_mul_pd(_mm_mul_pd(qq30,rinv30),_mm_sub_pd(rinvsq30,felec));
1146 cutoff_mask = _mm_cmplt_pd(rsq30,rcutoff2);
1150 fscal = _mm_and_pd(fscal,cutoff_mask);
1152 /* Update vectorial force */
1153 fix3 = _mm_macc_pd(dx30,fscal,fix3);
1154 fiy3 = _mm_macc_pd(dy30,fscal,fiy3);
1155 fiz3 = _mm_macc_pd(dz30,fscal,fiz3);
1157 fjx0 = _mm_macc_pd(dx30,fscal,fjx0);
1158 fjy0 = _mm_macc_pd(dy30,fscal,fjy0);
1159 fjz0 = _mm_macc_pd(dz30,fscal,fjz0);
1163 gmx_mm_decrement_1rvec_2ptr_swizzle_pd(f+j_coord_offsetA,f+j_coord_offsetB,fjx0,fjy0,fjz0);
1165 /* Inner loop uses 179 flops */
1168 if(jidx<j_index_end)
1172 j_coord_offsetA = DIM*jnrA;
1174 /* load j atom coordinates */
1175 gmx_mm_load_1rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
1178 /* Calculate displacement vector */
1179 dx00 = _mm_sub_pd(ix0,jx0);
1180 dy00 = _mm_sub_pd(iy0,jy0);
1181 dz00 = _mm_sub_pd(iz0,jz0);
1182 dx10 = _mm_sub_pd(ix1,jx0);
1183 dy10 = _mm_sub_pd(iy1,jy0);
1184 dz10 = _mm_sub_pd(iz1,jz0);
1185 dx20 = _mm_sub_pd(ix2,jx0);
1186 dy20 = _mm_sub_pd(iy2,jy0);
1187 dz20 = _mm_sub_pd(iz2,jz0);
1188 dx30 = _mm_sub_pd(ix3,jx0);
1189 dy30 = _mm_sub_pd(iy3,jy0);
1190 dz30 = _mm_sub_pd(iz3,jz0);
1192 /* Calculate squared distance and things based on it */
1193 rsq00 = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
1194 rsq10 = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
1195 rsq20 = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
1196 rsq30 = gmx_mm_calc_rsq_pd(dx30,dy30,dz30);
1198 rinv00 = avx128fma_invsqrt_d(rsq00);
1199 rinv10 = avx128fma_invsqrt_d(rsq10);
1200 rinv20 = avx128fma_invsqrt_d(rsq20);
1201 rinv30 = avx128fma_invsqrt_d(rsq30);
1203 rinvsq00 = _mm_mul_pd(rinv00,rinv00);
1204 rinvsq10 = _mm_mul_pd(rinv10,rinv10);
1205 rinvsq20 = _mm_mul_pd(rinv20,rinv20);
1206 rinvsq30 = _mm_mul_pd(rinv30,rinv30);
1208 /* Load parameters for j particles */
1209 jq0 = _mm_load_sd(charge+jnrA+0);
1210 vdwjidx0A = 2*vdwtype[jnrA+0];
1212 fjx0 = _mm_setzero_pd();
1213 fjy0 = _mm_setzero_pd();
1214 fjz0 = _mm_setzero_pd();
1216 /**************************
1217 * CALCULATE INTERACTIONS *
1218 **************************/
1220 if (gmx_mm_any_lt(rsq00,rcutoff2))
1223 r00 = _mm_mul_pd(rsq00,rinv00);
1225 /* Compute parameters for interactions between i and j atoms */
1226 gmx_mm_load_1pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,&c6_00,&c12_00);
1227 c6grid_00 = gmx_mm_load_1real_pd(vdwgridparam+vdwioffset0+vdwjidx0A);
1229 /* Analytical LJ-PME */
1230 rinvsix = _mm_mul_pd(_mm_mul_pd(rinvsq00,rinvsq00),rinvsq00);
1231 ewcljrsq = _mm_mul_pd(ewclj2,rsq00);
1232 ewclj6 = _mm_mul_pd(ewclj2,_mm_mul_pd(ewclj2,ewclj2));
1233 exponent = avx128fma_exp_d(ewcljrsq);
1234 /* poly = exp(-(beta*r)^2) * (1 + (beta*r)^2 + (beta*r)^4 /2) */
1235 poly = _mm_mul_pd(exponent,_mm_macc_pd(_mm_mul_pd(ewcljrsq,ewcljrsq),one_half,_mm_sub_pd(one,ewcljrsq)));
1236 /* f6A = 6 * C6grid * (1 - poly) */
1237 f6A = _mm_mul_pd(c6grid_00,_mm_sub_pd(one,poly));
1238 /* f6B = C6grid * exponent * beta^6 */
1239 f6B = _mm_mul_pd(_mm_mul_pd(c6grid_00,one_sixth),_mm_mul_pd(exponent,ewclj6));
1240 /* fvdw = 12*C12/r13 - ((6*C6 - f6A)/r6 + f6B)/r */
1241 fvdw = _mm_mul_pd(_mm_macc_pd(_mm_msub_pd(c12_00,rinvsix,_mm_sub_pd(c6_00,f6A)),rinvsix,f6B),rinvsq00);
1243 cutoff_mask = _mm_cmplt_pd(rsq00,rcutoff2);
1247 fscal = _mm_and_pd(fscal,cutoff_mask);
1249 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1251 /* Update vectorial force */
1252 fix0 = _mm_macc_pd(dx00,fscal,fix0);
1253 fiy0 = _mm_macc_pd(dy00,fscal,fiy0);
1254 fiz0 = _mm_macc_pd(dz00,fscal,fiz0);
1256 fjx0 = _mm_macc_pd(dx00,fscal,fjx0);
1257 fjy0 = _mm_macc_pd(dy00,fscal,fjy0);
1258 fjz0 = _mm_macc_pd(dz00,fscal,fjz0);
1262 /**************************
1263 * CALCULATE INTERACTIONS *
1264 **************************/
1266 if (gmx_mm_any_lt(rsq10,rcutoff2))
1269 r10 = _mm_mul_pd(rsq10,rinv10);
1271 /* Compute parameters for interactions between i and j atoms */
1272 qq10 = _mm_mul_pd(iq1,jq0);
1274 /* EWALD ELECTROSTATICS */
1276 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1277 ewrt = _mm_mul_pd(r10,ewtabscale);
1278 ewitab = _mm_cvttpd_epi32(ewrt);
1280 eweps = _mm_frcz_pd(ewrt);
1282 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1284 twoeweps = _mm_add_pd(eweps,eweps);
1285 gmx_mm_load_1pair_swizzle_pd(ewtab+_mm_extract_epi32(ewitab,0),&ewtabF,&ewtabFn);
1286 felec = _mm_macc_pd(eweps,ewtabFn,_mm_mul_pd( _mm_sub_pd(one,eweps),ewtabF));
1287 felec = _mm_mul_pd(_mm_mul_pd(qq10,rinv10),_mm_sub_pd(rinvsq10,felec));
1289 cutoff_mask = _mm_cmplt_pd(rsq10,rcutoff2);
1293 fscal = _mm_and_pd(fscal,cutoff_mask);
1295 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1297 /* Update vectorial force */
1298 fix1 = _mm_macc_pd(dx10,fscal,fix1);
1299 fiy1 = _mm_macc_pd(dy10,fscal,fiy1);
1300 fiz1 = _mm_macc_pd(dz10,fscal,fiz1);
1302 fjx0 = _mm_macc_pd(dx10,fscal,fjx0);
1303 fjy0 = _mm_macc_pd(dy10,fscal,fjy0);
1304 fjz0 = _mm_macc_pd(dz10,fscal,fjz0);
1308 /**************************
1309 * CALCULATE INTERACTIONS *
1310 **************************/
1312 if (gmx_mm_any_lt(rsq20,rcutoff2))
1315 r20 = _mm_mul_pd(rsq20,rinv20);
1317 /* Compute parameters for interactions between i and j atoms */
1318 qq20 = _mm_mul_pd(iq2,jq0);
1320 /* EWALD ELECTROSTATICS */
1322 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1323 ewrt = _mm_mul_pd(r20,ewtabscale);
1324 ewitab = _mm_cvttpd_epi32(ewrt);
1326 eweps = _mm_frcz_pd(ewrt);
1328 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1330 twoeweps = _mm_add_pd(eweps,eweps);
1331 gmx_mm_load_1pair_swizzle_pd(ewtab+_mm_extract_epi32(ewitab,0),&ewtabF,&ewtabFn);
1332 felec = _mm_macc_pd(eweps,ewtabFn,_mm_mul_pd( _mm_sub_pd(one,eweps),ewtabF));
1333 felec = _mm_mul_pd(_mm_mul_pd(qq20,rinv20),_mm_sub_pd(rinvsq20,felec));
1335 cutoff_mask = _mm_cmplt_pd(rsq20,rcutoff2);
1339 fscal = _mm_and_pd(fscal,cutoff_mask);
1341 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1343 /* Update vectorial force */
1344 fix2 = _mm_macc_pd(dx20,fscal,fix2);
1345 fiy2 = _mm_macc_pd(dy20,fscal,fiy2);
1346 fiz2 = _mm_macc_pd(dz20,fscal,fiz2);
1348 fjx0 = _mm_macc_pd(dx20,fscal,fjx0);
1349 fjy0 = _mm_macc_pd(dy20,fscal,fjy0);
1350 fjz0 = _mm_macc_pd(dz20,fscal,fjz0);
1354 /**************************
1355 * CALCULATE INTERACTIONS *
1356 **************************/
1358 if (gmx_mm_any_lt(rsq30,rcutoff2))
1361 r30 = _mm_mul_pd(rsq30,rinv30);
1363 /* Compute parameters for interactions between i and j atoms */
1364 qq30 = _mm_mul_pd(iq3,jq0);
1366 /* EWALD ELECTROSTATICS */
1368 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1369 ewrt = _mm_mul_pd(r30,ewtabscale);
1370 ewitab = _mm_cvttpd_epi32(ewrt);
1372 eweps = _mm_frcz_pd(ewrt);
1374 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1376 twoeweps = _mm_add_pd(eweps,eweps);
1377 gmx_mm_load_1pair_swizzle_pd(ewtab+_mm_extract_epi32(ewitab,0),&ewtabF,&ewtabFn);
1378 felec = _mm_macc_pd(eweps,ewtabFn,_mm_mul_pd( _mm_sub_pd(one,eweps),ewtabF));
1379 felec = _mm_mul_pd(_mm_mul_pd(qq30,rinv30),_mm_sub_pd(rinvsq30,felec));
1381 cutoff_mask = _mm_cmplt_pd(rsq30,rcutoff2);
1385 fscal = _mm_and_pd(fscal,cutoff_mask);
1387 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1389 /* Update vectorial force */
1390 fix3 = _mm_macc_pd(dx30,fscal,fix3);
1391 fiy3 = _mm_macc_pd(dy30,fscal,fiy3);
1392 fiz3 = _mm_macc_pd(dz30,fscal,fiz3);
1394 fjx0 = _mm_macc_pd(dx30,fscal,fjx0);
1395 fjy0 = _mm_macc_pd(dy30,fscal,fjy0);
1396 fjz0 = _mm_macc_pd(dz30,fscal,fjz0);
1400 gmx_mm_decrement_1rvec_1ptr_swizzle_pd(f+j_coord_offsetA,fjx0,fjy0,fjz0);
1402 /* Inner loop uses 179 flops */
1405 /* End of innermost loop */
1407 gmx_mm_update_iforce_4atom_swizzle_pd(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
1408 f+i_coord_offset,fshift+i_shift_offset);
1410 /* Increment number of inner iterations */
1411 inneriter += j_index_end - j_index_start;
1413 /* Outer loop uses 24 flops */
1416 /* Increment number of outer iterations */
1419 /* Update outer/inner flops */
1421 inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4_F,outeriter*24 + inneriter*179);