2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013,2014, 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 sse2_double kernel generator.
44 #include "../nb_kernel.h"
45 #include "gromacs/legacyheaders/types/simple.h"
46 #include "gromacs/math/vec.h"
47 #include "gromacs/legacyheaders/nrnb.h"
49 #include "gromacs/simd/math_x86_sse2_double.h"
50 #include "kernelutil_x86_sse2_double.h"
53 * Gromacs nonbonded kernel: nb_kernel_ElecRF_VdwCSTab_GeomW4W4_VF_sse2_double
54 * Electrostatics interaction: ReactionField
55 * VdW interaction: CubicSplineTable
56 * Geometry: Water4-Water4
57 * Calculate force/pot: PotentialAndForce
60 nb_kernel_ElecRF_VdwCSTab_GeomW4W4_VF_sse2_double
61 (t_nblist * gmx_restrict nlist,
62 rvec * gmx_restrict xx,
63 rvec * gmx_restrict ff,
64 t_forcerec * gmx_restrict fr,
65 t_mdatoms * gmx_restrict mdatoms,
66 nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
67 t_nrnb * gmx_restrict nrnb)
69 /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
70 * just 0 for non-waters.
71 * Suffixes A,B refer to j loop unrolling done with SSE double precision, e.g. for the two different
72 * jnr indices corresponding to data put in the four positions in the SIMD register.
74 int i_shift_offset,i_coord_offset,outeriter,inneriter;
75 int j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
77 int j_coord_offsetA,j_coord_offsetB;
78 int *iinr,*jindex,*jjnr,*shiftidx,*gid;
80 real *shiftvec,*fshift,*x,*f;
81 __m128d tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
83 __m128d ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
85 __m128d ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
87 __m128d ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
89 __m128d ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
90 int vdwjidx0A,vdwjidx0B;
91 __m128d jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
92 int vdwjidx1A,vdwjidx1B;
93 __m128d jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
94 int vdwjidx2A,vdwjidx2B;
95 __m128d jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
96 int vdwjidx3A,vdwjidx3B;
97 __m128d jx3,jy3,jz3,fjx3,fjy3,fjz3,jq3,isaj3;
98 __m128d dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
99 __m128d dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11;
100 __m128d dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12;
101 __m128d dx13,dy13,dz13,rsq13,rinv13,rinvsq13,r13,qq13,c6_13,c12_13;
102 __m128d dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21;
103 __m128d dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22;
104 __m128d dx23,dy23,dz23,rsq23,rinv23,rinvsq23,r23,qq23,c6_23,c12_23;
105 __m128d dx31,dy31,dz31,rsq31,rinv31,rinvsq31,r31,qq31,c6_31,c12_31;
106 __m128d dx32,dy32,dz32,rsq32,rinv32,rinvsq32,r32,qq32,c6_32,c12_32;
107 __m128d dx33,dy33,dz33,rsq33,rinv33,rinvsq33,r33,qq33,c6_33,c12_33;
108 __m128d velec,felec,velecsum,facel,crf,krf,krf2;
111 __m128d rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
114 __m128d one_sixth = _mm_set1_pd(1.0/6.0);
115 __m128d one_twelfth = _mm_set1_pd(1.0/12.0);
117 __m128i ifour = _mm_set1_epi32(4);
118 __m128d rt,vfeps,vftabscale,Y,F,G,H,Heps,Fp,VV,FF;
120 __m128d dummy_mask,cutoff_mask;
121 __m128d signbit = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
122 __m128d one = _mm_set1_pd(1.0);
123 __m128d two = _mm_set1_pd(2.0);
129 jindex = nlist->jindex;
131 shiftidx = nlist->shift;
133 shiftvec = fr->shift_vec[0];
134 fshift = fr->fshift[0];
135 facel = _mm_set1_pd(fr->epsfac);
136 charge = mdatoms->chargeA;
137 krf = _mm_set1_pd(fr->ic->k_rf);
138 krf2 = _mm_set1_pd(fr->ic->k_rf*2.0);
139 crf = _mm_set1_pd(fr->ic->c_rf);
140 nvdwtype = fr->ntype;
142 vdwtype = mdatoms->typeA;
144 vftab = kernel_data->table_vdw->data;
145 vftabscale = _mm_set1_pd(kernel_data->table_vdw->scale);
147 /* Setup water-specific parameters */
148 inr = nlist->iinr[0];
149 iq1 = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+1]));
150 iq2 = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+2]));
151 iq3 = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+3]));
152 vdwioffset0 = 2*nvdwtype*vdwtype[inr+0];
154 jq1 = _mm_set1_pd(charge[inr+1]);
155 jq2 = _mm_set1_pd(charge[inr+2]);
156 jq3 = _mm_set1_pd(charge[inr+3]);
157 vdwjidx0A = 2*vdwtype[inr+0];
158 c6_00 = _mm_set1_pd(vdwparam[vdwioffset0+vdwjidx0A]);
159 c12_00 = _mm_set1_pd(vdwparam[vdwioffset0+vdwjidx0A+1]);
160 qq11 = _mm_mul_pd(iq1,jq1);
161 qq12 = _mm_mul_pd(iq1,jq2);
162 qq13 = _mm_mul_pd(iq1,jq3);
163 qq21 = _mm_mul_pd(iq2,jq1);
164 qq22 = _mm_mul_pd(iq2,jq2);
165 qq23 = _mm_mul_pd(iq2,jq3);
166 qq31 = _mm_mul_pd(iq3,jq1);
167 qq32 = _mm_mul_pd(iq3,jq2);
168 qq33 = _mm_mul_pd(iq3,jq3);
170 /* Avoid stupid compiler warnings */
178 /* Start outer loop over neighborlists */
179 for(iidx=0; iidx<nri; iidx++)
181 /* Load shift vector for this list */
182 i_shift_offset = DIM*shiftidx[iidx];
184 /* Load limits for loop over neighbors */
185 j_index_start = jindex[iidx];
186 j_index_end = jindex[iidx+1];
188 /* Get outer coordinate index */
190 i_coord_offset = DIM*inr;
192 /* Load i particle coords and add shift vector */
193 gmx_mm_load_shift_and_4rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,
194 &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
196 fix0 = _mm_setzero_pd();
197 fiy0 = _mm_setzero_pd();
198 fiz0 = _mm_setzero_pd();
199 fix1 = _mm_setzero_pd();
200 fiy1 = _mm_setzero_pd();
201 fiz1 = _mm_setzero_pd();
202 fix2 = _mm_setzero_pd();
203 fiy2 = _mm_setzero_pd();
204 fiz2 = _mm_setzero_pd();
205 fix3 = _mm_setzero_pd();
206 fiy3 = _mm_setzero_pd();
207 fiz3 = _mm_setzero_pd();
209 /* Reset potential sums */
210 velecsum = _mm_setzero_pd();
211 vvdwsum = _mm_setzero_pd();
213 /* Start inner kernel loop */
214 for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
217 /* Get j neighbor index, and coordinate index */
220 j_coord_offsetA = DIM*jnrA;
221 j_coord_offsetB = DIM*jnrB;
223 /* load j atom coordinates */
224 gmx_mm_load_4rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
225 &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,
226 &jy2,&jz2,&jx3,&jy3,&jz3);
228 /* Calculate displacement vector */
229 dx00 = _mm_sub_pd(ix0,jx0);
230 dy00 = _mm_sub_pd(iy0,jy0);
231 dz00 = _mm_sub_pd(iz0,jz0);
232 dx11 = _mm_sub_pd(ix1,jx1);
233 dy11 = _mm_sub_pd(iy1,jy1);
234 dz11 = _mm_sub_pd(iz1,jz1);
235 dx12 = _mm_sub_pd(ix1,jx2);
236 dy12 = _mm_sub_pd(iy1,jy2);
237 dz12 = _mm_sub_pd(iz1,jz2);
238 dx13 = _mm_sub_pd(ix1,jx3);
239 dy13 = _mm_sub_pd(iy1,jy3);
240 dz13 = _mm_sub_pd(iz1,jz3);
241 dx21 = _mm_sub_pd(ix2,jx1);
242 dy21 = _mm_sub_pd(iy2,jy1);
243 dz21 = _mm_sub_pd(iz2,jz1);
244 dx22 = _mm_sub_pd(ix2,jx2);
245 dy22 = _mm_sub_pd(iy2,jy2);
246 dz22 = _mm_sub_pd(iz2,jz2);
247 dx23 = _mm_sub_pd(ix2,jx3);
248 dy23 = _mm_sub_pd(iy2,jy3);
249 dz23 = _mm_sub_pd(iz2,jz3);
250 dx31 = _mm_sub_pd(ix3,jx1);
251 dy31 = _mm_sub_pd(iy3,jy1);
252 dz31 = _mm_sub_pd(iz3,jz1);
253 dx32 = _mm_sub_pd(ix3,jx2);
254 dy32 = _mm_sub_pd(iy3,jy2);
255 dz32 = _mm_sub_pd(iz3,jz2);
256 dx33 = _mm_sub_pd(ix3,jx3);
257 dy33 = _mm_sub_pd(iy3,jy3);
258 dz33 = _mm_sub_pd(iz3,jz3);
260 /* Calculate squared distance and things based on it */
261 rsq00 = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
262 rsq11 = gmx_mm_calc_rsq_pd(dx11,dy11,dz11);
263 rsq12 = gmx_mm_calc_rsq_pd(dx12,dy12,dz12);
264 rsq13 = gmx_mm_calc_rsq_pd(dx13,dy13,dz13);
265 rsq21 = gmx_mm_calc_rsq_pd(dx21,dy21,dz21);
266 rsq22 = gmx_mm_calc_rsq_pd(dx22,dy22,dz22);
267 rsq23 = gmx_mm_calc_rsq_pd(dx23,dy23,dz23);
268 rsq31 = gmx_mm_calc_rsq_pd(dx31,dy31,dz31);
269 rsq32 = gmx_mm_calc_rsq_pd(dx32,dy32,dz32);
270 rsq33 = gmx_mm_calc_rsq_pd(dx33,dy33,dz33);
272 rinv00 = gmx_mm_invsqrt_pd(rsq00);
273 rinv11 = gmx_mm_invsqrt_pd(rsq11);
274 rinv12 = gmx_mm_invsqrt_pd(rsq12);
275 rinv13 = gmx_mm_invsqrt_pd(rsq13);
276 rinv21 = gmx_mm_invsqrt_pd(rsq21);
277 rinv22 = gmx_mm_invsqrt_pd(rsq22);
278 rinv23 = gmx_mm_invsqrt_pd(rsq23);
279 rinv31 = gmx_mm_invsqrt_pd(rsq31);
280 rinv32 = gmx_mm_invsqrt_pd(rsq32);
281 rinv33 = gmx_mm_invsqrt_pd(rsq33);
283 rinvsq11 = _mm_mul_pd(rinv11,rinv11);
284 rinvsq12 = _mm_mul_pd(rinv12,rinv12);
285 rinvsq13 = _mm_mul_pd(rinv13,rinv13);
286 rinvsq21 = _mm_mul_pd(rinv21,rinv21);
287 rinvsq22 = _mm_mul_pd(rinv22,rinv22);
288 rinvsq23 = _mm_mul_pd(rinv23,rinv23);
289 rinvsq31 = _mm_mul_pd(rinv31,rinv31);
290 rinvsq32 = _mm_mul_pd(rinv32,rinv32);
291 rinvsq33 = _mm_mul_pd(rinv33,rinv33);
293 fjx0 = _mm_setzero_pd();
294 fjy0 = _mm_setzero_pd();
295 fjz0 = _mm_setzero_pd();
296 fjx1 = _mm_setzero_pd();
297 fjy1 = _mm_setzero_pd();
298 fjz1 = _mm_setzero_pd();
299 fjx2 = _mm_setzero_pd();
300 fjy2 = _mm_setzero_pd();
301 fjz2 = _mm_setzero_pd();
302 fjx3 = _mm_setzero_pd();
303 fjy3 = _mm_setzero_pd();
304 fjz3 = _mm_setzero_pd();
306 /**************************
307 * CALCULATE INTERACTIONS *
308 **************************/
310 r00 = _mm_mul_pd(rsq00,rinv00);
312 /* Calculate table index by multiplying r with table scale and truncate to integer */
313 rt = _mm_mul_pd(r00,vftabscale);
314 vfitab = _mm_cvttpd_epi32(rt);
315 vfeps = _mm_sub_pd(rt,_mm_cvtepi32_pd(vfitab));
316 vfitab = _mm_slli_epi32(vfitab,3);
318 /* CUBIC SPLINE TABLE DISPERSION */
319 Y = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) );
320 F = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,1) );
321 GMX_MM_TRANSPOSE2_PD(Y,F);
322 G = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) +2);
323 H = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,1) +2);
324 GMX_MM_TRANSPOSE2_PD(G,H);
325 Heps = _mm_mul_pd(vfeps,H);
326 Fp = _mm_add_pd(F,_mm_mul_pd(vfeps,_mm_add_pd(G,Heps)));
327 VV = _mm_add_pd(Y,_mm_mul_pd(vfeps,Fp));
328 vvdw6 = _mm_mul_pd(c6_00,VV);
329 FF = _mm_add_pd(Fp,_mm_mul_pd(vfeps,_mm_add_pd(G,_mm_add_pd(Heps,Heps))));
330 fvdw6 = _mm_mul_pd(c6_00,FF);
332 /* CUBIC SPLINE TABLE REPULSION */
333 vfitab = _mm_add_epi32(vfitab,ifour);
334 Y = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) );
335 F = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,1) );
336 GMX_MM_TRANSPOSE2_PD(Y,F);
337 G = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) +2);
338 H = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,1) +2);
339 GMX_MM_TRANSPOSE2_PD(G,H);
340 Heps = _mm_mul_pd(vfeps,H);
341 Fp = _mm_add_pd(F,_mm_mul_pd(vfeps,_mm_add_pd(G,Heps)));
342 VV = _mm_add_pd(Y,_mm_mul_pd(vfeps,Fp));
343 vvdw12 = _mm_mul_pd(c12_00,VV);
344 FF = _mm_add_pd(Fp,_mm_mul_pd(vfeps,_mm_add_pd(G,_mm_add_pd(Heps,Heps))));
345 fvdw12 = _mm_mul_pd(c12_00,FF);
346 vvdw = _mm_add_pd(vvdw12,vvdw6);
347 fvdw = _mm_xor_pd(signbit,_mm_mul_pd(_mm_add_pd(fvdw6,fvdw12),_mm_mul_pd(vftabscale,rinv00)));
349 /* Update potential sum for this i atom from the interaction with this j atom. */
350 vvdwsum = _mm_add_pd(vvdwsum,vvdw);
354 /* Calculate temporary vectorial force */
355 tx = _mm_mul_pd(fscal,dx00);
356 ty = _mm_mul_pd(fscal,dy00);
357 tz = _mm_mul_pd(fscal,dz00);
359 /* Update vectorial force */
360 fix0 = _mm_add_pd(fix0,tx);
361 fiy0 = _mm_add_pd(fiy0,ty);
362 fiz0 = _mm_add_pd(fiz0,tz);
364 fjx0 = _mm_add_pd(fjx0,tx);
365 fjy0 = _mm_add_pd(fjy0,ty);
366 fjz0 = _mm_add_pd(fjz0,tz);
368 /**************************
369 * CALCULATE INTERACTIONS *
370 **************************/
372 /* REACTION-FIELD ELECTROSTATICS */
373 velec = _mm_mul_pd(qq11,_mm_sub_pd(_mm_add_pd(rinv11,_mm_mul_pd(krf,rsq11)),crf));
374 felec = _mm_mul_pd(qq11,_mm_sub_pd(_mm_mul_pd(rinv11,rinvsq11),krf2));
376 /* Update potential sum for this i atom from the interaction with this j atom. */
377 velecsum = _mm_add_pd(velecsum,velec);
381 /* Calculate temporary vectorial force */
382 tx = _mm_mul_pd(fscal,dx11);
383 ty = _mm_mul_pd(fscal,dy11);
384 tz = _mm_mul_pd(fscal,dz11);
386 /* Update vectorial force */
387 fix1 = _mm_add_pd(fix1,tx);
388 fiy1 = _mm_add_pd(fiy1,ty);
389 fiz1 = _mm_add_pd(fiz1,tz);
391 fjx1 = _mm_add_pd(fjx1,tx);
392 fjy1 = _mm_add_pd(fjy1,ty);
393 fjz1 = _mm_add_pd(fjz1,tz);
395 /**************************
396 * CALCULATE INTERACTIONS *
397 **************************/
399 /* REACTION-FIELD ELECTROSTATICS */
400 velec = _mm_mul_pd(qq12,_mm_sub_pd(_mm_add_pd(rinv12,_mm_mul_pd(krf,rsq12)),crf));
401 felec = _mm_mul_pd(qq12,_mm_sub_pd(_mm_mul_pd(rinv12,rinvsq12),krf2));
403 /* Update potential sum for this i atom from the interaction with this j atom. */
404 velecsum = _mm_add_pd(velecsum,velec);
408 /* Calculate temporary vectorial force */
409 tx = _mm_mul_pd(fscal,dx12);
410 ty = _mm_mul_pd(fscal,dy12);
411 tz = _mm_mul_pd(fscal,dz12);
413 /* Update vectorial force */
414 fix1 = _mm_add_pd(fix1,tx);
415 fiy1 = _mm_add_pd(fiy1,ty);
416 fiz1 = _mm_add_pd(fiz1,tz);
418 fjx2 = _mm_add_pd(fjx2,tx);
419 fjy2 = _mm_add_pd(fjy2,ty);
420 fjz2 = _mm_add_pd(fjz2,tz);
422 /**************************
423 * CALCULATE INTERACTIONS *
424 **************************/
426 /* REACTION-FIELD ELECTROSTATICS */
427 velec = _mm_mul_pd(qq13,_mm_sub_pd(_mm_add_pd(rinv13,_mm_mul_pd(krf,rsq13)),crf));
428 felec = _mm_mul_pd(qq13,_mm_sub_pd(_mm_mul_pd(rinv13,rinvsq13),krf2));
430 /* Update potential sum for this i atom from the interaction with this j atom. */
431 velecsum = _mm_add_pd(velecsum,velec);
435 /* Calculate temporary vectorial force */
436 tx = _mm_mul_pd(fscal,dx13);
437 ty = _mm_mul_pd(fscal,dy13);
438 tz = _mm_mul_pd(fscal,dz13);
440 /* Update vectorial force */
441 fix1 = _mm_add_pd(fix1,tx);
442 fiy1 = _mm_add_pd(fiy1,ty);
443 fiz1 = _mm_add_pd(fiz1,tz);
445 fjx3 = _mm_add_pd(fjx3,tx);
446 fjy3 = _mm_add_pd(fjy3,ty);
447 fjz3 = _mm_add_pd(fjz3,tz);
449 /**************************
450 * CALCULATE INTERACTIONS *
451 **************************/
453 /* REACTION-FIELD ELECTROSTATICS */
454 velec = _mm_mul_pd(qq21,_mm_sub_pd(_mm_add_pd(rinv21,_mm_mul_pd(krf,rsq21)),crf));
455 felec = _mm_mul_pd(qq21,_mm_sub_pd(_mm_mul_pd(rinv21,rinvsq21),krf2));
457 /* Update potential sum for this i atom from the interaction with this j atom. */
458 velecsum = _mm_add_pd(velecsum,velec);
462 /* Calculate temporary vectorial force */
463 tx = _mm_mul_pd(fscal,dx21);
464 ty = _mm_mul_pd(fscal,dy21);
465 tz = _mm_mul_pd(fscal,dz21);
467 /* Update vectorial force */
468 fix2 = _mm_add_pd(fix2,tx);
469 fiy2 = _mm_add_pd(fiy2,ty);
470 fiz2 = _mm_add_pd(fiz2,tz);
472 fjx1 = _mm_add_pd(fjx1,tx);
473 fjy1 = _mm_add_pd(fjy1,ty);
474 fjz1 = _mm_add_pd(fjz1,tz);
476 /**************************
477 * CALCULATE INTERACTIONS *
478 **************************/
480 /* REACTION-FIELD ELECTROSTATICS */
481 velec = _mm_mul_pd(qq22,_mm_sub_pd(_mm_add_pd(rinv22,_mm_mul_pd(krf,rsq22)),crf));
482 felec = _mm_mul_pd(qq22,_mm_sub_pd(_mm_mul_pd(rinv22,rinvsq22),krf2));
484 /* Update potential sum for this i atom from the interaction with this j atom. */
485 velecsum = _mm_add_pd(velecsum,velec);
489 /* Calculate temporary vectorial force */
490 tx = _mm_mul_pd(fscal,dx22);
491 ty = _mm_mul_pd(fscal,dy22);
492 tz = _mm_mul_pd(fscal,dz22);
494 /* Update vectorial force */
495 fix2 = _mm_add_pd(fix2,tx);
496 fiy2 = _mm_add_pd(fiy2,ty);
497 fiz2 = _mm_add_pd(fiz2,tz);
499 fjx2 = _mm_add_pd(fjx2,tx);
500 fjy2 = _mm_add_pd(fjy2,ty);
501 fjz2 = _mm_add_pd(fjz2,tz);
503 /**************************
504 * CALCULATE INTERACTIONS *
505 **************************/
507 /* REACTION-FIELD ELECTROSTATICS */
508 velec = _mm_mul_pd(qq23,_mm_sub_pd(_mm_add_pd(rinv23,_mm_mul_pd(krf,rsq23)),crf));
509 felec = _mm_mul_pd(qq23,_mm_sub_pd(_mm_mul_pd(rinv23,rinvsq23),krf2));
511 /* Update potential sum for this i atom from the interaction with this j atom. */
512 velecsum = _mm_add_pd(velecsum,velec);
516 /* Calculate temporary vectorial force */
517 tx = _mm_mul_pd(fscal,dx23);
518 ty = _mm_mul_pd(fscal,dy23);
519 tz = _mm_mul_pd(fscal,dz23);
521 /* Update vectorial force */
522 fix2 = _mm_add_pd(fix2,tx);
523 fiy2 = _mm_add_pd(fiy2,ty);
524 fiz2 = _mm_add_pd(fiz2,tz);
526 fjx3 = _mm_add_pd(fjx3,tx);
527 fjy3 = _mm_add_pd(fjy3,ty);
528 fjz3 = _mm_add_pd(fjz3,tz);
530 /**************************
531 * CALCULATE INTERACTIONS *
532 **************************/
534 /* REACTION-FIELD ELECTROSTATICS */
535 velec = _mm_mul_pd(qq31,_mm_sub_pd(_mm_add_pd(rinv31,_mm_mul_pd(krf,rsq31)),crf));
536 felec = _mm_mul_pd(qq31,_mm_sub_pd(_mm_mul_pd(rinv31,rinvsq31),krf2));
538 /* Update potential sum for this i atom from the interaction with this j atom. */
539 velecsum = _mm_add_pd(velecsum,velec);
543 /* Calculate temporary vectorial force */
544 tx = _mm_mul_pd(fscal,dx31);
545 ty = _mm_mul_pd(fscal,dy31);
546 tz = _mm_mul_pd(fscal,dz31);
548 /* Update vectorial force */
549 fix3 = _mm_add_pd(fix3,tx);
550 fiy3 = _mm_add_pd(fiy3,ty);
551 fiz3 = _mm_add_pd(fiz3,tz);
553 fjx1 = _mm_add_pd(fjx1,tx);
554 fjy1 = _mm_add_pd(fjy1,ty);
555 fjz1 = _mm_add_pd(fjz1,tz);
557 /**************************
558 * CALCULATE INTERACTIONS *
559 **************************/
561 /* REACTION-FIELD ELECTROSTATICS */
562 velec = _mm_mul_pd(qq32,_mm_sub_pd(_mm_add_pd(rinv32,_mm_mul_pd(krf,rsq32)),crf));
563 felec = _mm_mul_pd(qq32,_mm_sub_pd(_mm_mul_pd(rinv32,rinvsq32),krf2));
565 /* Update potential sum for this i atom from the interaction with this j atom. */
566 velecsum = _mm_add_pd(velecsum,velec);
570 /* Calculate temporary vectorial force */
571 tx = _mm_mul_pd(fscal,dx32);
572 ty = _mm_mul_pd(fscal,dy32);
573 tz = _mm_mul_pd(fscal,dz32);
575 /* Update vectorial force */
576 fix3 = _mm_add_pd(fix3,tx);
577 fiy3 = _mm_add_pd(fiy3,ty);
578 fiz3 = _mm_add_pd(fiz3,tz);
580 fjx2 = _mm_add_pd(fjx2,tx);
581 fjy2 = _mm_add_pd(fjy2,ty);
582 fjz2 = _mm_add_pd(fjz2,tz);
584 /**************************
585 * CALCULATE INTERACTIONS *
586 **************************/
588 /* REACTION-FIELD ELECTROSTATICS */
589 velec = _mm_mul_pd(qq33,_mm_sub_pd(_mm_add_pd(rinv33,_mm_mul_pd(krf,rsq33)),crf));
590 felec = _mm_mul_pd(qq33,_mm_sub_pd(_mm_mul_pd(rinv33,rinvsq33),krf2));
592 /* Update potential sum for this i atom from the interaction with this j atom. */
593 velecsum = _mm_add_pd(velecsum,velec);
597 /* Calculate temporary vectorial force */
598 tx = _mm_mul_pd(fscal,dx33);
599 ty = _mm_mul_pd(fscal,dy33);
600 tz = _mm_mul_pd(fscal,dz33);
602 /* Update vectorial force */
603 fix3 = _mm_add_pd(fix3,tx);
604 fiy3 = _mm_add_pd(fiy3,ty);
605 fiz3 = _mm_add_pd(fiz3,tz);
607 fjx3 = _mm_add_pd(fjx3,tx);
608 fjy3 = _mm_add_pd(fjy3,ty);
609 fjz3 = _mm_add_pd(fjz3,tz);
611 gmx_mm_decrement_4rvec_2ptr_swizzle_pd(f+j_coord_offsetA,f+j_coord_offsetB,fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,fjx2,fjy2,fjz2,fjx3,fjy3,fjz3);
613 /* Inner loop uses 347 flops */
620 j_coord_offsetA = DIM*jnrA;
622 /* load j atom coordinates */
623 gmx_mm_load_4rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
624 &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,
625 &jy2,&jz2,&jx3,&jy3,&jz3);
627 /* Calculate displacement vector */
628 dx00 = _mm_sub_pd(ix0,jx0);
629 dy00 = _mm_sub_pd(iy0,jy0);
630 dz00 = _mm_sub_pd(iz0,jz0);
631 dx11 = _mm_sub_pd(ix1,jx1);
632 dy11 = _mm_sub_pd(iy1,jy1);
633 dz11 = _mm_sub_pd(iz1,jz1);
634 dx12 = _mm_sub_pd(ix1,jx2);
635 dy12 = _mm_sub_pd(iy1,jy2);
636 dz12 = _mm_sub_pd(iz1,jz2);
637 dx13 = _mm_sub_pd(ix1,jx3);
638 dy13 = _mm_sub_pd(iy1,jy3);
639 dz13 = _mm_sub_pd(iz1,jz3);
640 dx21 = _mm_sub_pd(ix2,jx1);
641 dy21 = _mm_sub_pd(iy2,jy1);
642 dz21 = _mm_sub_pd(iz2,jz1);
643 dx22 = _mm_sub_pd(ix2,jx2);
644 dy22 = _mm_sub_pd(iy2,jy2);
645 dz22 = _mm_sub_pd(iz2,jz2);
646 dx23 = _mm_sub_pd(ix2,jx3);
647 dy23 = _mm_sub_pd(iy2,jy3);
648 dz23 = _mm_sub_pd(iz2,jz3);
649 dx31 = _mm_sub_pd(ix3,jx1);
650 dy31 = _mm_sub_pd(iy3,jy1);
651 dz31 = _mm_sub_pd(iz3,jz1);
652 dx32 = _mm_sub_pd(ix3,jx2);
653 dy32 = _mm_sub_pd(iy3,jy2);
654 dz32 = _mm_sub_pd(iz3,jz2);
655 dx33 = _mm_sub_pd(ix3,jx3);
656 dy33 = _mm_sub_pd(iy3,jy3);
657 dz33 = _mm_sub_pd(iz3,jz3);
659 /* Calculate squared distance and things based on it */
660 rsq00 = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
661 rsq11 = gmx_mm_calc_rsq_pd(dx11,dy11,dz11);
662 rsq12 = gmx_mm_calc_rsq_pd(dx12,dy12,dz12);
663 rsq13 = gmx_mm_calc_rsq_pd(dx13,dy13,dz13);
664 rsq21 = gmx_mm_calc_rsq_pd(dx21,dy21,dz21);
665 rsq22 = gmx_mm_calc_rsq_pd(dx22,dy22,dz22);
666 rsq23 = gmx_mm_calc_rsq_pd(dx23,dy23,dz23);
667 rsq31 = gmx_mm_calc_rsq_pd(dx31,dy31,dz31);
668 rsq32 = gmx_mm_calc_rsq_pd(dx32,dy32,dz32);
669 rsq33 = gmx_mm_calc_rsq_pd(dx33,dy33,dz33);
671 rinv00 = gmx_mm_invsqrt_pd(rsq00);
672 rinv11 = gmx_mm_invsqrt_pd(rsq11);
673 rinv12 = gmx_mm_invsqrt_pd(rsq12);
674 rinv13 = gmx_mm_invsqrt_pd(rsq13);
675 rinv21 = gmx_mm_invsqrt_pd(rsq21);
676 rinv22 = gmx_mm_invsqrt_pd(rsq22);
677 rinv23 = gmx_mm_invsqrt_pd(rsq23);
678 rinv31 = gmx_mm_invsqrt_pd(rsq31);
679 rinv32 = gmx_mm_invsqrt_pd(rsq32);
680 rinv33 = gmx_mm_invsqrt_pd(rsq33);
682 rinvsq11 = _mm_mul_pd(rinv11,rinv11);
683 rinvsq12 = _mm_mul_pd(rinv12,rinv12);
684 rinvsq13 = _mm_mul_pd(rinv13,rinv13);
685 rinvsq21 = _mm_mul_pd(rinv21,rinv21);
686 rinvsq22 = _mm_mul_pd(rinv22,rinv22);
687 rinvsq23 = _mm_mul_pd(rinv23,rinv23);
688 rinvsq31 = _mm_mul_pd(rinv31,rinv31);
689 rinvsq32 = _mm_mul_pd(rinv32,rinv32);
690 rinvsq33 = _mm_mul_pd(rinv33,rinv33);
692 fjx0 = _mm_setzero_pd();
693 fjy0 = _mm_setzero_pd();
694 fjz0 = _mm_setzero_pd();
695 fjx1 = _mm_setzero_pd();
696 fjy1 = _mm_setzero_pd();
697 fjz1 = _mm_setzero_pd();
698 fjx2 = _mm_setzero_pd();
699 fjy2 = _mm_setzero_pd();
700 fjz2 = _mm_setzero_pd();
701 fjx3 = _mm_setzero_pd();
702 fjy3 = _mm_setzero_pd();
703 fjz3 = _mm_setzero_pd();
705 /**************************
706 * CALCULATE INTERACTIONS *
707 **************************/
709 r00 = _mm_mul_pd(rsq00,rinv00);
711 /* Calculate table index by multiplying r with table scale and truncate to integer */
712 rt = _mm_mul_pd(r00,vftabscale);
713 vfitab = _mm_cvttpd_epi32(rt);
714 vfeps = _mm_sub_pd(rt,_mm_cvtepi32_pd(vfitab));
715 vfitab = _mm_slli_epi32(vfitab,3);
717 /* CUBIC SPLINE TABLE DISPERSION */
718 Y = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) );
719 F = _mm_setzero_pd();
720 GMX_MM_TRANSPOSE2_PD(Y,F);
721 G = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) +2);
722 H = _mm_setzero_pd();
723 GMX_MM_TRANSPOSE2_PD(G,H);
724 Heps = _mm_mul_pd(vfeps,H);
725 Fp = _mm_add_pd(F,_mm_mul_pd(vfeps,_mm_add_pd(G,Heps)));
726 VV = _mm_add_pd(Y,_mm_mul_pd(vfeps,Fp));
727 vvdw6 = _mm_mul_pd(c6_00,VV);
728 FF = _mm_add_pd(Fp,_mm_mul_pd(vfeps,_mm_add_pd(G,_mm_add_pd(Heps,Heps))));
729 fvdw6 = _mm_mul_pd(c6_00,FF);
731 /* CUBIC SPLINE TABLE REPULSION */
732 vfitab = _mm_add_epi32(vfitab,ifour);
733 Y = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) );
734 F = _mm_setzero_pd();
735 GMX_MM_TRANSPOSE2_PD(Y,F);
736 G = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) +2);
737 H = _mm_setzero_pd();
738 GMX_MM_TRANSPOSE2_PD(G,H);
739 Heps = _mm_mul_pd(vfeps,H);
740 Fp = _mm_add_pd(F,_mm_mul_pd(vfeps,_mm_add_pd(G,Heps)));
741 VV = _mm_add_pd(Y,_mm_mul_pd(vfeps,Fp));
742 vvdw12 = _mm_mul_pd(c12_00,VV);
743 FF = _mm_add_pd(Fp,_mm_mul_pd(vfeps,_mm_add_pd(G,_mm_add_pd(Heps,Heps))));
744 fvdw12 = _mm_mul_pd(c12_00,FF);
745 vvdw = _mm_add_pd(vvdw12,vvdw6);
746 fvdw = _mm_xor_pd(signbit,_mm_mul_pd(_mm_add_pd(fvdw6,fvdw12),_mm_mul_pd(vftabscale,rinv00)));
748 /* Update potential sum for this i atom from the interaction with this j atom. */
749 vvdw = _mm_unpacklo_pd(vvdw,_mm_setzero_pd());
750 vvdwsum = _mm_add_pd(vvdwsum,vvdw);
754 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
756 /* Calculate temporary vectorial force */
757 tx = _mm_mul_pd(fscal,dx00);
758 ty = _mm_mul_pd(fscal,dy00);
759 tz = _mm_mul_pd(fscal,dz00);
761 /* Update vectorial force */
762 fix0 = _mm_add_pd(fix0,tx);
763 fiy0 = _mm_add_pd(fiy0,ty);
764 fiz0 = _mm_add_pd(fiz0,tz);
766 fjx0 = _mm_add_pd(fjx0,tx);
767 fjy0 = _mm_add_pd(fjy0,ty);
768 fjz0 = _mm_add_pd(fjz0,tz);
770 /**************************
771 * CALCULATE INTERACTIONS *
772 **************************/
774 /* REACTION-FIELD ELECTROSTATICS */
775 velec = _mm_mul_pd(qq11,_mm_sub_pd(_mm_add_pd(rinv11,_mm_mul_pd(krf,rsq11)),crf));
776 felec = _mm_mul_pd(qq11,_mm_sub_pd(_mm_mul_pd(rinv11,rinvsq11),krf2));
778 /* Update potential sum for this i atom from the interaction with this j atom. */
779 velec = _mm_unpacklo_pd(velec,_mm_setzero_pd());
780 velecsum = _mm_add_pd(velecsum,velec);
784 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
786 /* Calculate temporary vectorial force */
787 tx = _mm_mul_pd(fscal,dx11);
788 ty = _mm_mul_pd(fscal,dy11);
789 tz = _mm_mul_pd(fscal,dz11);
791 /* Update vectorial force */
792 fix1 = _mm_add_pd(fix1,tx);
793 fiy1 = _mm_add_pd(fiy1,ty);
794 fiz1 = _mm_add_pd(fiz1,tz);
796 fjx1 = _mm_add_pd(fjx1,tx);
797 fjy1 = _mm_add_pd(fjy1,ty);
798 fjz1 = _mm_add_pd(fjz1,tz);
800 /**************************
801 * CALCULATE INTERACTIONS *
802 **************************/
804 /* REACTION-FIELD ELECTROSTATICS */
805 velec = _mm_mul_pd(qq12,_mm_sub_pd(_mm_add_pd(rinv12,_mm_mul_pd(krf,rsq12)),crf));
806 felec = _mm_mul_pd(qq12,_mm_sub_pd(_mm_mul_pd(rinv12,rinvsq12),krf2));
808 /* Update potential sum for this i atom from the interaction with this j atom. */
809 velec = _mm_unpacklo_pd(velec,_mm_setzero_pd());
810 velecsum = _mm_add_pd(velecsum,velec);
814 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
816 /* Calculate temporary vectorial force */
817 tx = _mm_mul_pd(fscal,dx12);
818 ty = _mm_mul_pd(fscal,dy12);
819 tz = _mm_mul_pd(fscal,dz12);
821 /* Update vectorial force */
822 fix1 = _mm_add_pd(fix1,tx);
823 fiy1 = _mm_add_pd(fiy1,ty);
824 fiz1 = _mm_add_pd(fiz1,tz);
826 fjx2 = _mm_add_pd(fjx2,tx);
827 fjy2 = _mm_add_pd(fjy2,ty);
828 fjz2 = _mm_add_pd(fjz2,tz);
830 /**************************
831 * CALCULATE INTERACTIONS *
832 **************************/
834 /* REACTION-FIELD ELECTROSTATICS */
835 velec = _mm_mul_pd(qq13,_mm_sub_pd(_mm_add_pd(rinv13,_mm_mul_pd(krf,rsq13)),crf));
836 felec = _mm_mul_pd(qq13,_mm_sub_pd(_mm_mul_pd(rinv13,rinvsq13),krf2));
838 /* Update potential sum for this i atom from the interaction with this j atom. */
839 velec = _mm_unpacklo_pd(velec,_mm_setzero_pd());
840 velecsum = _mm_add_pd(velecsum,velec);
844 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
846 /* Calculate temporary vectorial force */
847 tx = _mm_mul_pd(fscal,dx13);
848 ty = _mm_mul_pd(fscal,dy13);
849 tz = _mm_mul_pd(fscal,dz13);
851 /* Update vectorial force */
852 fix1 = _mm_add_pd(fix1,tx);
853 fiy1 = _mm_add_pd(fiy1,ty);
854 fiz1 = _mm_add_pd(fiz1,tz);
856 fjx3 = _mm_add_pd(fjx3,tx);
857 fjy3 = _mm_add_pd(fjy3,ty);
858 fjz3 = _mm_add_pd(fjz3,tz);
860 /**************************
861 * CALCULATE INTERACTIONS *
862 **************************/
864 /* REACTION-FIELD ELECTROSTATICS */
865 velec = _mm_mul_pd(qq21,_mm_sub_pd(_mm_add_pd(rinv21,_mm_mul_pd(krf,rsq21)),crf));
866 felec = _mm_mul_pd(qq21,_mm_sub_pd(_mm_mul_pd(rinv21,rinvsq21),krf2));
868 /* Update potential sum for this i atom from the interaction with this j atom. */
869 velec = _mm_unpacklo_pd(velec,_mm_setzero_pd());
870 velecsum = _mm_add_pd(velecsum,velec);
874 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
876 /* Calculate temporary vectorial force */
877 tx = _mm_mul_pd(fscal,dx21);
878 ty = _mm_mul_pd(fscal,dy21);
879 tz = _mm_mul_pd(fscal,dz21);
881 /* Update vectorial force */
882 fix2 = _mm_add_pd(fix2,tx);
883 fiy2 = _mm_add_pd(fiy2,ty);
884 fiz2 = _mm_add_pd(fiz2,tz);
886 fjx1 = _mm_add_pd(fjx1,tx);
887 fjy1 = _mm_add_pd(fjy1,ty);
888 fjz1 = _mm_add_pd(fjz1,tz);
890 /**************************
891 * CALCULATE INTERACTIONS *
892 **************************/
894 /* REACTION-FIELD ELECTROSTATICS */
895 velec = _mm_mul_pd(qq22,_mm_sub_pd(_mm_add_pd(rinv22,_mm_mul_pd(krf,rsq22)),crf));
896 felec = _mm_mul_pd(qq22,_mm_sub_pd(_mm_mul_pd(rinv22,rinvsq22),krf2));
898 /* Update potential sum for this i atom from the interaction with this j atom. */
899 velec = _mm_unpacklo_pd(velec,_mm_setzero_pd());
900 velecsum = _mm_add_pd(velecsum,velec);
904 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
906 /* Calculate temporary vectorial force */
907 tx = _mm_mul_pd(fscal,dx22);
908 ty = _mm_mul_pd(fscal,dy22);
909 tz = _mm_mul_pd(fscal,dz22);
911 /* Update vectorial force */
912 fix2 = _mm_add_pd(fix2,tx);
913 fiy2 = _mm_add_pd(fiy2,ty);
914 fiz2 = _mm_add_pd(fiz2,tz);
916 fjx2 = _mm_add_pd(fjx2,tx);
917 fjy2 = _mm_add_pd(fjy2,ty);
918 fjz2 = _mm_add_pd(fjz2,tz);
920 /**************************
921 * CALCULATE INTERACTIONS *
922 **************************/
924 /* REACTION-FIELD ELECTROSTATICS */
925 velec = _mm_mul_pd(qq23,_mm_sub_pd(_mm_add_pd(rinv23,_mm_mul_pd(krf,rsq23)),crf));
926 felec = _mm_mul_pd(qq23,_mm_sub_pd(_mm_mul_pd(rinv23,rinvsq23),krf2));
928 /* Update potential sum for this i atom from the interaction with this j atom. */
929 velec = _mm_unpacklo_pd(velec,_mm_setzero_pd());
930 velecsum = _mm_add_pd(velecsum,velec);
934 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
936 /* Calculate temporary vectorial force */
937 tx = _mm_mul_pd(fscal,dx23);
938 ty = _mm_mul_pd(fscal,dy23);
939 tz = _mm_mul_pd(fscal,dz23);
941 /* Update vectorial force */
942 fix2 = _mm_add_pd(fix2,tx);
943 fiy2 = _mm_add_pd(fiy2,ty);
944 fiz2 = _mm_add_pd(fiz2,tz);
946 fjx3 = _mm_add_pd(fjx3,tx);
947 fjy3 = _mm_add_pd(fjy3,ty);
948 fjz3 = _mm_add_pd(fjz3,tz);
950 /**************************
951 * CALCULATE INTERACTIONS *
952 **************************/
954 /* REACTION-FIELD ELECTROSTATICS */
955 velec = _mm_mul_pd(qq31,_mm_sub_pd(_mm_add_pd(rinv31,_mm_mul_pd(krf,rsq31)),crf));
956 felec = _mm_mul_pd(qq31,_mm_sub_pd(_mm_mul_pd(rinv31,rinvsq31),krf2));
958 /* Update potential sum for this i atom from the interaction with this j atom. */
959 velec = _mm_unpacklo_pd(velec,_mm_setzero_pd());
960 velecsum = _mm_add_pd(velecsum,velec);
964 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
966 /* Calculate temporary vectorial force */
967 tx = _mm_mul_pd(fscal,dx31);
968 ty = _mm_mul_pd(fscal,dy31);
969 tz = _mm_mul_pd(fscal,dz31);
971 /* Update vectorial force */
972 fix3 = _mm_add_pd(fix3,tx);
973 fiy3 = _mm_add_pd(fiy3,ty);
974 fiz3 = _mm_add_pd(fiz3,tz);
976 fjx1 = _mm_add_pd(fjx1,tx);
977 fjy1 = _mm_add_pd(fjy1,ty);
978 fjz1 = _mm_add_pd(fjz1,tz);
980 /**************************
981 * CALCULATE INTERACTIONS *
982 **************************/
984 /* REACTION-FIELD ELECTROSTATICS */
985 velec = _mm_mul_pd(qq32,_mm_sub_pd(_mm_add_pd(rinv32,_mm_mul_pd(krf,rsq32)),crf));
986 felec = _mm_mul_pd(qq32,_mm_sub_pd(_mm_mul_pd(rinv32,rinvsq32),krf2));
988 /* Update potential sum for this i atom from the interaction with this j atom. */
989 velec = _mm_unpacklo_pd(velec,_mm_setzero_pd());
990 velecsum = _mm_add_pd(velecsum,velec);
994 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
996 /* Calculate temporary vectorial force */
997 tx = _mm_mul_pd(fscal,dx32);
998 ty = _mm_mul_pd(fscal,dy32);
999 tz = _mm_mul_pd(fscal,dz32);
1001 /* Update vectorial force */
1002 fix3 = _mm_add_pd(fix3,tx);
1003 fiy3 = _mm_add_pd(fiy3,ty);
1004 fiz3 = _mm_add_pd(fiz3,tz);
1006 fjx2 = _mm_add_pd(fjx2,tx);
1007 fjy2 = _mm_add_pd(fjy2,ty);
1008 fjz2 = _mm_add_pd(fjz2,tz);
1010 /**************************
1011 * CALCULATE INTERACTIONS *
1012 **************************/
1014 /* REACTION-FIELD ELECTROSTATICS */
1015 velec = _mm_mul_pd(qq33,_mm_sub_pd(_mm_add_pd(rinv33,_mm_mul_pd(krf,rsq33)),crf));
1016 felec = _mm_mul_pd(qq33,_mm_sub_pd(_mm_mul_pd(rinv33,rinvsq33),krf2));
1018 /* Update potential sum for this i atom from the interaction with this j atom. */
1019 velec = _mm_unpacklo_pd(velec,_mm_setzero_pd());
1020 velecsum = _mm_add_pd(velecsum,velec);
1024 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1026 /* Calculate temporary vectorial force */
1027 tx = _mm_mul_pd(fscal,dx33);
1028 ty = _mm_mul_pd(fscal,dy33);
1029 tz = _mm_mul_pd(fscal,dz33);
1031 /* Update vectorial force */
1032 fix3 = _mm_add_pd(fix3,tx);
1033 fiy3 = _mm_add_pd(fiy3,ty);
1034 fiz3 = _mm_add_pd(fiz3,tz);
1036 fjx3 = _mm_add_pd(fjx3,tx);
1037 fjy3 = _mm_add_pd(fjy3,ty);
1038 fjz3 = _mm_add_pd(fjz3,tz);
1040 gmx_mm_decrement_4rvec_1ptr_swizzle_pd(f+j_coord_offsetA,fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,fjx2,fjy2,fjz2,fjx3,fjy3,fjz3);
1042 /* Inner loop uses 347 flops */
1045 /* End of innermost loop */
1047 gmx_mm_update_iforce_4atom_swizzle_pd(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
1048 f+i_coord_offset,fshift+i_shift_offset);
1051 /* Update potential energies */
1052 gmx_mm_update_1pot_pd(velecsum,kernel_data->energygrp_elec+ggid);
1053 gmx_mm_update_1pot_pd(vvdwsum,kernel_data->energygrp_vdw+ggid);
1055 /* Increment number of inner iterations */
1056 inneriter += j_index_end - j_index_start;
1058 /* Outer loop uses 26 flops */
1061 /* Increment number of outer iterations */
1064 /* Update outer/inner flops */
1066 inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4W4_VF,outeriter*26 + inneriter*347);
1069 * Gromacs nonbonded kernel: nb_kernel_ElecRF_VdwCSTab_GeomW4W4_F_sse2_double
1070 * Electrostatics interaction: ReactionField
1071 * VdW interaction: CubicSplineTable
1072 * Geometry: Water4-Water4
1073 * Calculate force/pot: Force
1076 nb_kernel_ElecRF_VdwCSTab_GeomW4W4_F_sse2_double
1077 (t_nblist * gmx_restrict nlist,
1078 rvec * gmx_restrict xx,
1079 rvec * gmx_restrict ff,
1080 t_forcerec * gmx_restrict fr,
1081 t_mdatoms * gmx_restrict mdatoms,
1082 nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
1083 t_nrnb * gmx_restrict nrnb)
1085 /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
1086 * just 0 for non-waters.
1087 * Suffixes A,B refer to j loop unrolling done with SSE double precision, e.g. for the two different
1088 * jnr indices corresponding to data put in the four positions in the SIMD register.
1090 int i_shift_offset,i_coord_offset,outeriter,inneriter;
1091 int j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
1093 int j_coord_offsetA,j_coord_offsetB;
1094 int *iinr,*jindex,*jjnr,*shiftidx,*gid;
1095 real rcutoff_scalar;
1096 real *shiftvec,*fshift,*x,*f;
1097 __m128d tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
1099 __m128d ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
1101 __m128d ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
1103 __m128d ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
1105 __m128d ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
1106 int vdwjidx0A,vdwjidx0B;
1107 __m128d jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
1108 int vdwjidx1A,vdwjidx1B;
1109 __m128d jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
1110 int vdwjidx2A,vdwjidx2B;
1111 __m128d jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
1112 int vdwjidx3A,vdwjidx3B;
1113 __m128d jx3,jy3,jz3,fjx3,fjy3,fjz3,jq3,isaj3;
1114 __m128d dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
1115 __m128d dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11;
1116 __m128d dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12;
1117 __m128d dx13,dy13,dz13,rsq13,rinv13,rinvsq13,r13,qq13,c6_13,c12_13;
1118 __m128d dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21;
1119 __m128d dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22;
1120 __m128d dx23,dy23,dz23,rsq23,rinv23,rinvsq23,r23,qq23,c6_23,c12_23;
1121 __m128d dx31,dy31,dz31,rsq31,rinv31,rinvsq31,r31,qq31,c6_31,c12_31;
1122 __m128d dx32,dy32,dz32,rsq32,rinv32,rinvsq32,r32,qq32,c6_32,c12_32;
1123 __m128d dx33,dy33,dz33,rsq33,rinv33,rinvsq33,r33,qq33,c6_33,c12_33;
1124 __m128d velec,felec,velecsum,facel,crf,krf,krf2;
1127 __m128d rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
1130 __m128d one_sixth = _mm_set1_pd(1.0/6.0);
1131 __m128d one_twelfth = _mm_set1_pd(1.0/12.0);
1133 __m128i ifour = _mm_set1_epi32(4);
1134 __m128d rt,vfeps,vftabscale,Y,F,G,H,Heps,Fp,VV,FF;
1136 __m128d dummy_mask,cutoff_mask;
1137 __m128d signbit = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
1138 __m128d one = _mm_set1_pd(1.0);
1139 __m128d two = _mm_set1_pd(2.0);
1145 jindex = nlist->jindex;
1147 shiftidx = nlist->shift;
1149 shiftvec = fr->shift_vec[0];
1150 fshift = fr->fshift[0];
1151 facel = _mm_set1_pd(fr->epsfac);
1152 charge = mdatoms->chargeA;
1153 krf = _mm_set1_pd(fr->ic->k_rf);
1154 krf2 = _mm_set1_pd(fr->ic->k_rf*2.0);
1155 crf = _mm_set1_pd(fr->ic->c_rf);
1156 nvdwtype = fr->ntype;
1157 vdwparam = fr->nbfp;
1158 vdwtype = mdatoms->typeA;
1160 vftab = kernel_data->table_vdw->data;
1161 vftabscale = _mm_set1_pd(kernel_data->table_vdw->scale);
1163 /* Setup water-specific parameters */
1164 inr = nlist->iinr[0];
1165 iq1 = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+1]));
1166 iq2 = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+2]));
1167 iq3 = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+3]));
1168 vdwioffset0 = 2*nvdwtype*vdwtype[inr+0];
1170 jq1 = _mm_set1_pd(charge[inr+1]);
1171 jq2 = _mm_set1_pd(charge[inr+2]);
1172 jq3 = _mm_set1_pd(charge[inr+3]);
1173 vdwjidx0A = 2*vdwtype[inr+0];
1174 c6_00 = _mm_set1_pd(vdwparam[vdwioffset0+vdwjidx0A]);
1175 c12_00 = _mm_set1_pd(vdwparam[vdwioffset0+vdwjidx0A+1]);
1176 qq11 = _mm_mul_pd(iq1,jq1);
1177 qq12 = _mm_mul_pd(iq1,jq2);
1178 qq13 = _mm_mul_pd(iq1,jq3);
1179 qq21 = _mm_mul_pd(iq2,jq1);
1180 qq22 = _mm_mul_pd(iq2,jq2);
1181 qq23 = _mm_mul_pd(iq2,jq3);
1182 qq31 = _mm_mul_pd(iq3,jq1);
1183 qq32 = _mm_mul_pd(iq3,jq2);
1184 qq33 = _mm_mul_pd(iq3,jq3);
1186 /* Avoid stupid compiler warnings */
1188 j_coord_offsetA = 0;
1189 j_coord_offsetB = 0;
1194 /* Start outer loop over neighborlists */
1195 for(iidx=0; iidx<nri; iidx++)
1197 /* Load shift vector for this list */
1198 i_shift_offset = DIM*shiftidx[iidx];
1200 /* Load limits for loop over neighbors */
1201 j_index_start = jindex[iidx];
1202 j_index_end = jindex[iidx+1];
1204 /* Get outer coordinate index */
1206 i_coord_offset = DIM*inr;
1208 /* Load i particle coords and add shift vector */
1209 gmx_mm_load_shift_and_4rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,
1210 &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
1212 fix0 = _mm_setzero_pd();
1213 fiy0 = _mm_setzero_pd();
1214 fiz0 = _mm_setzero_pd();
1215 fix1 = _mm_setzero_pd();
1216 fiy1 = _mm_setzero_pd();
1217 fiz1 = _mm_setzero_pd();
1218 fix2 = _mm_setzero_pd();
1219 fiy2 = _mm_setzero_pd();
1220 fiz2 = _mm_setzero_pd();
1221 fix3 = _mm_setzero_pd();
1222 fiy3 = _mm_setzero_pd();
1223 fiz3 = _mm_setzero_pd();
1225 /* Start inner kernel loop */
1226 for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
1229 /* Get j neighbor index, and coordinate index */
1231 jnrB = jjnr[jidx+1];
1232 j_coord_offsetA = DIM*jnrA;
1233 j_coord_offsetB = DIM*jnrB;
1235 /* load j atom coordinates */
1236 gmx_mm_load_4rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
1237 &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,
1238 &jy2,&jz2,&jx3,&jy3,&jz3);
1240 /* Calculate displacement vector */
1241 dx00 = _mm_sub_pd(ix0,jx0);
1242 dy00 = _mm_sub_pd(iy0,jy0);
1243 dz00 = _mm_sub_pd(iz0,jz0);
1244 dx11 = _mm_sub_pd(ix1,jx1);
1245 dy11 = _mm_sub_pd(iy1,jy1);
1246 dz11 = _mm_sub_pd(iz1,jz1);
1247 dx12 = _mm_sub_pd(ix1,jx2);
1248 dy12 = _mm_sub_pd(iy1,jy2);
1249 dz12 = _mm_sub_pd(iz1,jz2);
1250 dx13 = _mm_sub_pd(ix1,jx3);
1251 dy13 = _mm_sub_pd(iy1,jy3);
1252 dz13 = _mm_sub_pd(iz1,jz3);
1253 dx21 = _mm_sub_pd(ix2,jx1);
1254 dy21 = _mm_sub_pd(iy2,jy1);
1255 dz21 = _mm_sub_pd(iz2,jz1);
1256 dx22 = _mm_sub_pd(ix2,jx2);
1257 dy22 = _mm_sub_pd(iy2,jy2);
1258 dz22 = _mm_sub_pd(iz2,jz2);
1259 dx23 = _mm_sub_pd(ix2,jx3);
1260 dy23 = _mm_sub_pd(iy2,jy3);
1261 dz23 = _mm_sub_pd(iz2,jz3);
1262 dx31 = _mm_sub_pd(ix3,jx1);
1263 dy31 = _mm_sub_pd(iy3,jy1);
1264 dz31 = _mm_sub_pd(iz3,jz1);
1265 dx32 = _mm_sub_pd(ix3,jx2);
1266 dy32 = _mm_sub_pd(iy3,jy2);
1267 dz32 = _mm_sub_pd(iz3,jz2);
1268 dx33 = _mm_sub_pd(ix3,jx3);
1269 dy33 = _mm_sub_pd(iy3,jy3);
1270 dz33 = _mm_sub_pd(iz3,jz3);
1272 /* Calculate squared distance and things based on it */
1273 rsq00 = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
1274 rsq11 = gmx_mm_calc_rsq_pd(dx11,dy11,dz11);
1275 rsq12 = gmx_mm_calc_rsq_pd(dx12,dy12,dz12);
1276 rsq13 = gmx_mm_calc_rsq_pd(dx13,dy13,dz13);
1277 rsq21 = gmx_mm_calc_rsq_pd(dx21,dy21,dz21);
1278 rsq22 = gmx_mm_calc_rsq_pd(dx22,dy22,dz22);
1279 rsq23 = gmx_mm_calc_rsq_pd(dx23,dy23,dz23);
1280 rsq31 = gmx_mm_calc_rsq_pd(dx31,dy31,dz31);
1281 rsq32 = gmx_mm_calc_rsq_pd(dx32,dy32,dz32);
1282 rsq33 = gmx_mm_calc_rsq_pd(dx33,dy33,dz33);
1284 rinv00 = gmx_mm_invsqrt_pd(rsq00);
1285 rinv11 = gmx_mm_invsqrt_pd(rsq11);
1286 rinv12 = gmx_mm_invsqrt_pd(rsq12);
1287 rinv13 = gmx_mm_invsqrt_pd(rsq13);
1288 rinv21 = gmx_mm_invsqrt_pd(rsq21);
1289 rinv22 = gmx_mm_invsqrt_pd(rsq22);
1290 rinv23 = gmx_mm_invsqrt_pd(rsq23);
1291 rinv31 = gmx_mm_invsqrt_pd(rsq31);
1292 rinv32 = gmx_mm_invsqrt_pd(rsq32);
1293 rinv33 = gmx_mm_invsqrt_pd(rsq33);
1295 rinvsq11 = _mm_mul_pd(rinv11,rinv11);
1296 rinvsq12 = _mm_mul_pd(rinv12,rinv12);
1297 rinvsq13 = _mm_mul_pd(rinv13,rinv13);
1298 rinvsq21 = _mm_mul_pd(rinv21,rinv21);
1299 rinvsq22 = _mm_mul_pd(rinv22,rinv22);
1300 rinvsq23 = _mm_mul_pd(rinv23,rinv23);
1301 rinvsq31 = _mm_mul_pd(rinv31,rinv31);
1302 rinvsq32 = _mm_mul_pd(rinv32,rinv32);
1303 rinvsq33 = _mm_mul_pd(rinv33,rinv33);
1305 fjx0 = _mm_setzero_pd();
1306 fjy0 = _mm_setzero_pd();
1307 fjz0 = _mm_setzero_pd();
1308 fjx1 = _mm_setzero_pd();
1309 fjy1 = _mm_setzero_pd();
1310 fjz1 = _mm_setzero_pd();
1311 fjx2 = _mm_setzero_pd();
1312 fjy2 = _mm_setzero_pd();
1313 fjz2 = _mm_setzero_pd();
1314 fjx3 = _mm_setzero_pd();
1315 fjy3 = _mm_setzero_pd();
1316 fjz3 = _mm_setzero_pd();
1318 /**************************
1319 * CALCULATE INTERACTIONS *
1320 **************************/
1322 r00 = _mm_mul_pd(rsq00,rinv00);
1324 /* Calculate table index by multiplying r with table scale and truncate to integer */
1325 rt = _mm_mul_pd(r00,vftabscale);
1326 vfitab = _mm_cvttpd_epi32(rt);
1327 vfeps = _mm_sub_pd(rt,_mm_cvtepi32_pd(vfitab));
1328 vfitab = _mm_slli_epi32(vfitab,3);
1330 /* CUBIC SPLINE TABLE DISPERSION */
1331 Y = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) );
1332 F = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,1) );
1333 GMX_MM_TRANSPOSE2_PD(Y,F);
1334 G = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) +2);
1335 H = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,1) +2);
1336 GMX_MM_TRANSPOSE2_PD(G,H);
1337 Heps = _mm_mul_pd(vfeps,H);
1338 Fp = _mm_add_pd(F,_mm_mul_pd(vfeps,_mm_add_pd(G,Heps)));
1339 FF = _mm_add_pd(Fp,_mm_mul_pd(vfeps,_mm_add_pd(G,_mm_add_pd(Heps,Heps))));
1340 fvdw6 = _mm_mul_pd(c6_00,FF);
1342 /* CUBIC SPLINE TABLE REPULSION */
1343 vfitab = _mm_add_epi32(vfitab,ifour);
1344 Y = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) );
1345 F = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,1) );
1346 GMX_MM_TRANSPOSE2_PD(Y,F);
1347 G = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) +2);
1348 H = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,1) +2);
1349 GMX_MM_TRANSPOSE2_PD(G,H);
1350 Heps = _mm_mul_pd(vfeps,H);
1351 Fp = _mm_add_pd(F,_mm_mul_pd(vfeps,_mm_add_pd(G,Heps)));
1352 FF = _mm_add_pd(Fp,_mm_mul_pd(vfeps,_mm_add_pd(G,_mm_add_pd(Heps,Heps))));
1353 fvdw12 = _mm_mul_pd(c12_00,FF);
1354 fvdw = _mm_xor_pd(signbit,_mm_mul_pd(_mm_add_pd(fvdw6,fvdw12),_mm_mul_pd(vftabscale,rinv00)));
1358 /* Calculate temporary vectorial force */
1359 tx = _mm_mul_pd(fscal,dx00);
1360 ty = _mm_mul_pd(fscal,dy00);
1361 tz = _mm_mul_pd(fscal,dz00);
1363 /* Update vectorial force */
1364 fix0 = _mm_add_pd(fix0,tx);
1365 fiy0 = _mm_add_pd(fiy0,ty);
1366 fiz0 = _mm_add_pd(fiz0,tz);
1368 fjx0 = _mm_add_pd(fjx0,tx);
1369 fjy0 = _mm_add_pd(fjy0,ty);
1370 fjz0 = _mm_add_pd(fjz0,tz);
1372 /**************************
1373 * CALCULATE INTERACTIONS *
1374 **************************/
1376 /* REACTION-FIELD ELECTROSTATICS */
1377 felec = _mm_mul_pd(qq11,_mm_sub_pd(_mm_mul_pd(rinv11,rinvsq11),krf2));
1381 /* Calculate temporary vectorial force */
1382 tx = _mm_mul_pd(fscal,dx11);
1383 ty = _mm_mul_pd(fscal,dy11);
1384 tz = _mm_mul_pd(fscal,dz11);
1386 /* Update vectorial force */
1387 fix1 = _mm_add_pd(fix1,tx);
1388 fiy1 = _mm_add_pd(fiy1,ty);
1389 fiz1 = _mm_add_pd(fiz1,tz);
1391 fjx1 = _mm_add_pd(fjx1,tx);
1392 fjy1 = _mm_add_pd(fjy1,ty);
1393 fjz1 = _mm_add_pd(fjz1,tz);
1395 /**************************
1396 * CALCULATE INTERACTIONS *
1397 **************************/
1399 /* REACTION-FIELD ELECTROSTATICS */
1400 felec = _mm_mul_pd(qq12,_mm_sub_pd(_mm_mul_pd(rinv12,rinvsq12),krf2));
1404 /* Calculate temporary vectorial force */
1405 tx = _mm_mul_pd(fscal,dx12);
1406 ty = _mm_mul_pd(fscal,dy12);
1407 tz = _mm_mul_pd(fscal,dz12);
1409 /* Update vectorial force */
1410 fix1 = _mm_add_pd(fix1,tx);
1411 fiy1 = _mm_add_pd(fiy1,ty);
1412 fiz1 = _mm_add_pd(fiz1,tz);
1414 fjx2 = _mm_add_pd(fjx2,tx);
1415 fjy2 = _mm_add_pd(fjy2,ty);
1416 fjz2 = _mm_add_pd(fjz2,tz);
1418 /**************************
1419 * CALCULATE INTERACTIONS *
1420 **************************/
1422 /* REACTION-FIELD ELECTROSTATICS */
1423 felec = _mm_mul_pd(qq13,_mm_sub_pd(_mm_mul_pd(rinv13,rinvsq13),krf2));
1427 /* Calculate temporary vectorial force */
1428 tx = _mm_mul_pd(fscal,dx13);
1429 ty = _mm_mul_pd(fscal,dy13);
1430 tz = _mm_mul_pd(fscal,dz13);
1432 /* Update vectorial force */
1433 fix1 = _mm_add_pd(fix1,tx);
1434 fiy1 = _mm_add_pd(fiy1,ty);
1435 fiz1 = _mm_add_pd(fiz1,tz);
1437 fjx3 = _mm_add_pd(fjx3,tx);
1438 fjy3 = _mm_add_pd(fjy3,ty);
1439 fjz3 = _mm_add_pd(fjz3,tz);
1441 /**************************
1442 * CALCULATE INTERACTIONS *
1443 **************************/
1445 /* REACTION-FIELD ELECTROSTATICS */
1446 felec = _mm_mul_pd(qq21,_mm_sub_pd(_mm_mul_pd(rinv21,rinvsq21),krf2));
1450 /* Calculate temporary vectorial force */
1451 tx = _mm_mul_pd(fscal,dx21);
1452 ty = _mm_mul_pd(fscal,dy21);
1453 tz = _mm_mul_pd(fscal,dz21);
1455 /* Update vectorial force */
1456 fix2 = _mm_add_pd(fix2,tx);
1457 fiy2 = _mm_add_pd(fiy2,ty);
1458 fiz2 = _mm_add_pd(fiz2,tz);
1460 fjx1 = _mm_add_pd(fjx1,tx);
1461 fjy1 = _mm_add_pd(fjy1,ty);
1462 fjz1 = _mm_add_pd(fjz1,tz);
1464 /**************************
1465 * CALCULATE INTERACTIONS *
1466 **************************/
1468 /* REACTION-FIELD ELECTROSTATICS */
1469 felec = _mm_mul_pd(qq22,_mm_sub_pd(_mm_mul_pd(rinv22,rinvsq22),krf2));
1473 /* Calculate temporary vectorial force */
1474 tx = _mm_mul_pd(fscal,dx22);
1475 ty = _mm_mul_pd(fscal,dy22);
1476 tz = _mm_mul_pd(fscal,dz22);
1478 /* Update vectorial force */
1479 fix2 = _mm_add_pd(fix2,tx);
1480 fiy2 = _mm_add_pd(fiy2,ty);
1481 fiz2 = _mm_add_pd(fiz2,tz);
1483 fjx2 = _mm_add_pd(fjx2,tx);
1484 fjy2 = _mm_add_pd(fjy2,ty);
1485 fjz2 = _mm_add_pd(fjz2,tz);
1487 /**************************
1488 * CALCULATE INTERACTIONS *
1489 **************************/
1491 /* REACTION-FIELD ELECTROSTATICS */
1492 felec = _mm_mul_pd(qq23,_mm_sub_pd(_mm_mul_pd(rinv23,rinvsq23),krf2));
1496 /* Calculate temporary vectorial force */
1497 tx = _mm_mul_pd(fscal,dx23);
1498 ty = _mm_mul_pd(fscal,dy23);
1499 tz = _mm_mul_pd(fscal,dz23);
1501 /* Update vectorial force */
1502 fix2 = _mm_add_pd(fix2,tx);
1503 fiy2 = _mm_add_pd(fiy2,ty);
1504 fiz2 = _mm_add_pd(fiz2,tz);
1506 fjx3 = _mm_add_pd(fjx3,tx);
1507 fjy3 = _mm_add_pd(fjy3,ty);
1508 fjz3 = _mm_add_pd(fjz3,tz);
1510 /**************************
1511 * CALCULATE INTERACTIONS *
1512 **************************/
1514 /* REACTION-FIELD ELECTROSTATICS */
1515 felec = _mm_mul_pd(qq31,_mm_sub_pd(_mm_mul_pd(rinv31,rinvsq31),krf2));
1519 /* Calculate temporary vectorial force */
1520 tx = _mm_mul_pd(fscal,dx31);
1521 ty = _mm_mul_pd(fscal,dy31);
1522 tz = _mm_mul_pd(fscal,dz31);
1524 /* Update vectorial force */
1525 fix3 = _mm_add_pd(fix3,tx);
1526 fiy3 = _mm_add_pd(fiy3,ty);
1527 fiz3 = _mm_add_pd(fiz3,tz);
1529 fjx1 = _mm_add_pd(fjx1,tx);
1530 fjy1 = _mm_add_pd(fjy1,ty);
1531 fjz1 = _mm_add_pd(fjz1,tz);
1533 /**************************
1534 * CALCULATE INTERACTIONS *
1535 **************************/
1537 /* REACTION-FIELD ELECTROSTATICS */
1538 felec = _mm_mul_pd(qq32,_mm_sub_pd(_mm_mul_pd(rinv32,rinvsq32),krf2));
1542 /* Calculate temporary vectorial force */
1543 tx = _mm_mul_pd(fscal,dx32);
1544 ty = _mm_mul_pd(fscal,dy32);
1545 tz = _mm_mul_pd(fscal,dz32);
1547 /* Update vectorial force */
1548 fix3 = _mm_add_pd(fix3,tx);
1549 fiy3 = _mm_add_pd(fiy3,ty);
1550 fiz3 = _mm_add_pd(fiz3,tz);
1552 fjx2 = _mm_add_pd(fjx2,tx);
1553 fjy2 = _mm_add_pd(fjy2,ty);
1554 fjz2 = _mm_add_pd(fjz2,tz);
1556 /**************************
1557 * CALCULATE INTERACTIONS *
1558 **************************/
1560 /* REACTION-FIELD ELECTROSTATICS */
1561 felec = _mm_mul_pd(qq33,_mm_sub_pd(_mm_mul_pd(rinv33,rinvsq33),krf2));
1565 /* Calculate temporary vectorial force */
1566 tx = _mm_mul_pd(fscal,dx33);
1567 ty = _mm_mul_pd(fscal,dy33);
1568 tz = _mm_mul_pd(fscal,dz33);
1570 /* Update vectorial force */
1571 fix3 = _mm_add_pd(fix3,tx);
1572 fiy3 = _mm_add_pd(fiy3,ty);
1573 fiz3 = _mm_add_pd(fiz3,tz);
1575 fjx3 = _mm_add_pd(fjx3,tx);
1576 fjy3 = _mm_add_pd(fjy3,ty);
1577 fjz3 = _mm_add_pd(fjz3,tz);
1579 gmx_mm_decrement_4rvec_2ptr_swizzle_pd(f+j_coord_offsetA,f+j_coord_offsetB,fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,fjx2,fjy2,fjz2,fjx3,fjy3,fjz3);
1581 /* Inner loop uses 294 flops */
1584 if(jidx<j_index_end)
1588 j_coord_offsetA = DIM*jnrA;
1590 /* load j atom coordinates */
1591 gmx_mm_load_4rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
1592 &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,
1593 &jy2,&jz2,&jx3,&jy3,&jz3);
1595 /* Calculate displacement vector */
1596 dx00 = _mm_sub_pd(ix0,jx0);
1597 dy00 = _mm_sub_pd(iy0,jy0);
1598 dz00 = _mm_sub_pd(iz0,jz0);
1599 dx11 = _mm_sub_pd(ix1,jx1);
1600 dy11 = _mm_sub_pd(iy1,jy1);
1601 dz11 = _mm_sub_pd(iz1,jz1);
1602 dx12 = _mm_sub_pd(ix1,jx2);
1603 dy12 = _mm_sub_pd(iy1,jy2);
1604 dz12 = _mm_sub_pd(iz1,jz2);
1605 dx13 = _mm_sub_pd(ix1,jx3);
1606 dy13 = _mm_sub_pd(iy1,jy3);
1607 dz13 = _mm_sub_pd(iz1,jz3);
1608 dx21 = _mm_sub_pd(ix2,jx1);
1609 dy21 = _mm_sub_pd(iy2,jy1);
1610 dz21 = _mm_sub_pd(iz2,jz1);
1611 dx22 = _mm_sub_pd(ix2,jx2);
1612 dy22 = _mm_sub_pd(iy2,jy2);
1613 dz22 = _mm_sub_pd(iz2,jz2);
1614 dx23 = _mm_sub_pd(ix2,jx3);
1615 dy23 = _mm_sub_pd(iy2,jy3);
1616 dz23 = _mm_sub_pd(iz2,jz3);
1617 dx31 = _mm_sub_pd(ix3,jx1);
1618 dy31 = _mm_sub_pd(iy3,jy1);
1619 dz31 = _mm_sub_pd(iz3,jz1);
1620 dx32 = _mm_sub_pd(ix3,jx2);
1621 dy32 = _mm_sub_pd(iy3,jy2);
1622 dz32 = _mm_sub_pd(iz3,jz2);
1623 dx33 = _mm_sub_pd(ix3,jx3);
1624 dy33 = _mm_sub_pd(iy3,jy3);
1625 dz33 = _mm_sub_pd(iz3,jz3);
1627 /* Calculate squared distance and things based on it */
1628 rsq00 = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
1629 rsq11 = gmx_mm_calc_rsq_pd(dx11,dy11,dz11);
1630 rsq12 = gmx_mm_calc_rsq_pd(dx12,dy12,dz12);
1631 rsq13 = gmx_mm_calc_rsq_pd(dx13,dy13,dz13);
1632 rsq21 = gmx_mm_calc_rsq_pd(dx21,dy21,dz21);
1633 rsq22 = gmx_mm_calc_rsq_pd(dx22,dy22,dz22);
1634 rsq23 = gmx_mm_calc_rsq_pd(dx23,dy23,dz23);
1635 rsq31 = gmx_mm_calc_rsq_pd(dx31,dy31,dz31);
1636 rsq32 = gmx_mm_calc_rsq_pd(dx32,dy32,dz32);
1637 rsq33 = gmx_mm_calc_rsq_pd(dx33,dy33,dz33);
1639 rinv00 = gmx_mm_invsqrt_pd(rsq00);
1640 rinv11 = gmx_mm_invsqrt_pd(rsq11);
1641 rinv12 = gmx_mm_invsqrt_pd(rsq12);
1642 rinv13 = gmx_mm_invsqrt_pd(rsq13);
1643 rinv21 = gmx_mm_invsqrt_pd(rsq21);
1644 rinv22 = gmx_mm_invsqrt_pd(rsq22);
1645 rinv23 = gmx_mm_invsqrt_pd(rsq23);
1646 rinv31 = gmx_mm_invsqrt_pd(rsq31);
1647 rinv32 = gmx_mm_invsqrt_pd(rsq32);
1648 rinv33 = gmx_mm_invsqrt_pd(rsq33);
1650 rinvsq11 = _mm_mul_pd(rinv11,rinv11);
1651 rinvsq12 = _mm_mul_pd(rinv12,rinv12);
1652 rinvsq13 = _mm_mul_pd(rinv13,rinv13);
1653 rinvsq21 = _mm_mul_pd(rinv21,rinv21);
1654 rinvsq22 = _mm_mul_pd(rinv22,rinv22);
1655 rinvsq23 = _mm_mul_pd(rinv23,rinv23);
1656 rinvsq31 = _mm_mul_pd(rinv31,rinv31);
1657 rinvsq32 = _mm_mul_pd(rinv32,rinv32);
1658 rinvsq33 = _mm_mul_pd(rinv33,rinv33);
1660 fjx0 = _mm_setzero_pd();
1661 fjy0 = _mm_setzero_pd();
1662 fjz0 = _mm_setzero_pd();
1663 fjx1 = _mm_setzero_pd();
1664 fjy1 = _mm_setzero_pd();
1665 fjz1 = _mm_setzero_pd();
1666 fjx2 = _mm_setzero_pd();
1667 fjy2 = _mm_setzero_pd();
1668 fjz2 = _mm_setzero_pd();
1669 fjx3 = _mm_setzero_pd();
1670 fjy3 = _mm_setzero_pd();
1671 fjz3 = _mm_setzero_pd();
1673 /**************************
1674 * CALCULATE INTERACTIONS *
1675 **************************/
1677 r00 = _mm_mul_pd(rsq00,rinv00);
1679 /* Calculate table index by multiplying r with table scale and truncate to integer */
1680 rt = _mm_mul_pd(r00,vftabscale);
1681 vfitab = _mm_cvttpd_epi32(rt);
1682 vfeps = _mm_sub_pd(rt,_mm_cvtepi32_pd(vfitab));
1683 vfitab = _mm_slli_epi32(vfitab,3);
1685 /* CUBIC SPLINE TABLE DISPERSION */
1686 Y = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) );
1687 F = _mm_setzero_pd();
1688 GMX_MM_TRANSPOSE2_PD(Y,F);
1689 G = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) +2);
1690 H = _mm_setzero_pd();
1691 GMX_MM_TRANSPOSE2_PD(G,H);
1692 Heps = _mm_mul_pd(vfeps,H);
1693 Fp = _mm_add_pd(F,_mm_mul_pd(vfeps,_mm_add_pd(G,Heps)));
1694 FF = _mm_add_pd(Fp,_mm_mul_pd(vfeps,_mm_add_pd(G,_mm_add_pd(Heps,Heps))));
1695 fvdw6 = _mm_mul_pd(c6_00,FF);
1697 /* CUBIC SPLINE TABLE REPULSION */
1698 vfitab = _mm_add_epi32(vfitab,ifour);
1699 Y = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) );
1700 F = _mm_setzero_pd();
1701 GMX_MM_TRANSPOSE2_PD(Y,F);
1702 G = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) +2);
1703 H = _mm_setzero_pd();
1704 GMX_MM_TRANSPOSE2_PD(G,H);
1705 Heps = _mm_mul_pd(vfeps,H);
1706 Fp = _mm_add_pd(F,_mm_mul_pd(vfeps,_mm_add_pd(G,Heps)));
1707 FF = _mm_add_pd(Fp,_mm_mul_pd(vfeps,_mm_add_pd(G,_mm_add_pd(Heps,Heps))));
1708 fvdw12 = _mm_mul_pd(c12_00,FF);
1709 fvdw = _mm_xor_pd(signbit,_mm_mul_pd(_mm_add_pd(fvdw6,fvdw12),_mm_mul_pd(vftabscale,rinv00)));
1713 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1715 /* Calculate temporary vectorial force */
1716 tx = _mm_mul_pd(fscal,dx00);
1717 ty = _mm_mul_pd(fscal,dy00);
1718 tz = _mm_mul_pd(fscal,dz00);
1720 /* Update vectorial force */
1721 fix0 = _mm_add_pd(fix0,tx);
1722 fiy0 = _mm_add_pd(fiy0,ty);
1723 fiz0 = _mm_add_pd(fiz0,tz);
1725 fjx0 = _mm_add_pd(fjx0,tx);
1726 fjy0 = _mm_add_pd(fjy0,ty);
1727 fjz0 = _mm_add_pd(fjz0,tz);
1729 /**************************
1730 * CALCULATE INTERACTIONS *
1731 **************************/
1733 /* REACTION-FIELD ELECTROSTATICS */
1734 felec = _mm_mul_pd(qq11,_mm_sub_pd(_mm_mul_pd(rinv11,rinvsq11),krf2));
1738 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1740 /* Calculate temporary vectorial force */
1741 tx = _mm_mul_pd(fscal,dx11);
1742 ty = _mm_mul_pd(fscal,dy11);
1743 tz = _mm_mul_pd(fscal,dz11);
1745 /* Update vectorial force */
1746 fix1 = _mm_add_pd(fix1,tx);
1747 fiy1 = _mm_add_pd(fiy1,ty);
1748 fiz1 = _mm_add_pd(fiz1,tz);
1750 fjx1 = _mm_add_pd(fjx1,tx);
1751 fjy1 = _mm_add_pd(fjy1,ty);
1752 fjz1 = _mm_add_pd(fjz1,tz);
1754 /**************************
1755 * CALCULATE INTERACTIONS *
1756 **************************/
1758 /* REACTION-FIELD ELECTROSTATICS */
1759 felec = _mm_mul_pd(qq12,_mm_sub_pd(_mm_mul_pd(rinv12,rinvsq12),krf2));
1763 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1765 /* Calculate temporary vectorial force */
1766 tx = _mm_mul_pd(fscal,dx12);
1767 ty = _mm_mul_pd(fscal,dy12);
1768 tz = _mm_mul_pd(fscal,dz12);
1770 /* Update vectorial force */
1771 fix1 = _mm_add_pd(fix1,tx);
1772 fiy1 = _mm_add_pd(fiy1,ty);
1773 fiz1 = _mm_add_pd(fiz1,tz);
1775 fjx2 = _mm_add_pd(fjx2,tx);
1776 fjy2 = _mm_add_pd(fjy2,ty);
1777 fjz2 = _mm_add_pd(fjz2,tz);
1779 /**************************
1780 * CALCULATE INTERACTIONS *
1781 **************************/
1783 /* REACTION-FIELD ELECTROSTATICS */
1784 felec = _mm_mul_pd(qq13,_mm_sub_pd(_mm_mul_pd(rinv13,rinvsq13),krf2));
1788 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1790 /* Calculate temporary vectorial force */
1791 tx = _mm_mul_pd(fscal,dx13);
1792 ty = _mm_mul_pd(fscal,dy13);
1793 tz = _mm_mul_pd(fscal,dz13);
1795 /* Update vectorial force */
1796 fix1 = _mm_add_pd(fix1,tx);
1797 fiy1 = _mm_add_pd(fiy1,ty);
1798 fiz1 = _mm_add_pd(fiz1,tz);
1800 fjx3 = _mm_add_pd(fjx3,tx);
1801 fjy3 = _mm_add_pd(fjy3,ty);
1802 fjz3 = _mm_add_pd(fjz3,tz);
1804 /**************************
1805 * CALCULATE INTERACTIONS *
1806 **************************/
1808 /* REACTION-FIELD ELECTROSTATICS */
1809 felec = _mm_mul_pd(qq21,_mm_sub_pd(_mm_mul_pd(rinv21,rinvsq21),krf2));
1813 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1815 /* Calculate temporary vectorial force */
1816 tx = _mm_mul_pd(fscal,dx21);
1817 ty = _mm_mul_pd(fscal,dy21);
1818 tz = _mm_mul_pd(fscal,dz21);
1820 /* Update vectorial force */
1821 fix2 = _mm_add_pd(fix2,tx);
1822 fiy2 = _mm_add_pd(fiy2,ty);
1823 fiz2 = _mm_add_pd(fiz2,tz);
1825 fjx1 = _mm_add_pd(fjx1,tx);
1826 fjy1 = _mm_add_pd(fjy1,ty);
1827 fjz1 = _mm_add_pd(fjz1,tz);
1829 /**************************
1830 * CALCULATE INTERACTIONS *
1831 **************************/
1833 /* REACTION-FIELD ELECTROSTATICS */
1834 felec = _mm_mul_pd(qq22,_mm_sub_pd(_mm_mul_pd(rinv22,rinvsq22),krf2));
1838 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1840 /* Calculate temporary vectorial force */
1841 tx = _mm_mul_pd(fscal,dx22);
1842 ty = _mm_mul_pd(fscal,dy22);
1843 tz = _mm_mul_pd(fscal,dz22);
1845 /* Update vectorial force */
1846 fix2 = _mm_add_pd(fix2,tx);
1847 fiy2 = _mm_add_pd(fiy2,ty);
1848 fiz2 = _mm_add_pd(fiz2,tz);
1850 fjx2 = _mm_add_pd(fjx2,tx);
1851 fjy2 = _mm_add_pd(fjy2,ty);
1852 fjz2 = _mm_add_pd(fjz2,tz);
1854 /**************************
1855 * CALCULATE INTERACTIONS *
1856 **************************/
1858 /* REACTION-FIELD ELECTROSTATICS */
1859 felec = _mm_mul_pd(qq23,_mm_sub_pd(_mm_mul_pd(rinv23,rinvsq23),krf2));
1863 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1865 /* Calculate temporary vectorial force */
1866 tx = _mm_mul_pd(fscal,dx23);
1867 ty = _mm_mul_pd(fscal,dy23);
1868 tz = _mm_mul_pd(fscal,dz23);
1870 /* Update vectorial force */
1871 fix2 = _mm_add_pd(fix2,tx);
1872 fiy2 = _mm_add_pd(fiy2,ty);
1873 fiz2 = _mm_add_pd(fiz2,tz);
1875 fjx3 = _mm_add_pd(fjx3,tx);
1876 fjy3 = _mm_add_pd(fjy3,ty);
1877 fjz3 = _mm_add_pd(fjz3,tz);
1879 /**************************
1880 * CALCULATE INTERACTIONS *
1881 **************************/
1883 /* REACTION-FIELD ELECTROSTATICS */
1884 felec = _mm_mul_pd(qq31,_mm_sub_pd(_mm_mul_pd(rinv31,rinvsq31),krf2));
1888 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1890 /* Calculate temporary vectorial force */
1891 tx = _mm_mul_pd(fscal,dx31);
1892 ty = _mm_mul_pd(fscal,dy31);
1893 tz = _mm_mul_pd(fscal,dz31);
1895 /* Update vectorial force */
1896 fix3 = _mm_add_pd(fix3,tx);
1897 fiy3 = _mm_add_pd(fiy3,ty);
1898 fiz3 = _mm_add_pd(fiz3,tz);
1900 fjx1 = _mm_add_pd(fjx1,tx);
1901 fjy1 = _mm_add_pd(fjy1,ty);
1902 fjz1 = _mm_add_pd(fjz1,tz);
1904 /**************************
1905 * CALCULATE INTERACTIONS *
1906 **************************/
1908 /* REACTION-FIELD ELECTROSTATICS */
1909 felec = _mm_mul_pd(qq32,_mm_sub_pd(_mm_mul_pd(rinv32,rinvsq32),krf2));
1913 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1915 /* Calculate temporary vectorial force */
1916 tx = _mm_mul_pd(fscal,dx32);
1917 ty = _mm_mul_pd(fscal,dy32);
1918 tz = _mm_mul_pd(fscal,dz32);
1920 /* Update vectorial force */
1921 fix3 = _mm_add_pd(fix3,tx);
1922 fiy3 = _mm_add_pd(fiy3,ty);
1923 fiz3 = _mm_add_pd(fiz3,tz);
1925 fjx2 = _mm_add_pd(fjx2,tx);
1926 fjy2 = _mm_add_pd(fjy2,ty);
1927 fjz2 = _mm_add_pd(fjz2,tz);
1929 /**************************
1930 * CALCULATE INTERACTIONS *
1931 **************************/
1933 /* REACTION-FIELD ELECTROSTATICS */
1934 felec = _mm_mul_pd(qq33,_mm_sub_pd(_mm_mul_pd(rinv33,rinvsq33),krf2));
1938 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1940 /* Calculate temporary vectorial force */
1941 tx = _mm_mul_pd(fscal,dx33);
1942 ty = _mm_mul_pd(fscal,dy33);
1943 tz = _mm_mul_pd(fscal,dz33);
1945 /* Update vectorial force */
1946 fix3 = _mm_add_pd(fix3,tx);
1947 fiy3 = _mm_add_pd(fiy3,ty);
1948 fiz3 = _mm_add_pd(fiz3,tz);
1950 fjx3 = _mm_add_pd(fjx3,tx);
1951 fjy3 = _mm_add_pd(fjy3,ty);
1952 fjz3 = _mm_add_pd(fjz3,tz);
1954 gmx_mm_decrement_4rvec_1ptr_swizzle_pd(f+j_coord_offsetA,fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,fjx2,fjy2,fjz2,fjx3,fjy3,fjz3);
1956 /* Inner loop uses 294 flops */
1959 /* End of innermost loop */
1961 gmx_mm_update_iforce_4atom_swizzle_pd(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
1962 f+i_coord_offset,fshift+i_shift_offset);
1964 /* Increment number of inner iterations */
1965 inneriter += j_index_end - j_index_start;
1967 /* Outer loop uses 24 flops */
1970 /* Increment number of outer iterations */
1973 /* Update outer/inner flops */
1975 inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4W4_F,outeriter*24 + inneriter*294);