2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013, 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 sse4_1_single kernel generator.
44 #include "../nb_kernel.h"
45 #include "types/simple.h"
49 #include "gromacs/simd/math_x86_sse4_1_single.h"
50 #include "kernelutil_x86_sse4_1_single.h"
53 * Gromacs nonbonded kernel: nb_kernel_ElecEwSw_VdwNone_GeomW3W3_VF_sse4_1_single
54 * Electrostatics interaction: Ewald
55 * VdW interaction: None
56 * Geometry: Water3-Water3
57 * Calculate force/pot: PotentialAndForce
60 nb_kernel_ElecEwSw_VdwNone_GeomW3W3_VF_sse4_1_single
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,C,D refer to j loop unrolling done with SSE, e.g. for the four 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;
76 int jnrA,jnrB,jnrC,jnrD;
77 int jnrlistA,jnrlistB,jnrlistC,jnrlistD;
78 int j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
79 int *iinr,*jindex,*jjnr,*shiftidx,*gid;
81 real *shiftvec,*fshift,*x,*f;
82 real *fjptrA,*fjptrB,*fjptrC,*fjptrD;
84 __m128 tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
86 __m128 ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
88 __m128 ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
90 __m128 ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
91 int vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
92 __m128 jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
93 int vdwjidx1A,vdwjidx1B,vdwjidx1C,vdwjidx1D;
94 __m128 jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
95 int vdwjidx2A,vdwjidx2B,vdwjidx2C,vdwjidx2D;
96 __m128 jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
97 __m128 dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
98 __m128 dx01,dy01,dz01,rsq01,rinv01,rinvsq01,r01,qq01,c6_01,c12_01;
99 __m128 dx02,dy02,dz02,rsq02,rinv02,rinvsq02,r02,qq02,c6_02,c12_02;
100 __m128 dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
101 __m128 dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11;
102 __m128 dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12;
103 __m128 dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
104 __m128 dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21;
105 __m128 dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22;
106 __m128 velec,felec,velecsum,facel,crf,krf,krf2;
109 __m128 ewtabscale,eweps,sh_ewald,ewrt,ewtabhalfspace,ewtabF,ewtabFn,ewtabD,ewtabV;
111 __m128 rswitch,swV3,swV4,swV5,swF2,swF3,swF4,d,d2,sw,dsw;
112 real rswitch_scalar,d_scalar;
113 __m128 dummy_mask,cutoff_mask;
114 __m128 signbit = _mm_castsi128_ps( _mm_set1_epi32(0x80000000) );
115 __m128 one = _mm_set1_ps(1.0);
116 __m128 two = _mm_set1_ps(2.0);
122 jindex = nlist->jindex;
124 shiftidx = nlist->shift;
126 shiftvec = fr->shift_vec[0];
127 fshift = fr->fshift[0];
128 facel = _mm_set1_ps(fr->epsfac);
129 charge = mdatoms->chargeA;
131 sh_ewald = _mm_set1_ps(fr->ic->sh_ewald);
132 ewtab = fr->ic->tabq_coul_FDV0;
133 ewtabscale = _mm_set1_ps(fr->ic->tabq_scale);
134 ewtabhalfspace = _mm_set1_ps(0.5/fr->ic->tabq_scale);
136 /* Setup water-specific parameters */
137 inr = nlist->iinr[0];
138 iq0 = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+0]));
139 iq1 = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+1]));
140 iq2 = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+2]));
142 jq0 = _mm_set1_ps(charge[inr+0]);
143 jq1 = _mm_set1_ps(charge[inr+1]);
144 jq2 = _mm_set1_ps(charge[inr+2]);
145 qq00 = _mm_mul_ps(iq0,jq0);
146 qq01 = _mm_mul_ps(iq0,jq1);
147 qq02 = _mm_mul_ps(iq0,jq2);
148 qq10 = _mm_mul_ps(iq1,jq0);
149 qq11 = _mm_mul_ps(iq1,jq1);
150 qq12 = _mm_mul_ps(iq1,jq2);
151 qq20 = _mm_mul_ps(iq2,jq0);
152 qq21 = _mm_mul_ps(iq2,jq1);
153 qq22 = _mm_mul_ps(iq2,jq2);
155 /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
156 rcutoff_scalar = fr->rcoulomb;
157 rcutoff = _mm_set1_ps(rcutoff_scalar);
158 rcutoff2 = _mm_mul_ps(rcutoff,rcutoff);
160 rswitch_scalar = fr->rcoulomb_switch;
161 rswitch = _mm_set1_ps(rswitch_scalar);
162 /* Setup switch parameters */
163 d_scalar = rcutoff_scalar-rswitch_scalar;
164 d = _mm_set1_ps(d_scalar);
165 swV3 = _mm_set1_ps(-10.0/(d_scalar*d_scalar*d_scalar));
166 swV4 = _mm_set1_ps( 15.0/(d_scalar*d_scalar*d_scalar*d_scalar));
167 swV5 = _mm_set1_ps( -6.0/(d_scalar*d_scalar*d_scalar*d_scalar*d_scalar));
168 swF2 = _mm_set1_ps(-30.0/(d_scalar*d_scalar*d_scalar));
169 swF3 = _mm_set1_ps( 60.0/(d_scalar*d_scalar*d_scalar*d_scalar));
170 swF4 = _mm_set1_ps(-30.0/(d_scalar*d_scalar*d_scalar*d_scalar*d_scalar));
172 /* Avoid stupid compiler warnings */
173 jnrA = jnrB = jnrC = jnrD = 0;
182 for(iidx=0;iidx<4*DIM;iidx++)
187 /* Start outer loop over neighborlists */
188 for(iidx=0; iidx<nri; iidx++)
190 /* Load shift vector for this list */
191 i_shift_offset = DIM*shiftidx[iidx];
193 /* Load limits for loop over neighbors */
194 j_index_start = jindex[iidx];
195 j_index_end = jindex[iidx+1];
197 /* Get outer coordinate index */
199 i_coord_offset = DIM*inr;
201 /* Load i particle coords and add shift vector */
202 gmx_mm_load_shift_and_3rvec_broadcast_ps(shiftvec+i_shift_offset,x+i_coord_offset,
203 &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2);
205 fix0 = _mm_setzero_ps();
206 fiy0 = _mm_setzero_ps();
207 fiz0 = _mm_setzero_ps();
208 fix1 = _mm_setzero_ps();
209 fiy1 = _mm_setzero_ps();
210 fiz1 = _mm_setzero_ps();
211 fix2 = _mm_setzero_ps();
212 fiy2 = _mm_setzero_ps();
213 fiz2 = _mm_setzero_ps();
215 /* Reset potential sums */
216 velecsum = _mm_setzero_ps();
218 /* Start inner kernel loop */
219 for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
222 /* Get j neighbor index, and coordinate index */
227 j_coord_offsetA = DIM*jnrA;
228 j_coord_offsetB = DIM*jnrB;
229 j_coord_offsetC = DIM*jnrC;
230 j_coord_offsetD = DIM*jnrD;
232 /* load j atom coordinates */
233 gmx_mm_load_3rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
234 x+j_coord_offsetC,x+j_coord_offsetD,
235 &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,&jy2,&jz2);
237 /* Calculate displacement vector */
238 dx00 = _mm_sub_ps(ix0,jx0);
239 dy00 = _mm_sub_ps(iy0,jy0);
240 dz00 = _mm_sub_ps(iz0,jz0);
241 dx01 = _mm_sub_ps(ix0,jx1);
242 dy01 = _mm_sub_ps(iy0,jy1);
243 dz01 = _mm_sub_ps(iz0,jz1);
244 dx02 = _mm_sub_ps(ix0,jx2);
245 dy02 = _mm_sub_ps(iy0,jy2);
246 dz02 = _mm_sub_ps(iz0,jz2);
247 dx10 = _mm_sub_ps(ix1,jx0);
248 dy10 = _mm_sub_ps(iy1,jy0);
249 dz10 = _mm_sub_ps(iz1,jz0);
250 dx11 = _mm_sub_ps(ix1,jx1);
251 dy11 = _mm_sub_ps(iy1,jy1);
252 dz11 = _mm_sub_ps(iz1,jz1);
253 dx12 = _mm_sub_ps(ix1,jx2);
254 dy12 = _mm_sub_ps(iy1,jy2);
255 dz12 = _mm_sub_ps(iz1,jz2);
256 dx20 = _mm_sub_ps(ix2,jx0);
257 dy20 = _mm_sub_ps(iy2,jy0);
258 dz20 = _mm_sub_ps(iz2,jz0);
259 dx21 = _mm_sub_ps(ix2,jx1);
260 dy21 = _mm_sub_ps(iy2,jy1);
261 dz21 = _mm_sub_ps(iz2,jz1);
262 dx22 = _mm_sub_ps(ix2,jx2);
263 dy22 = _mm_sub_ps(iy2,jy2);
264 dz22 = _mm_sub_ps(iz2,jz2);
266 /* Calculate squared distance and things based on it */
267 rsq00 = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
268 rsq01 = gmx_mm_calc_rsq_ps(dx01,dy01,dz01);
269 rsq02 = gmx_mm_calc_rsq_ps(dx02,dy02,dz02);
270 rsq10 = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
271 rsq11 = gmx_mm_calc_rsq_ps(dx11,dy11,dz11);
272 rsq12 = gmx_mm_calc_rsq_ps(dx12,dy12,dz12);
273 rsq20 = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
274 rsq21 = gmx_mm_calc_rsq_ps(dx21,dy21,dz21);
275 rsq22 = gmx_mm_calc_rsq_ps(dx22,dy22,dz22);
277 rinv00 = gmx_mm_invsqrt_ps(rsq00);
278 rinv01 = gmx_mm_invsqrt_ps(rsq01);
279 rinv02 = gmx_mm_invsqrt_ps(rsq02);
280 rinv10 = gmx_mm_invsqrt_ps(rsq10);
281 rinv11 = gmx_mm_invsqrt_ps(rsq11);
282 rinv12 = gmx_mm_invsqrt_ps(rsq12);
283 rinv20 = gmx_mm_invsqrt_ps(rsq20);
284 rinv21 = gmx_mm_invsqrt_ps(rsq21);
285 rinv22 = gmx_mm_invsqrt_ps(rsq22);
287 rinvsq00 = _mm_mul_ps(rinv00,rinv00);
288 rinvsq01 = _mm_mul_ps(rinv01,rinv01);
289 rinvsq02 = _mm_mul_ps(rinv02,rinv02);
290 rinvsq10 = _mm_mul_ps(rinv10,rinv10);
291 rinvsq11 = _mm_mul_ps(rinv11,rinv11);
292 rinvsq12 = _mm_mul_ps(rinv12,rinv12);
293 rinvsq20 = _mm_mul_ps(rinv20,rinv20);
294 rinvsq21 = _mm_mul_ps(rinv21,rinv21);
295 rinvsq22 = _mm_mul_ps(rinv22,rinv22);
297 fjx0 = _mm_setzero_ps();
298 fjy0 = _mm_setzero_ps();
299 fjz0 = _mm_setzero_ps();
300 fjx1 = _mm_setzero_ps();
301 fjy1 = _mm_setzero_ps();
302 fjz1 = _mm_setzero_ps();
303 fjx2 = _mm_setzero_ps();
304 fjy2 = _mm_setzero_ps();
305 fjz2 = _mm_setzero_ps();
307 /**************************
308 * CALCULATE INTERACTIONS *
309 **************************/
311 if (gmx_mm_any_lt(rsq00,rcutoff2))
314 r00 = _mm_mul_ps(rsq00,rinv00);
316 /* EWALD ELECTROSTATICS */
318 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
319 ewrt = _mm_mul_ps(r00,ewtabscale);
320 ewitab = _mm_cvttps_epi32(ewrt);
321 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
322 ewitab = _mm_slli_epi32(ewitab,2);
323 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
324 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
325 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
326 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
327 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
328 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
329 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
330 velec = _mm_mul_ps(qq00,_mm_sub_ps(rinv00,velec));
331 felec = _mm_mul_ps(_mm_mul_ps(qq00,rinv00),_mm_sub_ps(rinvsq00,felec));
333 d = _mm_sub_ps(r00,rswitch);
334 d = _mm_max_ps(d,_mm_setzero_ps());
335 d2 = _mm_mul_ps(d,d);
336 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
338 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
340 /* Evaluate switch function */
341 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
342 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv00,_mm_mul_ps(velec,dsw)) );
343 velec = _mm_mul_ps(velec,sw);
344 cutoff_mask = _mm_cmplt_ps(rsq00,rcutoff2);
346 /* Update potential sum for this i atom from the interaction with this j atom. */
347 velec = _mm_and_ps(velec,cutoff_mask);
348 velecsum = _mm_add_ps(velecsum,velec);
352 fscal = _mm_and_ps(fscal,cutoff_mask);
354 /* Calculate temporary vectorial force */
355 tx = _mm_mul_ps(fscal,dx00);
356 ty = _mm_mul_ps(fscal,dy00);
357 tz = _mm_mul_ps(fscal,dz00);
359 /* Update vectorial force */
360 fix0 = _mm_add_ps(fix0,tx);
361 fiy0 = _mm_add_ps(fiy0,ty);
362 fiz0 = _mm_add_ps(fiz0,tz);
364 fjx0 = _mm_add_ps(fjx0,tx);
365 fjy0 = _mm_add_ps(fjy0,ty);
366 fjz0 = _mm_add_ps(fjz0,tz);
370 /**************************
371 * CALCULATE INTERACTIONS *
372 **************************/
374 if (gmx_mm_any_lt(rsq01,rcutoff2))
377 r01 = _mm_mul_ps(rsq01,rinv01);
379 /* EWALD ELECTROSTATICS */
381 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
382 ewrt = _mm_mul_ps(r01,ewtabscale);
383 ewitab = _mm_cvttps_epi32(ewrt);
384 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
385 ewitab = _mm_slli_epi32(ewitab,2);
386 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
387 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
388 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
389 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
390 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
391 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
392 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
393 velec = _mm_mul_ps(qq01,_mm_sub_ps(rinv01,velec));
394 felec = _mm_mul_ps(_mm_mul_ps(qq01,rinv01),_mm_sub_ps(rinvsq01,felec));
396 d = _mm_sub_ps(r01,rswitch);
397 d = _mm_max_ps(d,_mm_setzero_ps());
398 d2 = _mm_mul_ps(d,d);
399 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
401 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
403 /* Evaluate switch function */
404 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
405 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv01,_mm_mul_ps(velec,dsw)) );
406 velec = _mm_mul_ps(velec,sw);
407 cutoff_mask = _mm_cmplt_ps(rsq01,rcutoff2);
409 /* Update potential sum for this i atom from the interaction with this j atom. */
410 velec = _mm_and_ps(velec,cutoff_mask);
411 velecsum = _mm_add_ps(velecsum,velec);
415 fscal = _mm_and_ps(fscal,cutoff_mask);
417 /* Calculate temporary vectorial force */
418 tx = _mm_mul_ps(fscal,dx01);
419 ty = _mm_mul_ps(fscal,dy01);
420 tz = _mm_mul_ps(fscal,dz01);
422 /* Update vectorial force */
423 fix0 = _mm_add_ps(fix0,tx);
424 fiy0 = _mm_add_ps(fiy0,ty);
425 fiz0 = _mm_add_ps(fiz0,tz);
427 fjx1 = _mm_add_ps(fjx1,tx);
428 fjy1 = _mm_add_ps(fjy1,ty);
429 fjz1 = _mm_add_ps(fjz1,tz);
433 /**************************
434 * CALCULATE INTERACTIONS *
435 **************************/
437 if (gmx_mm_any_lt(rsq02,rcutoff2))
440 r02 = _mm_mul_ps(rsq02,rinv02);
442 /* EWALD ELECTROSTATICS */
444 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
445 ewrt = _mm_mul_ps(r02,ewtabscale);
446 ewitab = _mm_cvttps_epi32(ewrt);
447 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
448 ewitab = _mm_slli_epi32(ewitab,2);
449 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
450 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
451 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
452 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
453 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
454 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
455 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
456 velec = _mm_mul_ps(qq02,_mm_sub_ps(rinv02,velec));
457 felec = _mm_mul_ps(_mm_mul_ps(qq02,rinv02),_mm_sub_ps(rinvsq02,felec));
459 d = _mm_sub_ps(r02,rswitch);
460 d = _mm_max_ps(d,_mm_setzero_ps());
461 d2 = _mm_mul_ps(d,d);
462 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
464 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
466 /* Evaluate switch function */
467 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
468 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv02,_mm_mul_ps(velec,dsw)) );
469 velec = _mm_mul_ps(velec,sw);
470 cutoff_mask = _mm_cmplt_ps(rsq02,rcutoff2);
472 /* Update potential sum for this i atom from the interaction with this j atom. */
473 velec = _mm_and_ps(velec,cutoff_mask);
474 velecsum = _mm_add_ps(velecsum,velec);
478 fscal = _mm_and_ps(fscal,cutoff_mask);
480 /* Calculate temporary vectorial force */
481 tx = _mm_mul_ps(fscal,dx02);
482 ty = _mm_mul_ps(fscal,dy02);
483 tz = _mm_mul_ps(fscal,dz02);
485 /* Update vectorial force */
486 fix0 = _mm_add_ps(fix0,tx);
487 fiy0 = _mm_add_ps(fiy0,ty);
488 fiz0 = _mm_add_ps(fiz0,tz);
490 fjx2 = _mm_add_ps(fjx2,tx);
491 fjy2 = _mm_add_ps(fjy2,ty);
492 fjz2 = _mm_add_ps(fjz2,tz);
496 /**************************
497 * CALCULATE INTERACTIONS *
498 **************************/
500 if (gmx_mm_any_lt(rsq10,rcutoff2))
503 r10 = _mm_mul_ps(rsq10,rinv10);
505 /* EWALD ELECTROSTATICS */
507 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
508 ewrt = _mm_mul_ps(r10,ewtabscale);
509 ewitab = _mm_cvttps_epi32(ewrt);
510 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
511 ewitab = _mm_slli_epi32(ewitab,2);
512 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
513 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
514 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
515 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
516 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
517 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
518 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
519 velec = _mm_mul_ps(qq10,_mm_sub_ps(rinv10,velec));
520 felec = _mm_mul_ps(_mm_mul_ps(qq10,rinv10),_mm_sub_ps(rinvsq10,felec));
522 d = _mm_sub_ps(r10,rswitch);
523 d = _mm_max_ps(d,_mm_setzero_ps());
524 d2 = _mm_mul_ps(d,d);
525 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
527 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
529 /* Evaluate switch function */
530 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
531 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv10,_mm_mul_ps(velec,dsw)) );
532 velec = _mm_mul_ps(velec,sw);
533 cutoff_mask = _mm_cmplt_ps(rsq10,rcutoff2);
535 /* Update potential sum for this i atom from the interaction with this j atom. */
536 velec = _mm_and_ps(velec,cutoff_mask);
537 velecsum = _mm_add_ps(velecsum,velec);
541 fscal = _mm_and_ps(fscal,cutoff_mask);
543 /* Calculate temporary vectorial force */
544 tx = _mm_mul_ps(fscal,dx10);
545 ty = _mm_mul_ps(fscal,dy10);
546 tz = _mm_mul_ps(fscal,dz10);
548 /* Update vectorial force */
549 fix1 = _mm_add_ps(fix1,tx);
550 fiy1 = _mm_add_ps(fiy1,ty);
551 fiz1 = _mm_add_ps(fiz1,tz);
553 fjx0 = _mm_add_ps(fjx0,tx);
554 fjy0 = _mm_add_ps(fjy0,ty);
555 fjz0 = _mm_add_ps(fjz0,tz);
559 /**************************
560 * CALCULATE INTERACTIONS *
561 **************************/
563 if (gmx_mm_any_lt(rsq11,rcutoff2))
566 r11 = _mm_mul_ps(rsq11,rinv11);
568 /* EWALD ELECTROSTATICS */
570 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
571 ewrt = _mm_mul_ps(r11,ewtabscale);
572 ewitab = _mm_cvttps_epi32(ewrt);
573 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
574 ewitab = _mm_slli_epi32(ewitab,2);
575 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
576 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
577 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
578 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
579 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
580 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
581 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
582 velec = _mm_mul_ps(qq11,_mm_sub_ps(rinv11,velec));
583 felec = _mm_mul_ps(_mm_mul_ps(qq11,rinv11),_mm_sub_ps(rinvsq11,felec));
585 d = _mm_sub_ps(r11,rswitch);
586 d = _mm_max_ps(d,_mm_setzero_ps());
587 d2 = _mm_mul_ps(d,d);
588 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
590 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
592 /* Evaluate switch function */
593 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
594 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv11,_mm_mul_ps(velec,dsw)) );
595 velec = _mm_mul_ps(velec,sw);
596 cutoff_mask = _mm_cmplt_ps(rsq11,rcutoff2);
598 /* Update potential sum for this i atom from the interaction with this j atom. */
599 velec = _mm_and_ps(velec,cutoff_mask);
600 velecsum = _mm_add_ps(velecsum,velec);
604 fscal = _mm_and_ps(fscal,cutoff_mask);
606 /* Calculate temporary vectorial force */
607 tx = _mm_mul_ps(fscal,dx11);
608 ty = _mm_mul_ps(fscal,dy11);
609 tz = _mm_mul_ps(fscal,dz11);
611 /* Update vectorial force */
612 fix1 = _mm_add_ps(fix1,tx);
613 fiy1 = _mm_add_ps(fiy1,ty);
614 fiz1 = _mm_add_ps(fiz1,tz);
616 fjx1 = _mm_add_ps(fjx1,tx);
617 fjy1 = _mm_add_ps(fjy1,ty);
618 fjz1 = _mm_add_ps(fjz1,tz);
622 /**************************
623 * CALCULATE INTERACTIONS *
624 **************************/
626 if (gmx_mm_any_lt(rsq12,rcutoff2))
629 r12 = _mm_mul_ps(rsq12,rinv12);
631 /* EWALD ELECTROSTATICS */
633 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
634 ewrt = _mm_mul_ps(r12,ewtabscale);
635 ewitab = _mm_cvttps_epi32(ewrt);
636 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
637 ewitab = _mm_slli_epi32(ewitab,2);
638 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
639 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
640 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
641 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
642 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
643 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
644 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
645 velec = _mm_mul_ps(qq12,_mm_sub_ps(rinv12,velec));
646 felec = _mm_mul_ps(_mm_mul_ps(qq12,rinv12),_mm_sub_ps(rinvsq12,felec));
648 d = _mm_sub_ps(r12,rswitch);
649 d = _mm_max_ps(d,_mm_setzero_ps());
650 d2 = _mm_mul_ps(d,d);
651 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
653 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
655 /* Evaluate switch function */
656 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
657 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv12,_mm_mul_ps(velec,dsw)) );
658 velec = _mm_mul_ps(velec,sw);
659 cutoff_mask = _mm_cmplt_ps(rsq12,rcutoff2);
661 /* Update potential sum for this i atom from the interaction with this j atom. */
662 velec = _mm_and_ps(velec,cutoff_mask);
663 velecsum = _mm_add_ps(velecsum,velec);
667 fscal = _mm_and_ps(fscal,cutoff_mask);
669 /* Calculate temporary vectorial force */
670 tx = _mm_mul_ps(fscal,dx12);
671 ty = _mm_mul_ps(fscal,dy12);
672 tz = _mm_mul_ps(fscal,dz12);
674 /* Update vectorial force */
675 fix1 = _mm_add_ps(fix1,tx);
676 fiy1 = _mm_add_ps(fiy1,ty);
677 fiz1 = _mm_add_ps(fiz1,tz);
679 fjx2 = _mm_add_ps(fjx2,tx);
680 fjy2 = _mm_add_ps(fjy2,ty);
681 fjz2 = _mm_add_ps(fjz2,tz);
685 /**************************
686 * CALCULATE INTERACTIONS *
687 **************************/
689 if (gmx_mm_any_lt(rsq20,rcutoff2))
692 r20 = _mm_mul_ps(rsq20,rinv20);
694 /* EWALD ELECTROSTATICS */
696 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
697 ewrt = _mm_mul_ps(r20,ewtabscale);
698 ewitab = _mm_cvttps_epi32(ewrt);
699 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
700 ewitab = _mm_slli_epi32(ewitab,2);
701 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
702 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
703 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
704 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
705 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
706 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
707 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
708 velec = _mm_mul_ps(qq20,_mm_sub_ps(rinv20,velec));
709 felec = _mm_mul_ps(_mm_mul_ps(qq20,rinv20),_mm_sub_ps(rinvsq20,felec));
711 d = _mm_sub_ps(r20,rswitch);
712 d = _mm_max_ps(d,_mm_setzero_ps());
713 d2 = _mm_mul_ps(d,d);
714 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
716 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
718 /* Evaluate switch function */
719 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
720 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv20,_mm_mul_ps(velec,dsw)) );
721 velec = _mm_mul_ps(velec,sw);
722 cutoff_mask = _mm_cmplt_ps(rsq20,rcutoff2);
724 /* Update potential sum for this i atom from the interaction with this j atom. */
725 velec = _mm_and_ps(velec,cutoff_mask);
726 velecsum = _mm_add_ps(velecsum,velec);
730 fscal = _mm_and_ps(fscal,cutoff_mask);
732 /* Calculate temporary vectorial force */
733 tx = _mm_mul_ps(fscal,dx20);
734 ty = _mm_mul_ps(fscal,dy20);
735 tz = _mm_mul_ps(fscal,dz20);
737 /* Update vectorial force */
738 fix2 = _mm_add_ps(fix2,tx);
739 fiy2 = _mm_add_ps(fiy2,ty);
740 fiz2 = _mm_add_ps(fiz2,tz);
742 fjx0 = _mm_add_ps(fjx0,tx);
743 fjy0 = _mm_add_ps(fjy0,ty);
744 fjz0 = _mm_add_ps(fjz0,tz);
748 /**************************
749 * CALCULATE INTERACTIONS *
750 **************************/
752 if (gmx_mm_any_lt(rsq21,rcutoff2))
755 r21 = _mm_mul_ps(rsq21,rinv21);
757 /* EWALD ELECTROSTATICS */
759 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
760 ewrt = _mm_mul_ps(r21,ewtabscale);
761 ewitab = _mm_cvttps_epi32(ewrt);
762 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
763 ewitab = _mm_slli_epi32(ewitab,2);
764 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
765 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
766 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
767 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
768 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
769 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
770 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
771 velec = _mm_mul_ps(qq21,_mm_sub_ps(rinv21,velec));
772 felec = _mm_mul_ps(_mm_mul_ps(qq21,rinv21),_mm_sub_ps(rinvsq21,felec));
774 d = _mm_sub_ps(r21,rswitch);
775 d = _mm_max_ps(d,_mm_setzero_ps());
776 d2 = _mm_mul_ps(d,d);
777 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
779 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
781 /* Evaluate switch function */
782 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
783 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv21,_mm_mul_ps(velec,dsw)) );
784 velec = _mm_mul_ps(velec,sw);
785 cutoff_mask = _mm_cmplt_ps(rsq21,rcutoff2);
787 /* Update potential sum for this i atom from the interaction with this j atom. */
788 velec = _mm_and_ps(velec,cutoff_mask);
789 velecsum = _mm_add_ps(velecsum,velec);
793 fscal = _mm_and_ps(fscal,cutoff_mask);
795 /* Calculate temporary vectorial force */
796 tx = _mm_mul_ps(fscal,dx21);
797 ty = _mm_mul_ps(fscal,dy21);
798 tz = _mm_mul_ps(fscal,dz21);
800 /* Update vectorial force */
801 fix2 = _mm_add_ps(fix2,tx);
802 fiy2 = _mm_add_ps(fiy2,ty);
803 fiz2 = _mm_add_ps(fiz2,tz);
805 fjx1 = _mm_add_ps(fjx1,tx);
806 fjy1 = _mm_add_ps(fjy1,ty);
807 fjz1 = _mm_add_ps(fjz1,tz);
811 /**************************
812 * CALCULATE INTERACTIONS *
813 **************************/
815 if (gmx_mm_any_lt(rsq22,rcutoff2))
818 r22 = _mm_mul_ps(rsq22,rinv22);
820 /* EWALD ELECTROSTATICS */
822 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
823 ewrt = _mm_mul_ps(r22,ewtabscale);
824 ewitab = _mm_cvttps_epi32(ewrt);
825 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
826 ewitab = _mm_slli_epi32(ewitab,2);
827 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
828 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
829 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
830 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
831 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
832 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
833 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
834 velec = _mm_mul_ps(qq22,_mm_sub_ps(rinv22,velec));
835 felec = _mm_mul_ps(_mm_mul_ps(qq22,rinv22),_mm_sub_ps(rinvsq22,felec));
837 d = _mm_sub_ps(r22,rswitch);
838 d = _mm_max_ps(d,_mm_setzero_ps());
839 d2 = _mm_mul_ps(d,d);
840 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
842 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
844 /* Evaluate switch function */
845 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
846 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv22,_mm_mul_ps(velec,dsw)) );
847 velec = _mm_mul_ps(velec,sw);
848 cutoff_mask = _mm_cmplt_ps(rsq22,rcutoff2);
850 /* Update potential sum for this i atom from the interaction with this j atom. */
851 velec = _mm_and_ps(velec,cutoff_mask);
852 velecsum = _mm_add_ps(velecsum,velec);
856 fscal = _mm_and_ps(fscal,cutoff_mask);
858 /* Calculate temporary vectorial force */
859 tx = _mm_mul_ps(fscal,dx22);
860 ty = _mm_mul_ps(fscal,dy22);
861 tz = _mm_mul_ps(fscal,dz22);
863 /* Update vectorial force */
864 fix2 = _mm_add_ps(fix2,tx);
865 fiy2 = _mm_add_ps(fiy2,ty);
866 fiz2 = _mm_add_ps(fiz2,tz);
868 fjx2 = _mm_add_ps(fjx2,tx);
869 fjy2 = _mm_add_ps(fjy2,ty);
870 fjz2 = _mm_add_ps(fjz2,tz);
874 fjptrA = f+j_coord_offsetA;
875 fjptrB = f+j_coord_offsetB;
876 fjptrC = f+j_coord_offsetC;
877 fjptrD = f+j_coord_offsetD;
879 gmx_mm_decrement_3rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,
880 fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,fjx2,fjy2,fjz2);
882 /* Inner loop uses 585 flops */
888 /* Get j neighbor index, and coordinate index */
889 jnrlistA = jjnr[jidx];
890 jnrlistB = jjnr[jidx+1];
891 jnrlistC = jjnr[jidx+2];
892 jnrlistD = jjnr[jidx+3];
893 /* Sign of each element will be negative for non-real atoms.
894 * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
895 * so use it as val = _mm_andnot_ps(mask,val) to clear dummy entries.
897 dummy_mask = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
898 jnrA = (jnrlistA>=0) ? jnrlistA : 0;
899 jnrB = (jnrlistB>=0) ? jnrlistB : 0;
900 jnrC = (jnrlistC>=0) ? jnrlistC : 0;
901 jnrD = (jnrlistD>=0) ? jnrlistD : 0;
902 j_coord_offsetA = DIM*jnrA;
903 j_coord_offsetB = DIM*jnrB;
904 j_coord_offsetC = DIM*jnrC;
905 j_coord_offsetD = DIM*jnrD;
907 /* load j atom coordinates */
908 gmx_mm_load_3rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
909 x+j_coord_offsetC,x+j_coord_offsetD,
910 &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,&jy2,&jz2);
912 /* Calculate displacement vector */
913 dx00 = _mm_sub_ps(ix0,jx0);
914 dy00 = _mm_sub_ps(iy0,jy0);
915 dz00 = _mm_sub_ps(iz0,jz0);
916 dx01 = _mm_sub_ps(ix0,jx1);
917 dy01 = _mm_sub_ps(iy0,jy1);
918 dz01 = _mm_sub_ps(iz0,jz1);
919 dx02 = _mm_sub_ps(ix0,jx2);
920 dy02 = _mm_sub_ps(iy0,jy2);
921 dz02 = _mm_sub_ps(iz0,jz2);
922 dx10 = _mm_sub_ps(ix1,jx0);
923 dy10 = _mm_sub_ps(iy1,jy0);
924 dz10 = _mm_sub_ps(iz1,jz0);
925 dx11 = _mm_sub_ps(ix1,jx1);
926 dy11 = _mm_sub_ps(iy1,jy1);
927 dz11 = _mm_sub_ps(iz1,jz1);
928 dx12 = _mm_sub_ps(ix1,jx2);
929 dy12 = _mm_sub_ps(iy1,jy2);
930 dz12 = _mm_sub_ps(iz1,jz2);
931 dx20 = _mm_sub_ps(ix2,jx0);
932 dy20 = _mm_sub_ps(iy2,jy0);
933 dz20 = _mm_sub_ps(iz2,jz0);
934 dx21 = _mm_sub_ps(ix2,jx1);
935 dy21 = _mm_sub_ps(iy2,jy1);
936 dz21 = _mm_sub_ps(iz2,jz1);
937 dx22 = _mm_sub_ps(ix2,jx2);
938 dy22 = _mm_sub_ps(iy2,jy2);
939 dz22 = _mm_sub_ps(iz2,jz2);
941 /* Calculate squared distance and things based on it */
942 rsq00 = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
943 rsq01 = gmx_mm_calc_rsq_ps(dx01,dy01,dz01);
944 rsq02 = gmx_mm_calc_rsq_ps(dx02,dy02,dz02);
945 rsq10 = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
946 rsq11 = gmx_mm_calc_rsq_ps(dx11,dy11,dz11);
947 rsq12 = gmx_mm_calc_rsq_ps(dx12,dy12,dz12);
948 rsq20 = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
949 rsq21 = gmx_mm_calc_rsq_ps(dx21,dy21,dz21);
950 rsq22 = gmx_mm_calc_rsq_ps(dx22,dy22,dz22);
952 rinv00 = gmx_mm_invsqrt_ps(rsq00);
953 rinv01 = gmx_mm_invsqrt_ps(rsq01);
954 rinv02 = gmx_mm_invsqrt_ps(rsq02);
955 rinv10 = gmx_mm_invsqrt_ps(rsq10);
956 rinv11 = gmx_mm_invsqrt_ps(rsq11);
957 rinv12 = gmx_mm_invsqrt_ps(rsq12);
958 rinv20 = gmx_mm_invsqrt_ps(rsq20);
959 rinv21 = gmx_mm_invsqrt_ps(rsq21);
960 rinv22 = gmx_mm_invsqrt_ps(rsq22);
962 rinvsq00 = _mm_mul_ps(rinv00,rinv00);
963 rinvsq01 = _mm_mul_ps(rinv01,rinv01);
964 rinvsq02 = _mm_mul_ps(rinv02,rinv02);
965 rinvsq10 = _mm_mul_ps(rinv10,rinv10);
966 rinvsq11 = _mm_mul_ps(rinv11,rinv11);
967 rinvsq12 = _mm_mul_ps(rinv12,rinv12);
968 rinvsq20 = _mm_mul_ps(rinv20,rinv20);
969 rinvsq21 = _mm_mul_ps(rinv21,rinv21);
970 rinvsq22 = _mm_mul_ps(rinv22,rinv22);
972 fjx0 = _mm_setzero_ps();
973 fjy0 = _mm_setzero_ps();
974 fjz0 = _mm_setzero_ps();
975 fjx1 = _mm_setzero_ps();
976 fjy1 = _mm_setzero_ps();
977 fjz1 = _mm_setzero_ps();
978 fjx2 = _mm_setzero_ps();
979 fjy2 = _mm_setzero_ps();
980 fjz2 = _mm_setzero_ps();
982 /**************************
983 * CALCULATE INTERACTIONS *
984 **************************/
986 if (gmx_mm_any_lt(rsq00,rcutoff2))
989 r00 = _mm_mul_ps(rsq00,rinv00);
990 r00 = _mm_andnot_ps(dummy_mask,r00);
992 /* EWALD ELECTROSTATICS */
994 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
995 ewrt = _mm_mul_ps(r00,ewtabscale);
996 ewitab = _mm_cvttps_epi32(ewrt);
997 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
998 ewitab = _mm_slli_epi32(ewitab,2);
999 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1000 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
1001 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
1002 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
1003 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
1004 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
1005 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
1006 velec = _mm_mul_ps(qq00,_mm_sub_ps(rinv00,velec));
1007 felec = _mm_mul_ps(_mm_mul_ps(qq00,rinv00),_mm_sub_ps(rinvsq00,felec));
1009 d = _mm_sub_ps(r00,rswitch);
1010 d = _mm_max_ps(d,_mm_setzero_ps());
1011 d2 = _mm_mul_ps(d,d);
1012 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
1014 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
1016 /* Evaluate switch function */
1017 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1018 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv00,_mm_mul_ps(velec,dsw)) );
1019 velec = _mm_mul_ps(velec,sw);
1020 cutoff_mask = _mm_cmplt_ps(rsq00,rcutoff2);
1022 /* Update potential sum for this i atom from the interaction with this j atom. */
1023 velec = _mm_and_ps(velec,cutoff_mask);
1024 velec = _mm_andnot_ps(dummy_mask,velec);
1025 velecsum = _mm_add_ps(velecsum,velec);
1029 fscal = _mm_and_ps(fscal,cutoff_mask);
1031 fscal = _mm_andnot_ps(dummy_mask,fscal);
1033 /* Calculate temporary vectorial force */
1034 tx = _mm_mul_ps(fscal,dx00);
1035 ty = _mm_mul_ps(fscal,dy00);
1036 tz = _mm_mul_ps(fscal,dz00);
1038 /* Update vectorial force */
1039 fix0 = _mm_add_ps(fix0,tx);
1040 fiy0 = _mm_add_ps(fiy0,ty);
1041 fiz0 = _mm_add_ps(fiz0,tz);
1043 fjx0 = _mm_add_ps(fjx0,tx);
1044 fjy0 = _mm_add_ps(fjy0,ty);
1045 fjz0 = _mm_add_ps(fjz0,tz);
1049 /**************************
1050 * CALCULATE INTERACTIONS *
1051 **************************/
1053 if (gmx_mm_any_lt(rsq01,rcutoff2))
1056 r01 = _mm_mul_ps(rsq01,rinv01);
1057 r01 = _mm_andnot_ps(dummy_mask,r01);
1059 /* EWALD ELECTROSTATICS */
1061 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1062 ewrt = _mm_mul_ps(r01,ewtabscale);
1063 ewitab = _mm_cvttps_epi32(ewrt);
1064 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
1065 ewitab = _mm_slli_epi32(ewitab,2);
1066 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1067 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
1068 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
1069 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
1070 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
1071 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
1072 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
1073 velec = _mm_mul_ps(qq01,_mm_sub_ps(rinv01,velec));
1074 felec = _mm_mul_ps(_mm_mul_ps(qq01,rinv01),_mm_sub_ps(rinvsq01,felec));
1076 d = _mm_sub_ps(r01,rswitch);
1077 d = _mm_max_ps(d,_mm_setzero_ps());
1078 d2 = _mm_mul_ps(d,d);
1079 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
1081 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
1083 /* Evaluate switch function */
1084 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1085 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv01,_mm_mul_ps(velec,dsw)) );
1086 velec = _mm_mul_ps(velec,sw);
1087 cutoff_mask = _mm_cmplt_ps(rsq01,rcutoff2);
1089 /* Update potential sum for this i atom from the interaction with this j atom. */
1090 velec = _mm_and_ps(velec,cutoff_mask);
1091 velec = _mm_andnot_ps(dummy_mask,velec);
1092 velecsum = _mm_add_ps(velecsum,velec);
1096 fscal = _mm_and_ps(fscal,cutoff_mask);
1098 fscal = _mm_andnot_ps(dummy_mask,fscal);
1100 /* Calculate temporary vectorial force */
1101 tx = _mm_mul_ps(fscal,dx01);
1102 ty = _mm_mul_ps(fscal,dy01);
1103 tz = _mm_mul_ps(fscal,dz01);
1105 /* Update vectorial force */
1106 fix0 = _mm_add_ps(fix0,tx);
1107 fiy0 = _mm_add_ps(fiy0,ty);
1108 fiz0 = _mm_add_ps(fiz0,tz);
1110 fjx1 = _mm_add_ps(fjx1,tx);
1111 fjy1 = _mm_add_ps(fjy1,ty);
1112 fjz1 = _mm_add_ps(fjz1,tz);
1116 /**************************
1117 * CALCULATE INTERACTIONS *
1118 **************************/
1120 if (gmx_mm_any_lt(rsq02,rcutoff2))
1123 r02 = _mm_mul_ps(rsq02,rinv02);
1124 r02 = _mm_andnot_ps(dummy_mask,r02);
1126 /* EWALD ELECTROSTATICS */
1128 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1129 ewrt = _mm_mul_ps(r02,ewtabscale);
1130 ewitab = _mm_cvttps_epi32(ewrt);
1131 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
1132 ewitab = _mm_slli_epi32(ewitab,2);
1133 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1134 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
1135 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
1136 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
1137 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
1138 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
1139 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
1140 velec = _mm_mul_ps(qq02,_mm_sub_ps(rinv02,velec));
1141 felec = _mm_mul_ps(_mm_mul_ps(qq02,rinv02),_mm_sub_ps(rinvsq02,felec));
1143 d = _mm_sub_ps(r02,rswitch);
1144 d = _mm_max_ps(d,_mm_setzero_ps());
1145 d2 = _mm_mul_ps(d,d);
1146 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
1148 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
1150 /* Evaluate switch function */
1151 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1152 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv02,_mm_mul_ps(velec,dsw)) );
1153 velec = _mm_mul_ps(velec,sw);
1154 cutoff_mask = _mm_cmplt_ps(rsq02,rcutoff2);
1156 /* Update potential sum for this i atom from the interaction with this j atom. */
1157 velec = _mm_and_ps(velec,cutoff_mask);
1158 velec = _mm_andnot_ps(dummy_mask,velec);
1159 velecsum = _mm_add_ps(velecsum,velec);
1163 fscal = _mm_and_ps(fscal,cutoff_mask);
1165 fscal = _mm_andnot_ps(dummy_mask,fscal);
1167 /* Calculate temporary vectorial force */
1168 tx = _mm_mul_ps(fscal,dx02);
1169 ty = _mm_mul_ps(fscal,dy02);
1170 tz = _mm_mul_ps(fscal,dz02);
1172 /* Update vectorial force */
1173 fix0 = _mm_add_ps(fix0,tx);
1174 fiy0 = _mm_add_ps(fiy0,ty);
1175 fiz0 = _mm_add_ps(fiz0,tz);
1177 fjx2 = _mm_add_ps(fjx2,tx);
1178 fjy2 = _mm_add_ps(fjy2,ty);
1179 fjz2 = _mm_add_ps(fjz2,tz);
1183 /**************************
1184 * CALCULATE INTERACTIONS *
1185 **************************/
1187 if (gmx_mm_any_lt(rsq10,rcutoff2))
1190 r10 = _mm_mul_ps(rsq10,rinv10);
1191 r10 = _mm_andnot_ps(dummy_mask,r10);
1193 /* EWALD ELECTROSTATICS */
1195 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1196 ewrt = _mm_mul_ps(r10,ewtabscale);
1197 ewitab = _mm_cvttps_epi32(ewrt);
1198 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
1199 ewitab = _mm_slli_epi32(ewitab,2);
1200 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1201 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
1202 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
1203 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
1204 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
1205 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
1206 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
1207 velec = _mm_mul_ps(qq10,_mm_sub_ps(rinv10,velec));
1208 felec = _mm_mul_ps(_mm_mul_ps(qq10,rinv10),_mm_sub_ps(rinvsq10,felec));
1210 d = _mm_sub_ps(r10,rswitch);
1211 d = _mm_max_ps(d,_mm_setzero_ps());
1212 d2 = _mm_mul_ps(d,d);
1213 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
1215 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
1217 /* Evaluate switch function */
1218 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1219 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv10,_mm_mul_ps(velec,dsw)) );
1220 velec = _mm_mul_ps(velec,sw);
1221 cutoff_mask = _mm_cmplt_ps(rsq10,rcutoff2);
1223 /* Update potential sum for this i atom from the interaction with this j atom. */
1224 velec = _mm_and_ps(velec,cutoff_mask);
1225 velec = _mm_andnot_ps(dummy_mask,velec);
1226 velecsum = _mm_add_ps(velecsum,velec);
1230 fscal = _mm_and_ps(fscal,cutoff_mask);
1232 fscal = _mm_andnot_ps(dummy_mask,fscal);
1234 /* Calculate temporary vectorial force */
1235 tx = _mm_mul_ps(fscal,dx10);
1236 ty = _mm_mul_ps(fscal,dy10);
1237 tz = _mm_mul_ps(fscal,dz10);
1239 /* Update vectorial force */
1240 fix1 = _mm_add_ps(fix1,tx);
1241 fiy1 = _mm_add_ps(fiy1,ty);
1242 fiz1 = _mm_add_ps(fiz1,tz);
1244 fjx0 = _mm_add_ps(fjx0,tx);
1245 fjy0 = _mm_add_ps(fjy0,ty);
1246 fjz0 = _mm_add_ps(fjz0,tz);
1250 /**************************
1251 * CALCULATE INTERACTIONS *
1252 **************************/
1254 if (gmx_mm_any_lt(rsq11,rcutoff2))
1257 r11 = _mm_mul_ps(rsq11,rinv11);
1258 r11 = _mm_andnot_ps(dummy_mask,r11);
1260 /* EWALD ELECTROSTATICS */
1262 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1263 ewrt = _mm_mul_ps(r11,ewtabscale);
1264 ewitab = _mm_cvttps_epi32(ewrt);
1265 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
1266 ewitab = _mm_slli_epi32(ewitab,2);
1267 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1268 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
1269 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
1270 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
1271 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
1272 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
1273 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
1274 velec = _mm_mul_ps(qq11,_mm_sub_ps(rinv11,velec));
1275 felec = _mm_mul_ps(_mm_mul_ps(qq11,rinv11),_mm_sub_ps(rinvsq11,felec));
1277 d = _mm_sub_ps(r11,rswitch);
1278 d = _mm_max_ps(d,_mm_setzero_ps());
1279 d2 = _mm_mul_ps(d,d);
1280 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
1282 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
1284 /* Evaluate switch function */
1285 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1286 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv11,_mm_mul_ps(velec,dsw)) );
1287 velec = _mm_mul_ps(velec,sw);
1288 cutoff_mask = _mm_cmplt_ps(rsq11,rcutoff2);
1290 /* Update potential sum for this i atom from the interaction with this j atom. */
1291 velec = _mm_and_ps(velec,cutoff_mask);
1292 velec = _mm_andnot_ps(dummy_mask,velec);
1293 velecsum = _mm_add_ps(velecsum,velec);
1297 fscal = _mm_and_ps(fscal,cutoff_mask);
1299 fscal = _mm_andnot_ps(dummy_mask,fscal);
1301 /* Calculate temporary vectorial force */
1302 tx = _mm_mul_ps(fscal,dx11);
1303 ty = _mm_mul_ps(fscal,dy11);
1304 tz = _mm_mul_ps(fscal,dz11);
1306 /* Update vectorial force */
1307 fix1 = _mm_add_ps(fix1,tx);
1308 fiy1 = _mm_add_ps(fiy1,ty);
1309 fiz1 = _mm_add_ps(fiz1,tz);
1311 fjx1 = _mm_add_ps(fjx1,tx);
1312 fjy1 = _mm_add_ps(fjy1,ty);
1313 fjz1 = _mm_add_ps(fjz1,tz);
1317 /**************************
1318 * CALCULATE INTERACTIONS *
1319 **************************/
1321 if (gmx_mm_any_lt(rsq12,rcutoff2))
1324 r12 = _mm_mul_ps(rsq12,rinv12);
1325 r12 = _mm_andnot_ps(dummy_mask,r12);
1327 /* EWALD ELECTROSTATICS */
1329 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1330 ewrt = _mm_mul_ps(r12,ewtabscale);
1331 ewitab = _mm_cvttps_epi32(ewrt);
1332 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
1333 ewitab = _mm_slli_epi32(ewitab,2);
1334 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1335 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
1336 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
1337 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
1338 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
1339 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
1340 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
1341 velec = _mm_mul_ps(qq12,_mm_sub_ps(rinv12,velec));
1342 felec = _mm_mul_ps(_mm_mul_ps(qq12,rinv12),_mm_sub_ps(rinvsq12,felec));
1344 d = _mm_sub_ps(r12,rswitch);
1345 d = _mm_max_ps(d,_mm_setzero_ps());
1346 d2 = _mm_mul_ps(d,d);
1347 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
1349 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
1351 /* Evaluate switch function */
1352 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1353 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv12,_mm_mul_ps(velec,dsw)) );
1354 velec = _mm_mul_ps(velec,sw);
1355 cutoff_mask = _mm_cmplt_ps(rsq12,rcutoff2);
1357 /* Update potential sum for this i atom from the interaction with this j atom. */
1358 velec = _mm_and_ps(velec,cutoff_mask);
1359 velec = _mm_andnot_ps(dummy_mask,velec);
1360 velecsum = _mm_add_ps(velecsum,velec);
1364 fscal = _mm_and_ps(fscal,cutoff_mask);
1366 fscal = _mm_andnot_ps(dummy_mask,fscal);
1368 /* Calculate temporary vectorial force */
1369 tx = _mm_mul_ps(fscal,dx12);
1370 ty = _mm_mul_ps(fscal,dy12);
1371 tz = _mm_mul_ps(fscal,dz12);
1373 /* Update vectorial force */
1374 fix1 = _mm_add_ps(fix1,tx);
1375 fiy1 = _mm_add_ps(fiy1,ty);
1376 fiz1 = _mm_add_ps(fiz1,tz);
1378 fjx2 = _mm_add_ps(fjx2,tx);
1379 fjy2 = _mm_add_ps(fjy2,ty);
1380 fjz2 = _mm_add_ps(fjz2,tz);
1384 /**************************
1385 * CALCULATE INTERACTIONS *
1386 **************************/
1388 if (gmx_mm_any_lt(rsq20,rcutoff2))
1391 r20 = _mm_mul_ps(rsq20,rinv20);
1392 r20 = _mm_andnot_ps(dummy_mask,r20);
1394 /* EWALD ELECTROSTATICS */
1396 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1397 ewrt = _mm_mul_ps(r20,ewtabscale);
1398 ewitab = _mm_cvttps_epi32(ewrt);
1399 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
1400 ewitab = _mm_slli_epi32(ewitab,2);
1401 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1402 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
1403 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
1404 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
1405 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
1406 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
1407 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
1408 velec = _mm_mul_ps(qq20,_mm_sub_ps(rinv20,velec));
1409 felec = _mm_mul_ps(_mm_mul_ps(qq20,rinv20),_mm_sub_ps(rinvsq20,felec));
1411 d = _mm_sub_ps(r20,rswitch);
1412 d = _mm_max_ps(d,_mm_setzero_ps());
1413 d2 = _mm_mul_ps(d,d);
1414 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
1416 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
1418 /* Evaluate switch function */
1419 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1420 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv20,_mm_mul_ps(velec,dsw)) );
1421 velec = _mm_mul_ps(velec,sw);
1422 cutoff_mask = _mm_cmplt_ps(rsq20,rcutoff2);
1424 /* Update potential sum for this i atom from the interaction with this j atom. */
1425 velec = _mm_and_ps(velec,cutoff_mask);
1426 velec = _mm_andnot_ps(dummy_mask,velec);
1427 velecsum = _mm_add_ps(velecsum,velec);
1431 fscal = _mm_and_ps(fscal,cutoff_mask);
1433 fscal = _mm_andnot_ps(dummy_mask,fscal);
1435 /* Calculate temporary vectorial force */
1436 tx = _mm_mul_ps(fscal,dx20);
1437 ty = _mm_mul_ps(fscal,dy20);
1438 tz = _mm_mul_ps(fscal,dz20);
1440 /* Update vectorial force */
1441 fix2 = _mm_add_ps(fix2,tx);
1442 fiy2 = _mm_add_ps(fiy2,ty);
1443 fiz2 = _mm_add_ps(fiz2,tz);
1445 fjx0 = _mm_add_ps(fjx0,tx);
1446 fjy0 = _mm_add_ps(fjy0,ty);
1447 fjz0 = _mm_add_ps(fjz0,tz);
1451 /**************************
1452 * CALCULATE INTERACTIONS *
1453 **************************/
1455 if (gmx_mm_any_lt(rsq21,rcutoff2))
1458 r21 = _mm_mul_ps(rsq21,rinv21);
1459 r21 = _mm_andnot_ps(dummy_mask,r21);
1461 /* EWALD ELECTROSTATICS */
1463 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1464 ewrt = _mm_mul_ps(r21,ewtabscale);
1465 ewitab = _mm_cvttps_epi32(ewrt);
1466 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
1467 ewitab = _mm_slli_epi32(ewitab,2);
1468 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1469 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
1470 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
1471 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
1472 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
1473 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
1474 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
1475 velec = _mm_mul_ps(qq21,_mm_sub_ps(rinv21,velec));
1476 felec = _mm_mul_ps(_mm_mul_ps(qq21,rinv21),_mm_sub_ps(rinvsq21,felec));
1478 d = _mm_sub_ps(r21,rswitch);
1479 d = _mm_max_ps(d,_mm_setzero_ps());
1480 d2 = _mm_mul_ps(d,d);
1481 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
1483 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
1485 /* Evaluate switch function */
1486 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1487 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv21,_mm_mul_ps(velec,dsw)) );
1488 velec = _mm_mul_ps(velec,sw);
1489 cutoff_mask = _mm_cmplt_ps(rsq21,rcutoff2);
1491 /* Update potential sum for this i atom from the interaction with this j atom. */
1492 velec = _mm_and_ps(velec,cutoff_mask);
1493 velec = _mm_andnot_ps(dummy_mask,velec);
1494 velecsum = _mm_add_ps(velecsum,velec);
1498 fscal = _mm_and_ps(fscal,cutoff_mask);
1500 fscal = _mm_andnot_ps(dummy_mask,fscal);
1502 /* Calculate temporary vectorial force */
1503 tx = _mm_mul_ps(fscal,dx21);
1504 ty = _mm_mul_ps(fscal,dy21);
1505 tz = _mm_mul_ps(fscal,dz21);
1507 /* Update vectorial force */
1508 fix2 = _mm_add_ps(fix2,tx);
1509 fiy2 = _mm_add_ps(fiy2,ty);
1510 fiz2 = _mm_add_ps(fiz2,tz);
1512 fjx1 = _mm_add_ps(fjx1,tx);
1513 fjy1 = _mm_add_ps(fjy1,ty);
1514 fjz1 = _mm_add_ps(fjz1,tz);
1518 /**************************
1519 * CALCULATE INTERACTIONS *
1520 **************************/
1522 if (gmx_mm_any_lt(rsq22,rcutoff2))
1525 r22 = _mm_mul_ps(rsq22,rinv22);
1526 r22 = _mm_andnot_ps(dummy_mask,r22);
1528 /* EWALD ELECTROSTATICS */
1530 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1531 ewrt = _mm_mul_ps(r22,ewtabscale);
1532 ewitab = _mm_cvttps_epi32(ewrt);
1533 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
1534 ewitab = _mm_slli_epi32(ewitab,2);
1535 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1536 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
1537 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
1538 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
1539 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
1540 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
1541 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
1542 velec = _mm_mul_ps(qq22,_mm_sub_ps(rinv22,velec));
1543 felec = _mm_mul_ps(_mm_mul_ps(qq22,rinv22),_mm_sub_ps(rinvsq22,felec));
1545 d = _mm_sub_ps(r22,rswitch);
1546 d = _mm_max_ps(d,_mm_setzero_ps());
1547 d2 = _mm_mul_ps(d,d);
1548 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
1550 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
1552 /* Evaluate switch function */
1553 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1554 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv22,_mm_mul_ps(velec,dsw)) );
1555 velec = _mm_mul_ps(velec,sw);
1556 cutoff_mask = _mm_cmplt_ps(rsq22,rcutoff2);
1558 /* Update potential sum for this i atom from the interaction with this j atom. */
1559 velec = _mm_and_ps(velec,cutoff_mask);
1560 velec = _mm_andnot_ps(dummy_mask,velec);
1561 velecsum = _mm_add_ps(velecsum,velec);
1565 fscal = _mm_and_ps(fscal,cutoff_mask);
1567 fscal = _mm_andnot_ps(dummy_mask,fscal);
1569 /* Calculate temporary vectorial force */
1570 tx = _mm_mul_ps(fscal,dx22);
1571 ty = _mm_mul_ps(fscal,dy22);
1572 tz = _mm_mul_ps(fscal,dz22);
1574 /* Update vectorial force */
1575 fix2 = _mm_add_ps(fix2,tx);
1576 fiy2 = _mm_add_ps(fiy2,ty);
1577 fiz2 = _mm_add_ps(fiz2,tz);
1579 fjx2 = _mm_add_ps(fjx2,tx);
1580 fjy2 = _mm_add_ps(fjy2,ty);
1581 fjz2 = _mm_add_ps(fjz2,tz);
1585 fjptrA = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
1586 fjptrB = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
1587 fjptrC = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
1588 fjptrD = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
1590 gmx_mm_decrement_3rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,
1591 fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,fjx2,fjy2,fjz2);
1593 /* Inner loop uses 594 flops */
1596 /* End of innermost loop */
1598 gmx_mm_update_iforce_3atom_swizzle_ps(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,
1599 f+i_coord_offset,fshift+i_shift_offset);
1602 /* Update potential energies */
1603 gmx_mm_update_1pot_ps(velecsum,kernel_data->energygrp_elec+ggid);
1605 /* Increment number of inner iterations */
1606 inneriter += j_index_end - j_index_start;
1608 /* Outer loop uses 19 flops */
1611 /* Increment number of outer iterations */
1614 /* Update outer/inner flops */
1616 inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_W3W3_VF,outeriter*19 + inneriter*594);
1619 * Gromacs nonbonded kernel: nb_kernel_ElecEwSw_VdwNone_GeomW3W3_F_sse4_1_single
1620 * Electrostatics interaction: Ewald
1621 * VdW interaction: None
1622 * Geometry: Water3-Water3
1623 * Calculate force/pot: Force
1626 nb_kernel_ElecEwSw_VdwNone_GeomW3W3_F_sse4_1_single
1627 (t_nblist * gmx_restrict nlist,
1628 rvec * gmx_restrict xx,
1629 rvec * gmx_restrict ff,
1630 t_forcerec * gmx_restrict fr,
1631 t_mdatoms * gmx_restrict mdatoms,
1632 nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
1633 t_nrnb * gmx_restrict nrnb)
1635 /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
1636 * just 0 for non-waters.
1637 * Suffixes A,B,C,D refer to j loop unrolling done with SSE, e.g. for the four different
1638 * jnr indices corresponding to data put in the four positions in the SIMD register.
1640 int i_shift_offset,i_coord_offset,outeriter,inneriter;
1641 int j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
1642 int jnrA,jnrB,jnrC,jnrD;
1643 int jnrlistA,jnrlistB,jnrlistC,jnrlistD;
1644 int j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
1645 int *iinr,*jindex,*jjnr,*shiftidx,*gid;
1646 real rcutoff_scalar;
1647 real *shiftvec,*fshift,*x,*f;
1648 real *fjptrA,*fjptrB,*fjptrC,*fjptrD;
1649 real scratch[4*DIM];
1650 __m128 tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
1652 __m128 ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
1654 __m128 ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
1656 __m128 ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
1657 int vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
1658 __m128 jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
1659 int vdwjidx1A,vdwjidx1B,vdwjidx1C,vdwjidx1D;
1660 __m128 jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
1661 int vdwjidx2A,vdwjidx2B,vdwjidx2C,vdwjidx2D;
1662 __m128 jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
1663 __m128 dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
1664 __m128 dx01,dy01,dz01,rsq01,rinv01,rinvsq01,r01,qq01,c6_01,c12_01;
1665 __m128 dx02,dy02,dz02,rsq02,rinv02,rinvsq02,r02,qq02,c6_02,c12_02;
1666 __m128 dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
1667 __m128 dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11;
1668 __m128 dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12;
1669 __m128 dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
1670 __m128 dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21;
1671 __m128 dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22;
1672 __m128 velec,felec,velecsum,facel,crf,krf,krf2;
1675 __m128 ewtabscale,eweps,sh_ewald,ewrt,ewtabhalfspace,ewtabF,ewtabFn,ewtabD,ewtabV;
1677 __m128 rswitch,swV3,swV4,swV5,swF2,swF3,swF4,d,d2,sw,dsw;
1678 real rswitch_scalar,d_scalar;
1679 __m128 dummy_mask,cutoff_mask;
1680 __m128 signbit = _mm_castsi128_ps( _mm_set1_epi32(0x80000000) );
1681 __m128 one = _mm_set1_ps(1.0);
1682 __m128 two = _mm_set1_ps(2.0);
1688 jindex = nlist->jindex;
1690 shiftidx = nlist->shift;
1692 shiftvec = fr->shift_vec[0];
1693 fshift = fr->fshift[0];
1694 facel = _mm_set1_ps(fr->epsfac);
1695 charge = mdatoms->chargeA;
1697 sh_ewald = _mm_set1_ps(fr->ic->sh_ewald);
1698 ewtab = fr->ic->tabq_coul_FDV0;
1699 ewtabscale = _mm_set1_ps(fr->ic->tabq_scale);
1700 ewtabhalfspace = _mm_set1_ps(0.5/fr->ic->tabq_scale);
1702 /* Setup water-specific parameters */
1703 inr = nlist->iinr[0];
1704 iq0 = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+0]));
1705 iq1 = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+1]));
1706 iq2 = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+2]));
1708 jq0 = _mm_set1_ps(charge[inr+0]);
1709 jq1 = _mm_set1_ps(charge[inr+1]);
1710 jq2 = _mm_set1_ps(charge[inr+2]);
1711 qq00 = _mm_mul_ps(iq0,jq0);
1712 qq01 = _mm_mul_ps(iq0,jq1);
1713 qq02 = _mm_mul_ps(iq0,jq2);
1714 qq10 = _mm_mul_ps(iq1,jq0);
1715 qq11 = _mm_mul_ps(iq1,jq1);
1716 qq12 = _mm_mul_ps(iq1,jq2);
1717 qq20 = _mm_mul_ps(iq2,jq0);
1718 qq21 = _mm_mul_ps(iq2,jq1);
1719 qq22 = _mm_mul_ps(iq2,jq2);
1721 /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
1722 rcutoff_scalar = fr->rcoulomb;
1723 rcutoff = _mm_set1_ps(rcutoff_scalar);
1724 rcutoff2 = _mm_mul_ps(rcutoff,rcutoff);
1726 rswitch_scalar = fr->rcoulomb_switch;
1727 rswitch = _mm_set1_ps(rswitch_scalar);
1728 /* Setup switch parameters */
1729 d_scalar = rcutoff_scalar-rswitch_scalar;
1730 d = _mm_set1_ps(d_scalar);
1731 swV3 = _mm_set1_ps(-10.0/(d_scalar*d_scalar*d_scalar));
1732 swV4 = _mm_set1_ps( 15.0/(d_scalar*d_scalar*d_scalar*d_scalar));
1733 swV5 = _mm_set1_ps( -6.0/(d_scalar*d_scalar*d_scalar*d_scalar*d_scalar));
1734 swF2 = _mm_set1_ps(-30.0/(d_scalar*d_scalar*d_scalar));
1735 swF3 = _mm_set1_ps( 60.0/(d_scalar*d_scalar*d_scalar*d_scalar));
1736 swF4 = _mm_set1_ps(-30.0/(d_scalar*d_scalar*d_scalar*d_scalar*d_scalar));
1738 /* Avoid stupid compiler warnings */
1739 jnrA = jnrB = jnrC = jnrD = 0;
1740 j_coord_offsetA = 0;
1741 j_coord_offsetB = 0;
1742 j_coord_offsetC = 0;
1743 j_coord_offsetD = 0;
1748 for(iidx=0;iidx<4*DIM;iidx++)
1750 scratch[iidx] = 0.0;
1753 /* Start outer loop over neighborlists */
1754 for(iidx=0; iidx<nri; iidx++)
1756 /* Load shift vector for this list */
1757 i_shift_offset = DIM*shiftidx[iidx];
1759 /* Load limits for loop over neighbors */
1760 j_index_start = jindex[iidx];
1761 j_index_end = jindex[iidx+1];
1763 /* Get outer coordinate index */
1765 i_coord_offset = DIM*inr;
1767 /* Load i particle coords and add shift vector */
1768 gmx_mm_load_shift_and_3rvec_broadcast_ps(shiftvec+i_shift_offset,x+i_coord_offset,
1769 &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2);
1771 fix0 = _mm_setzero_ps();
1772 fiy0 = _mm_setzero_ps();
1773 fiz0 = _mm_setzero_ps();
1774 fix1 = _mm_setzero_ps();
1775 fiy1 = _mm_setzero_ps();
1776 fiz1 = _mm_setzero_ps();
1777 fix2 = _mm_setzero_ps();
1778 fiy2 = _mm_setzero_ps();
1779 fiz2 = _mm_setzero_ps();
1781 /* Start inner kernel loop */
1782 for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
1785 /* Get j neighbor index, and coordinate index */
1787 jnrB = jjnr[jidx+1];
1788 jnrC = jjnr[jidx+2];
1789 jnrD = jjnr[jidx+3];
1790 j_coord_offsetA = DIM*jnrA;
1791 j_coord_offsetB = DIM*jnrB;
1792 j_coord_offsetC = DIM*jnrC;
1793 j_coord_offsetD = DIM*jnrD;
1795 /* load j atom coordinates */
1796 gmx_mm_load_3rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
1797 x+j_coord_offsetC,x+j_coord_offsetD,
1798 &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,&jy2,&jz2);
1800 /* Calculate displacement vector */
1801 dx00 = _mm_sub_ps(ix0,jx0);
1802 dy00 = _mm_sub_ps(iy0,jy0);
1803 dz00 = _mm_sub_ps(iz0,jz0);
1804 dx01 = _mm_sub_ps(ix0,jx1);
1805 dy01 = _mm_sub_ps(iy0,jy1);
1806 dz01 = _mm_sub_ps(iz0,jz1);
1807 dx02 = _mm_sub_ps(ix0,jx2);
1808 dy02 = _mm_sub_ps(iy0,jy2);
1809 dz02 = _mm_sub_ps(iz0,jz2);
1810 dx10 = _mm_sub_ps(ix1,jx0);
1811 dy10 = _mm_sub_ps(iy1,jy0);
1812 dz10 = _mm_sub_ps(iz1,jz0);
1813 dx11 = _mm_sub_ps(ix1,jx1);
1814 dy11 = _mm_sub_ps(iy1,jy1);
1815 dz11 = _mm_sub_ps(iz1,jz1);
1816 dx12 = _mm_sub_ps(ix1,jx2);
1817 dy12 = _mm_sub_ps(iy1,jy2);
1818 dz12 = _mm_sub_ps(iz1,jz2);
1819 dx20 = _mm_sub_ps(ix2,jx0);
1820 dy20 = _mm_sub_ps(iy2,jy0);
1821 dz20 = _mm_sub_ps(iz2,jz0);
1822 dx21 = _mm_sub_ps(ix2,jx1);
1823 dy21 = _mm_sub_ps(iy2,jy1);
1824 dz21 = _mm_sub_ps(iz2,jz1);
1825 dx22 = _mm_sub_ps(ix2,jx2);
1826 dy22 = _mm_sub_ps(iy2,jy2);
1827 dz22 = _mm_sub_ps(iz2,jz2);
1829 /* Calculate squared distance and things based on it */
1830 rsq00 = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
1831 rsq01 = gmx_mm_calc_rsq_ps(dx01,dy01,dz01);
1832 rsq02 = gmx_mm_calc_rsq_ps(dx02,dy02,dz02);
1833 rsq10 = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
1834 rsq11 = gmx_mm_calc_rsq_ps(dx11,dy11,dz11);
1835 rsq12 = gmx_mm_calc_rsq_ps(dx12,dy12,dz12);
1836 rsq20 = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
1837 rsq21 = gmx_mm_calc_rsq_ps(dx21,dy21,dz21);
1838 rsq22 = gmx_mm_calc_rsq_ps(dx22,dy22,dz22);
1840 rinv00 = gmx_mm_invsqrt_ps(rsq00);
1841 rinv01 = gmx_mm_invsqrt_ps(rsq01);
1842 rinv02 = gmx_mm_invsqrt_ps(rsq02);
1843 rinv10 = gmx_mm_invsqrt_ps(rsq10);
1844 rinv11 = gmx_mm_invsqrt_ps(rsq11);
1845 rinv12 = gmx_mm_invsqrt_ps(rsq12);
1846 rinv20 = gmx_mm_invsqrt_ps(rsq20);
1847 rinv21 = gmx_mm_invsqrt_ps(rsq21);
1848 rinv22 = gmx_mm_invsqrt_ps(rsq22);
1850 rinvsq00 = _mm_mul_ps(rinv00,rinv00);
1851 rinvsq01 = _mm_mul_ps(rinv01,rinv01);
1852 rinvsq02 = _mm_mul_ps(rinv02,rinv02);
1853 rinvsq10 = _mm_mul_ps(rinv10,rinv10);
1854 rinvsq11 = _mm_mul_ps(rinv11,rinv11);
1855 rinvsq12 = _mm_mul_ps(rinv12,rinv12);
1856 rinvsq20 = _mm_mul_ps(rinv20,rinv20);
1857 rinvsq21 = _mm_mul_ps(rinv21,rinv21);
1858 rinvsq22 = _mm_mul_ps(rinv22,rinv22);
1860 fjx0 = _mm_setzero_ps();
1861 fjy0 = _mm_setzero_ps();
1862 fjz0 = _mm_setzero_ps();
1863 fjx1 = _mm_setzero_ps();
1864 fjy1 = _mm_setzero_ps();
1865 fjz1 = _mm_setzero_ps();
1866 fjx2 = _mm_setzero_ps();
1867 fjy2 = _mm_setzero_ps();
1868 fjz2 = _mm_setzero_ps();
1870 /**************************
1871 * CALCULATE INTERACTIONS *
1872 **************************/
1874 if (gmx_mm_any_lt(rsq00,rcutoff2))
1877 r00 = _mm_mul_ps(rsq00,rinv00);
1879 /* EWALD ELECTROSTATICS */
1881 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1882 ewrt = _mm_mul_ps(r00,ewtabscale);
1883 ewitab = _mm_cvttps_epi32(ewrt);
1884 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
1885 ewitab = _mm_slli_epi32(ewitab,2);
1886 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1887 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
1888 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
1889 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
1890 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
1891 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
1892 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
1893 velec = _mm_mul_ps(qq00,_mm_sub_ps(rinv00,velec));
1894 felec = _mm_mul_ps(_mm_mul_ps(qq00,rinv00),_mm_sub_ps(rinvsq00,felec));
1896 d = _mm_sub_ps(r00,rswitch);
1897 d = _mm_max_ps(d,_mm_setzero_ps());
1898 d2 = _mm_mul_ps(d,d);
1899 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
1901 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
1903 /* Evaluate switch function */
1904 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1905 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv00,_mm_mul_ps(velec,dsw)) );
1906 cutoff_mask = _mm_cmplt_ps(rsq00,rcutoff2);
1910 fscal = _mm_and_ps(fscal,cutoff_mask);
1912 /* Calculate temporary vectorial force */
1913 tx = _mm_mul_ps(fscal,dx00);
1914 ty = _mm_mul_ps(fscal,dy00);
1915 tz = _mm_mul_ps(fscal,dz00);
1917 /* Update vectorial force */
1918 fix0 = _mm_add_ps(fix0,tx);
1919 fiy0 = _mm_add_ps(fiy0,ty);
1920 fiz0 = _mm_add_ps(fiz0,tz);
1922 fjx0 = _mm_add_ps(fjx0,tx);
1923 fjy0 = _mm_add_ps(fjy0,ty);
1924 fjz0 = _mm_add_ps(fjz0,tz);
1928 /**************************
1929 * CALCULATE INTERACTIONS *
1930 **************************/
1932 if (gmx_mm_any_lt(rsq01,rcutoff2))
1935 r01 = _mm_mul_ps(rsq01,rinv01);
1937 /* EWALD ELECTROSTATICS */
1939 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1940 ewrt = _mm_mul_ps(r01,ewtabscale);
1941 ewitab = _mm_cvttps_epi32(ewrt);
1942 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
1943 ewitab = _mm_slli_epi32(ewitab,2);
1944 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
1945 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
1946 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
1947 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
1948 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
1949 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
1950 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
1951 velec = _mm_mul_ps(qq01,_mm_sub_ps(rinv01,velec));
1952 felec = _mm_mul_ps(_mm_mul_ps(qq01,rinv01),_mm_sub_ps(rinvsq01,felec));
1954 d = _mm_sub_ps(r01,rswitch);
1955 d = _mm_max_ps(d,_mm_setzero_ps());
1956 d2 = _mm_mul_ps(d,d);
1957 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
1959 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
1961 /* Evaluate switch function */
1962 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1963 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv01,_mm_mul_ps(velec,dsw)) );
1964 cutoff_mask = _mm_cmplt_ps(rsq01,rcutoff2);
1968 fscal = _mm_and_ps(fscal,cutoff_mask);
1970 /* Calculate temporary vectorial force */
1971 tx = _mm_mul_ps(fscal,dx01);
1972 ty = _mm_mul_ps(fscal,dy01);
1973 tz = _mm_mul_ps(fscal,dz01);
1975 /* Update vectorial force */
1976 fix0 = _mm_add_ps(fix0,tx);
1977 fiy0 = _mm_add_ps(fiy0,ty);
1978 fiz0 = _mm_add_ps(fiz0,tz);
1980 fjx1 = _mm_add_ps(fjx1,tx);
1981 fjy1 = _mm_add_ps(fjy1,ty);
1982 fjz1 = _mm_add_ps(fjz1,tz);
1986 /**************************
1987 * CALCULATE INTERACTIONS *
1988 **************************/
1990 if (gmx_mm_any_lt(rsq02,rcutoff2))
1993 r02 = _mm_mul_ps(rsq02,rinv02);
1995 /* EWALD ELECTROSTATICS */
1997 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1998 ewrt = _mm_mul_ps(r02,ewtabscale);
1999 ewitab = _mm_cvttps_epi32(ewrt);
2000 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
2001 ewitab = _mm_slli_epi32(ewitab,2);
2002 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2003 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2004 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
2005 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
2006 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
2007 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
2008 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
2009 velec = _mm_mul_ps(qq02,_mm_sub_ps(rinv02,velec));
2010 felec = _mm_mul_ps(_mm_mul_ps(qq02,rinv02),_mm_sub_ps(rinvsq02,felec));
2012 d = _mm_sub_ps(r02,rswitch);
2013 d = _mm_max_ps(d,_mm_setzero_ps());
2014 d2 = _mm_mul_ps(d,d);
2015 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
2017 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
2019 /* Evaluate switch function */
2020 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2021 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv02,_mm_mul_ps(velec,dsw)) );
2022 cutoff_mask = _mm_cmplt_ps(rsq02,rcutoff2);
2026 fscal = _mm_and_ps(fscal,cutoff_mask);
2028 /* Calculate temporary vectorial force */
2029 tx = _mm_mul_ps(fscal,dx02);
2030 ty = _mm_mul_ps(fscal,dy02);
2031 tz = _mm_mul_ps(fscal,dz02);
2033 /* Update vectorial force */
2034 fix0 = _mm_add_ps(fix0,tx);
2035 fiy0 = _mm_add_ps(fiy0,ty);
2036 fiz0 = _mm_add_ps(fiz0,tz);
2038 fjx2 = _mm_add_ps(fjx2,tx);
2039 fjy2 = _mm_add_ps(fjy2,ty);
2040 fjz2 = _mm_add_ps(fjz2,tz);
2044 /**************************
2045 * CALCULATE INTERACTIONS *
2046 **************************/
2048 if (gmx_mm_any_lt(rsq10,rcutoff2))
2051 r10 = _mm_mul_ps(rsq10,rinv10);
2053 /* EWALD ELECTROSTATICS */
2055 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2056 ewrt = _mm_mul_ps(r10,ewtabscale);
2057 ewitab = _mm_cvttps_epi32(ewrt);
2058 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
2059 ewitab = _mm_slli_epi32(ewitab,2);
2060 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2061 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2062 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
2063 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
2064 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
2065 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
2066 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
2067 velec = _mm_mul_ps(qq10,_mm_sub_ps(rinv10,velec));
2068 felec = _mm_mul_ps(_mm_mul_ps(qq10,rinv10),_mm_sub_ps(rinvsq10,felec));
2070 d = _mm_sub_ps(r10,rswitch);
2071 d = _mm_max_ps(d,_mm_setzero_ps());
2072 d2 = _mm_mul_ps(d,d);
2073 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
2075 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
2077 /* Evaluate switch function */
2078 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2079 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv10,_mm_mul_ps(velec,dsw)) );
2080 cutoff_mask = _mm_cmplt_ps(rsq10,rcutoff2);
2084 fscal = _mm_and_ps(fscal,cutoff_mask);
2086 /* Calculate temporary vectorial force */
2087 tx = _mm_mul_ps(fscal,dx10);
2088 ty = _mm_mul_ps(fscal,dy10);
2089 tz = _mm_mul_ps(fscal,dz10);
2091 /* Update vectorial force */
2092 fix1 = _mm_add_ps(fix1,tx);
2093 fiy1 = _mm_add_ps(fiy1,ty);
2094 fiz1 = _mm_add_ps(fiz1,tz);
2096 fjx0 = _mm_add_ps(fjx0,tx);
2097 fjy0 = _mm_add_ps(fjy0,ty);
2098 fjz0 = _mm_add_ps(fjz0,tz);
2102 /**************************
2103 * CALCULATE INTERACTIONS *
2104 **************************/
2106 if (gmx_mm_any_lt(rsq11,rcutoff2))
2109 r11 = _mm_mul_ps(rsq11,rinv11);
2111 /* EWALD ELECTROSTATICS */
2113 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2114 ewrt = _mm_mul_ps(r11,ewtabscale);
2115 ewitab = _mm_cvttps_epi32(ewrt);
2116 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
2117 ewitab = _mm_slli_epi32(ewitab,2);
2118 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2119 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2120 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
2121 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
2122 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
2123 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
2124 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
2125 velec = _mm_mul_ps(qq11,_mm_sub_ps(rinv11,velec));
2126 felec = _mm_mul_ps(_mm_mul_ps(qq11,rinv11),_mm_sub_ps(rinvsq11,felec));
2128 d = _mm_sub_ps(r11,rswitch);
2129 d = _mm_max_ps(d,_mm_setzero_ps());
2130 d2 = _mm_mul_ps(d,d);
2131 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
2133 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
2135 /* Evaluate switch function */
2136 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2137 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv11,_mm_mul_ps(velec,dsw)) );
2138 cutoff_mask = _mm_cmplt_ps(rsq11,rcutoff2);
2142 fscal = _mm_and_ps(fscal,cutoff_mask);
2144 /* Calculate temporary vectorial force */
2145 tx = _mm_mul_ps(fscal,dx11);
2146 ty = _mm_mul_ps(fscal,dy11);
2147 tz = _mm_mul_ps(fscal,dz11);
2149 /* Update vectorial force */
2150 fix1 = _mm_add_ps(fix1,tx);
2151 fiy1 = _mm_add_ps(fiy1,ty);
2152 fiz1 = _mm_add_ps(fiz1,tz);
2154 fjx1 = _mm_add_ps(fjx1,tx);
2155 fjy1 = _mm_add_ps(fjy1,ty);
2156 fjz1 = _mm_add_ps(fjz1,tz);
2160 /**************************
2161 * CALCULATE INTERACTIONS *
2162 **************************/
2164 if (gmx_mm_any_lt(rsq12,rcutoff2))
2167 r12 = _mm_mul_ps(rsq12,rinv12);
2169 /* EWALD ELECTROSTATICS */
2171 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2172 ewrt = _mm_mul_ps(r12,ewtabscale);
2173 ewitab = _mm_cvttps_epi32(ewrt);
2174 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
2175 ewitab = _mm_slli_epi32(ewitab,2);
2176 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2177 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2178 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
2179 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
2180 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
2181 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
2182 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
2183 velec = _mm_mul_ps(qq12,_mm_sub_ps(rinv12,velec));
2184 felec = _mm_mul_ps(_mm_mul_ps(qq12,rinv12),_mm_sub_ps(rinvsq12,felec));
2186 d = _mm_sub_ps(r12,rswitch);
2187 d = _mm_max_ps(d,_mm_setzero_ps());
2188 d2 = _mm_mul_ps(d,d);
2189 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
2191 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
2193 /* Evaluate switch function */
2194 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2195 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv12,_mm_mul_ps(velec,dsw)) );
2196 cutoff_mask = _mm_cmplt_ps(rsq12,rcutoff2);
2200 fscal = _mm_and_ps(fscal,cutoff_mask);
2202 /* Calculate temporary vectorial force */
2203 tx = _mm_mul_ps(fscal,dx12);
2204 ty = _mm_mul_ps(fscal,dy12);
2205 tz = _mm_mul_ps(fscal,dz12);
2207 /* Update vectorial force */
2208 fix1 = _mm_add_ps(fix1,tx);
2209 fiy1 = _mm_add_ps(fiy1,ty);
2210 fiz1 = _mm_add_ps(fiz1,tz);
2212 fjx2 = _mm_add_ps(fjx2,tx);
2213 fjy2 = _mm_add_ps(fjy2,ty);
2214 fjz2 = _mm_add_ps(fjz2,tz);
2218 /**************************
2219 * CALCULATE INTERACTIONS *
2220 **************************/
2222 if (gmx_mm_any_lt(rsq20,rcutoff2))
2225 r20 = _mm_mul_ps(rsq20,rinv20);
2227 /* EWALD ELECTROSTATICS */
2229 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2230 ewrt = _mm_mul_ps(r20,ewtabscale);
2231 ewitab = _mm_cvttps_epi32(ewrt);
2232 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
2233 ewitab = _mm_slli_epi32(ewitab,2);
2234 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2235 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2236 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
2237 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
2238 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
2239 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
2240 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
2241 velec = _mm_mul_ps(qq20,_mm_sub_ps(rinv20,velec));
2242 felec = _mm_mul_ps(_mm_mul_ps(qq20,rinv20),_mm_sub_ps(rinvsq20,felec));
2244 d = _mm_sub_ps(r20,rswitch);
2245 d = _mm_max_ps(d,_mm_setzero_ps());
2246 d2 = _mm_mul_ps(d,d);
2247 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
2249 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
2251 /* Evaluate switch function */
2252 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2253 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv20,_mm_mul_ps(velec,dsw)) );
2254 cutoff_mask = _mm_cmplt_ps(rsq20,rcutoff2);
2258 fscal = _mm_and_ps(fscal,cutoff_mask);
2260 /* Calculate temporary vectorial force */
2261 tx = _mm_mul_ps(fscal,dx20);
2262 ty = _mm_mul_ps(fscal,dy20);
2263 tz = _mm_mul_ps(fscal,dz20);
2265 /* Update vectorial force */
2266 fix2 = _mm_add_ps(fix2,tx);
2267 fiy2 = _mm_add_ps(fiy2,ty);
2268 fiz2 = _mm_add_ps(fiz2,tz);
2270 fjx0 = _mm_add_ps(fjx0,tx);
2271 fjy0 = _mm_add_ps(fjy0,ty);
2272 fjz0 = _mm_add_ps(fjz0,tz);
2276 /**************************
2277 * CALCULATE INTERACTIONS *
2278 **************************/
2280 if (gmx_mm_any_lt(rsq21,rcutoff2))
2283 r21 = _mm_mul_ps(rsq21,rinv21);
2285 /* EWALD ELECTROSTATICS */
2287 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2288 ewrt = _mm_mul_ps(r21,ewtabscale);
2289 ewitab = _mm_cvttps_epi32(ewrt);
2290 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
2291 ewitab = _mm_slli_epi32(ewitab,2);
2292 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2293 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2294 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
2295 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
2296 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
2297 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
2298 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
2299 velec = _mm_mul_ps(qq21,_mm_sub_ps(rinv21,velec));
2300 felec = _mm_mul_ps(_mm_mul_ps(qq21,rinv21),_mm_sub_ps(rinvsq21,felec));
2302 d = _mm_sub_ps(r21,rswitch);
2303 d = _mm_max_ps(d,_mm_setzero_ps());
2304 d2 = _mm_mul_ps(d,d);
2305 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
2307 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
2309 /* Evaluate switch function */
2310 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2311 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv21,_mm_mul_ps(velec,dsw)) );
2312 cutoff_mask = _mm_cmplt_ps(rsq21,rcutoff2);
2316 fscal = _mm_and_ps(fscal,cutoff_mask);
2318 /* Calculate temporary vectorial force */
2319 tx = _mm_mul_ps(fscal,dx21);
2320 ty = _mm_mul_ps(fscal,dy21);
2321 tz = _mm_mul_ps(fscal,dz21);
2323 /* Update vectorial force */
2324 fix2 = _mm_add_ps(fix2,tx);
2325 fiy2 = _mm_add_ps(fiy2,ty);
2326 fiz2 = _mm_add_ps(fiz2,tz);
2328 fjx1 = _mm_add_ps(fjx1,tx);
2329 fjy1 = _mm_add_ps(fjy1,ty);
2330 fjz1 = _mm_add_ps(fjz1,tz);
2334 /**************************
2335 * CALCULATE INTERACTIONS *
2336 **************************/
2338 if (gmx_mm_any_lt(rsq22,rcutoff2))
2341 r22 = _mm_mul_ps(rsq22,rinv22);
2343 /* EWALD ELECTROSTATICS */
2345 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2346 ewrt = _mm_mul_ps(r22,ewtabscale);
2347 ewitab = _mm_cvttps_epi32(ewrt);
2348 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
2349 ewitab = _mm_slli_epi32(ewitab,2);
2350 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2351 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2352 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
2353 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
2354 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
2355 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
2356 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
2357 velec = _mm_mul_ps(qq22,_mm_sub_ps(rinv22,velec));
2358 felec = _mm_mul_ps(_mm_mul_ps(qq22,rinv22),_mm_sub_ps(rinvsq22,felec));
2360 d = _mm_sub_ps(r22,rswitch);
2361 d = _mm_max_ps(d,_mm_setzero_ps());
2362 d2 = _mm_mul_ps(d,d);
2363 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
2365 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
2367 /* Evaluate switch function */
2368 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2369 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv22,_mm_mul_ps(velec,dsw)) );
2370 cutoff_mask = _mm_cmplt_ps(rsq22,rcutoff2);
2374 fscal = _mm_and_ps(fscal,cutoff_mask);
2376 /* Calculate temporary vectorial force */
2377 tx = _mm_mul_ps(fscal,dx22);
2378 ty = _mm_mul_ps(fscal,dy22);
2379 tz = _mm_mul_ps(fscal,dz22);
2381 /* Update vectorial force */
2382 fix2 = _mm_add_ps(fix2,tx);
2383 fiy2 = _mm_add_ps(fiy2,ty);
2384 fiz2 = _mm_add_ps(fiz2,tz);
2386 fjx2 = _mm_add_ps(fjx2,tx);
2387 fjy2 = _mm_add_ps(fjy2,ty);
2388 fjz2 = _mm_add_ps(fjz2,tz);
2392 fjptrA = f+j_coord_offsetA;
2393 fjptrB = f+j_coord_offsetB;
2394 fjptrC = f+j_coord_offsetC;
2395 fjptrD = f+j_coord_offsetD;
2397 gmx_mm_decrement_3rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,
2398 fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,fjx2,fjy2,fjz2);
2400 /* Inner loop uses 558 flops */
2403 if(jidx<j_index_end)
2406 /* Get j neighbor index, and coordinate index */
2407 jnrlistA = jjnr[jidx];
2408 jnrlistB = jjnr[jidx+1];
2409 jnrlistC = jjnr[jidx+2];
2410 jnrlistD = jjnr[jidx+3];
2411 /* Sign of each element will be negative for non-real atoms.
2412 * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
2413 * so use it as val = _mm_andnot_ps(mask,val) to clear dummy entries.
2415 dummy_mask = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
2416 jnrA = (jnrlistA>=0) ? jnrlistA : 0;
2417 jnrB = (jnrlistB>=0) ? jnrlistB : 0;
2418 jnrC = (jnrlistC>=0) ? jnrlistC : 0;
2419 jnrD = (jnrlistD>=0) ? jnrlistD : 0;
2420 j_coord_offsetA = DIM*jnrA;
2421 j_coord_offsetB = DIM*jnrB;
2422 j_coord_offsetC = DIM*jnrC;
2423 j_coord_offsetD = DIM*jnrD;
2425 /* load j atom coordinates */
2426 gmx_mm_load_3rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
2427 x+j_coord_offsetC,x+j_coord_offsetD,
2428 &jx0,&jy0,&jz0,&jx1,&jy1,&jz1,&jx2,&jy2,&jz2);
2430 /* Calculate displacement vector */
2431 dx00 = _mm_sub_ps(ix0,jx0);
2432 dy00 = _mm_sub_ps(iy0,jy0);
2433 dz00 = _mm_sub_ps(iz0,jz0);
2434 dx01 = _mm_sub_ps(ix0,jx1);
2435 dy01 = _mm_sub_ps(iy0,jy1);
2436 dz01 = _mm_sub_ps(iz0,jz1);
2437 dx02 = _mm_sub_ps(ix0,jx2);
2438 dy02 = _mm_sub_ps(iy0,jy2);
2439 dz02 = _mm_sub_ps(iz0,jz2);
2440 dx10 = _mm_sub_ps(ix1,jx0);
2441 dy10 = _mm_sub_ps(iy1,jy0);
2442 dz10 = _mm_sub_ps(iz1,jz0);
2443 dx11 = _mm_sub_ps(ix1,jx1);
2444 dy11 = _mm_sub_ps(iy1,jy1);
2445 dz11 = _mm_sub_ps(iz1,jz1);
2446 dx12 = _mm_sub_ps(ix1,jx2);
2447 dy12 = _mm_sub_ps(iy1,jy2);
2448 dz12 = _mm_sub_ps(iz1,jz2);
2449 dx20 = _mm_sub_ps(ix2,jx0);
2450 dy20 = _mm_sub_ps(iy2,jy0);
2451 dz20 = _mm_sub_ps(iz2,jz0);
2452 dx21 = _mm_sub_ps(ix2,jx1);
2453 dy21 = _mm_sub_ps(iy2,jy1);
2454 dz21 = _mm_sub_ps(iz2,jz1);
2455 dx22 = _mm_sub_ps(ix2,jx2);
2456 dy22 = _mm_sub_ps(iy2,jy2);
2457 dz22 = _mm_sub_ps(iz2,jz2);
2459 /* Calculate squared distance and things based on it */
2460 rsq00 = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
2461 rsq01 = gmx_mm_calc_rsq_ps(dx01,dy01,dz01);
2462 rsq02 = gmx_mm_calc_rsq_ps(dx02,dy02,dz02);
2463 rsq10 = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
2464 rsq11 = gmx_mm_calc_rsq_ps(dx11,dy11,dz11);
2465 rsq12 = gmx_mm_calc_rsq_ps(dx12,dy12,dz12);
2466 rsq20 = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
2467 rsq21 = gmx_mm_calc_rsq_ps(dx21,dy21,dz21);
2468 rsq22 = gmx_mm_calc_rsq_ps(dx22,dy22,dz22);
2470 rinv00 = gmx_mm_invsqrt_ps(rsq00);
2471 rinv01 = gmx_mm_invsqrt_ps(rsq01);
2472 rinv02 = gmx_mm_invsqrt_ps(rsq02);
2473 rinv10 = gmx_mm_invsqrt_ps(rsq10);
2474 rinv11 = gmx_mm_invsqrt_ps(rsq11);
2475 rinv12 = gmx_mm_invsqrt_ps(rsq12);
2476 rinv20 = gmx_mm_invsqrt_ps(rsq20);
2477 rinv21 = gmx_mm_invsqrt_ps(rsq21);
2478 rinv22 = gmx_mm_invsqrt_ps(rsq22);
2480 rinvsq00 = _mm_mul_ps(rinv00,rinv00);
2481 rinvsq01 = _mm_mul_ps(rinv01,rinv01);
2482 rinvsq02 = _mm_mul_ps(rinv02,rinv02);
2483 rinvsq10 = _mm_mul_ps(rinv10,rinv10);
2484 rinvsq11 = _mm_mul_ps(rinv11,rinv11);
2485 rinvsq12 = _mm_mul_ps(rinv12,rinv12);
2486 rinvsq20 = _mm_mul_ps(rinv20,rinv20);
2487 rinvsq21 = _mm_mul_ps(rinv21,rinv21);
2488 rinvsq22 = _mm_mul_ps(rinv22,rinv22);
2490 fjx0 = _mm_setzero_ps();
2491 fjy0 = _mm_setzero_ps();
2492 fjz0 = _mm_setzero_ps();
2493 fjx1 = _mm_setzero_ps();
2494 fjy1 = _mm_setzero_ps();
2495 fjz1 = _mm_setzero_ps();
2496 fjx2 = _mm_setzero_ps();
2497 fjy2 = _mm_setzero_ps();
2498 fjz2 = _mm_setzero_ps();
2500 /**************************
2501 * CALCULATE INTERACTIONS *
2502 **************************/
2504 if (gmx_mm_any_lt(rsq00,rcutoff2))
2507 r00 = _mm_mul_ps(rsq00,rinv00);
2508 r00 = _mm_andnot_ps(dummy_mask,r00);
2510 /* EWALD ELECTROSTATICS */
2512 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2513 ewrt = _mm_mul_ps(r00,ewtabscale);
2514 ewitab = _mm_cvttps_epi32(ewrt);
2515 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
2516 ewitab = _mm_slli_epi32(ewitab,2);
2517 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2518 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2519 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
2520 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
2521 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
2522 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
2523 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
2524 velec = _mm_mul_ps(qq00,_mm_sub_ps(rinv00,velec));
2525 felec = _mm_mul_ps(_mm_mul_ps(qq00,rinv00),_mm_sub_ps(rinvsq00,felec));
2527 d = _mm_sub_ps(r00,rswitch);
2528 d = _mm_max_ps(d,_mm_setzero_ps());
2529 d2 = _mm_mul_ps(d,d);
2530 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
2532 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
2534 /* Evaluate switch function */
2535 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2536 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv00,_mm_mul_ps(velec,dsw)) );
2537 cutoff_mask = _mm_cmplt_ps(rsq00,rcutoff2);
2541 fscal = _mm_and_ps(fscal,cutoff_mask);
2543 fscal = _mm_andnot_ps(dummy_mask,fscal);
2545 /* Calculate temporary vectorial force */
2546 tx = _mm_mul_ps(fscal,dx00);
2547 ty = _mm_mul_ps(fscal,dy00);
2548 tz = _mm_mul_ps(fscal,dz00);
2550 /* Update vectorial force */
2551 fix0 = _mm_add_ps(fix0,tx);
2552 fiy0 = _mm_add_ps(fiy0,ty);
2553 fiz0 = _mm_add_ps(fiz0,tz);
2555 fjx0 = _mm_add_ps(fjx0,tx);
2556 fjy0 = _mm_add_ps(fjy0,ty);
2557 fjz0 = _mm_add_ps(fjz0,tz);
2561 /**************************
2562 * CALCULATE INTERACTIONS *
2563 **************************/
2565 if (gmx_mm_any_lt(rsq01,rcutoff2))
2568 r01 = _mm_mul_ps(rsq01,rinv01);
2569 r01 = _mm_andnot_ps(dummy_mask,r01);
2571 /* EWALD ELECTROSTATICS */
2573 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2574 ewrt = _mm_mul_ps(r01,ewtabscale);
2575 ewitab = _mm_cvttps_epi32(ewrt);
2576 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
2577 ewitab = _mm_slli_epi32(ewitab,2);
2578 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2579 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2580 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
2581 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
2582 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
2583 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
2584 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
2585 velec = _mm_mul_ps(qq01,_mm_sub_ps(rinv01,velec));
2586 felec = _mm_mul_ps(_mm_mul_ps(qq01,rinv01),_mm_sub_ps(rinvsq01,felec));
2588 d = _mm_sub_ps(r01,rswitch);
2589 d = _mm_max_ps(d,_mm_setzero_ps());
2590 d2 = _mm_mul_ps(d,d);
2591 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
2593 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
2595 /* Evaluate switch function */
2596 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2597 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv01,_mm_mul_ps(velec,dsw)) );
2598 cutoff_mask = _mm_cmplt_ps(rsq01,rcutoff2);
2602 fscal = _mm_and_ps(fscal,cutoff_mask);
2604 fscal = _mm_andnot_ps(dummy_mask,fscal);
2606 /* Calculate temporary vectorial force */
2607 tx = _mm_mul_ps(fscal,dx01);
2608 ty = _mm_mul_ps(fscal,dy01);
2609 tz = _mm_mul_ps(fscal,dz01);
2611 /* Update vectorial force */
2612 fix0 = _mm_add_ps(fix0,tx);
2613 fiy0 = _mm_add_ps(fiy0,ty);
2614 fiz0 = _mm_add_ps(fiz0,tz);
2616 fjx1 = _mm_add_ps(fjx1,tx);
2617 fjy1 = _mm_add_ps(fjy1,ty);
2618 fjz1 = _mm_add_ps(fjz1,tz);
2622 /**************************
2623 * CALCULATE INTERACTIONS *
2624 **************************/
2626 if (gmx_mm_any_lt(rsq02,rcutoff2))
2629 r02 = _mm_mul_ps(rsq02,rinv02);
2630 r02 = _mm_andnot_ps(dummy_mask,r02);
2632 /* EWALD ELECTROSTATICS */
2634 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2635 ewrt = _mm_mul_ps(r02,ewtabscale);
2636 ewitab = _mm_cvttps_epi32(ewrt);
2637 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
2638 ewitab = _mm_slli_epi32(ewitab,2);
2639 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2640 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2641 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
2642 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
2643 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
2644 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
2645 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
2646 velec = _mm_mul_ps(qq02,_mm_sub_ps(rinv02,velec));
2647 felec = _mm_mul_ps(_mm_mul_ps(qq02,rinv02),_mm_sub_ps(rinvsq02,felec));
2649 d = _mm_sub_ps(r02,rswitch);
2650 d = _mm_max_ps(d,_mm_setzero_ps());
2651 d2 = _mm_mul_ps(d,d);
2652 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
2654 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
2656 /* Evaluate switch function */
2657 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2658 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv02,_mm_mul_ps(velec,dsw)) );
2659 cutoff_mask = _mm_cmplt_ps(rsq02,rcutoff2);
2663 fscal = _mm_and_ps(fscal,cutoff_mask);
2665 fscal = _mm_andnot_ps(dummy_mask,fscal);
2667 /* Calculate temporary vectorial force */
2668 tx = _mm_mul_ps(fscal,dx02);
2669 ty = _mm_mul_ps(fscal,dy02);
2670 tz = _mm_mul_ps(fscal,dz02);
2672 /* Update vectorial force */
2673 fix0 = _mm_add_ps(fix0,tx);
2674 fiy0 = _mm_add_ps(fiy0,ty);
2675 fiz0 = _mm_add_ps(fiz0,tz);
2677 fjx2 = _mm_add_ps(fjx2,tx);
2678 fjy2 = _mm_add_ps(fjy2,ty);
2679 fjz2 = _mm_add_ps(fjz2,tz);
2683 /**************************
2684 * CALCULATE INTERACTIONS *
2685 **************************/
2687 if (gmx_mm_any_lt(rsq10,rcutoff2))
2690 r10 = _mm_mul_ps(rsq10,rinv10);
2691 r10 = _mm_andnot_ps(dummy_mask,r10);
2693 /* EWALD ELECTROSTATICS */
2695 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2696 ewrt = _mm_mul_ps(r10,ewtabscale);
2697 ewitab = _mm_cvttps_epi32(ewrt);
2698 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
2699 ewitab = _mm_slli_epi32(ewitab,2);
2700 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2701 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2702 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
2703 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
2704 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
2705 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
2706 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
2707 velec = _mm_mul_ps(qq10,_mm_sub_ps(rinv10,velec));
2708 felec = _mm_mul_ps(_mm_mul_ps(qq10,rinv10),_mm_sub_ps(rinvsq10,felec));
2710 d = _mm_sub_ps(r10,rswitch);
2711 d = _mm_max_ps(d,_mm_setzero_ps());
2712 d2 = _mm_mul_ps(d,d);
2713 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
2715 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
2717 /* Evaluate switch function */
2718 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2719 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv10,_mm_mul_ps(velec,dsw)) );
2720 cutoff_mask = _mm_cmplt_ps(rsq10,rcutoff2);
2724 fscal = _mm_and_ps(fscal,cutoff_mask);
2726 fscal = _mm_andnot_ps(dummy_mask,fscal);
2728 /* Calculate temporary vectorial force */
2729 tx = _mm_mul_ps(fscal,dx10);
2730 ty = _mm_mul_ps(fscal,dy10);
2731 tz = _mm_mul_ps(fscal,dz10);
2733 /* Update vectorial force */
2734 fix1 = _mm_add_ps(fix1,tx);
2735 fiy1 = _mm_add_ps(fiy1,ty);
2736 fiz1 = _mm_add_ps(fiz1,tz);
2738 fjx0 = _mm_add_ps(fjx0,tx);
2739 fjy0 = _mm_add_ps(fjy0,ty);
2740 fjz0 = _mm_add_ps(fjz0,tz);
2744 /**************************
2745 * CALCULATE INTERACTIONS *
2746 **************************/
2748 if (gmx_mm_any_lt(rsq11,rcutoff2))
2751 r11 = _mm_mul_ps(rsq11,rinv11);
2752 r11 = _mm_andnot_ps(dummy_mask,r11);
2754 /* EWALD ELECTROSTATICS */
2756 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2757 ewrt = _mm_mul_ps(r11,ewtabscale);
2758 ewitab = _mm_cvttps_epi32(ewrt);
2759 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
2760 ewitab = _mm_slli_epi32(ewitab,2);
2761 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2762 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2763 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
2764 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
2765 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
2766 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
2767 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
2768 velec = _mm_mul_ps(qq11,_mm_sub_ps(rinv11,velec));
2769 felec = _mm_mul_ps(_mm_mul_ps(qq11,rinv11),_mm_sub_ps(rinvsq11,felec));
2771 d = _mm_sub_ps(r11,rswitch);
2772 d = _mm_max_ps(d,_mm_setzero_ps());
2773 d2 = _mm_mul_ps(d,d);
2774 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
2776 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
2778 /* Evaluate switch function */
2779 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2780 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv11,_mm_mul_ps(velec,dsw)) );
2781 cutoff_mask = _mm_cmplt_ps(rsq11,rcutoff2);
2785 fscal = _mm_and_ps(fscal,cutoff_mask);
2787 fscal = _mm_andnot_ps(dummy_mask,fscal);
2789 /* Calculate temporary vectorial force */
2790 tx = _mm_mul_ps(fscal,dx11);
2791 ty = _mm_mul_ps(fscal,dy11);
2792 tz = _mm_mul_ps(fscal,dz11);
2794 /* Update vectorial force */
2795 fix1 = _mm_add_ps(fix1,tx);
2796 fiy1 = _mm_add_ps(fiy1,ty);
2797 fiz1 = _mm_add_ps(fiz1,tz);
2799 fjx1 = _mm_add_ps(fjx1,tx);
2800 fjy1 = _mm_add_ps(fjy1,ty);
2801 fjz1 = _mm_add_ps(fjz1,tz);
2805 /**************************
2806 * CALCULATE INTERACTIONS *
2807 **************************/
2809 if (gmx_mm_any_lt(rsq12,rcutoff2))
2812 r12 = _mm_mul_ps(rsq12,rinv12);
2813 r12 = _mm_andnot_ps(dummy_mask,r12);
2815 /* EWALD ELECTROSTATICS */
2817 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2818 ewrt = _mm_mul_ps(r12,ewtabscale);
2819 ewitab = _mm_cvttps_epi32(ewrt);
2820 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
2821 ewitab = _mm_slli_epi32(ewitab,2);
2822 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2823 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2824 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
2825 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
2826 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
2827 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
2828 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
2829 velec = _mm_mul_ps(qq12,_mm_sub_ps(rinv12,velec));
2830 felec = _mm_mul_ps(_mm_mul_ps(qq12,rinv12),_mm_sub_ps(rinvsq12,felec));
2832 d = _mm_sub_ps(r12,rswitch);
2833 d = _mm_max_ps(d,_mm_setzero_ps());
2834 d2 = _mm_mul_ps(d,d);
2835 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
2837 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
2839 /* Evaluate switch function */
2840 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2841 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv12,_mm_mul_ps(velec,dsw)) );
2842 cutoff_mask = _mm_cmplt_ps(rsq12,rcutoff2);
2846 fscal = _mm_and_ps(fscal,cutoff_mask);
2848 fscal = _mm_andnot_ps(dummy_mask,fscal);
2850 /* Calculate temporary vectorial force */
2851 tx = _mm_mul_ps(fscal,dx12);
2852 ty = _mm_mul_ps(fscal,dy12);
2853 tz = _mm_mul_ps(fscal,dz12);
2855 /* Update vectorial force */
2856 fix1 = _mm_add_ps(fix1,tx);
2857 fiy1 = _mm_add_ps(fiy1,ty);
2858 fiz1 = _mm_add_ps(fiz1,tz);
2860 fjx2 = _mm_add_ps(fjx2,tx);
2861 fjy2 = _mm_add_ps(fjy2,ty);
2862 fjz2 = _mm_add_ps(fjz2,tz);
2866 /**************************
2867 * CALCULATE INTERACTIONS *
2868 **************************/
2870 if (gmx_mm_any_lt(rsq20,rcutoff2))
2873 r20 = _mm_mul_ps(rsq20,rinv20);
2874 r20 = _mm_andnot_ps(dummy_mask,r20);
2876 /* EWALD ELECTROSTATICS */
2878 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2879 ewrt = _mm_mul_ps(r20,ewtabscale);
2880 ewitab = _mm_cvttps_epi32(ewrt);
2881 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
2882 ewitab = _mm_slli_epi32(ewitab,2);
2883 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2884 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2885 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
2886 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
2887 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
2888 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
2889 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
2890 velec = _mm_mul_ps(qq20,_mm_sub_ps(rinv20,velec));
2891 felec = _mm_mul_ps(_mm_mul_ps(qq20,rinv20),_mm_sub_ps(rinvsq20,felec));
2893 d = _mm_sub_ps(r20,rswitch);
2894 d = _mm_max_ps(d,_mm_setzero_ps());
2895 d2 = _mm_mul_ps(d,d);
2896 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
2898 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
2900 /* Evaluate switch function */
2901 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2902 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv20,_mm_mul_ps(velec,dsw)) );
2903 cutoff_mask = _mm_cmplt_ps(rsq20,rcutoff2);
2907 fscal = _mm_and_ps(fscal,cutoff_mask);
2909 fscal = _mm_andnot_ps(dummy_mask,fscal);
2911 /* Calculate temporary vectorial force */
2912 tx = _mm_mul_ps(fscal,dx20);
2913 ty = _mm_mul_ps(fscal,dy20);
2914 tz = _mm_mul_ps(fscal,dz20);
2916 /* Update vectorial force */
2917 fix2 = _mm_add_ps(fix2,tx);
2918 fiy2 = _mm_add_ps(fiy2,ty);
2919 fiz2 = _mm_add_ps(fiz2,tz);
2921 fjx0 = _mm_add_ps(fjx0,tx);
2922 fjy0 = _mm_add_ps(fjy0,ty);
2923 fjz0 = _mm_add_ps(fjz0,tz);
2927 /**************************
2928 * CALCULATE INTERACTIONS *
2929 **************************/
2931 if (gmx_mm_any_lt(rsq21,rcutoff2))
2934 r21 = _mm_mul_ps(rsq21,rinv21);
2935 r21 = _mm_andnot_ps(dummy_mask,r21);
2937 /* EWALD ELECTROSTATICS */
2939 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
2940 ewrt = _mm_mul_ps(r21,ewtabscale);
2941 ewitab = _mm_cvttps_epi32(ewrt);
2942 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
2943 ewitab = _mm_slli_epi32(ewitab,2);
2944 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
2945 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
2946 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
2947 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
2948 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
2949 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
2950 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
2951 velec = _mm_mul_ps(qq21,_mm_sub_ps(rinv21,velec));
2952 felec = _mm_mul_ps(_mm_mul_ps(qq21,rinv21),_mm_sub_ps(rinvsq21,felec));
2954 d = _mm_sub_ps(r21,rswitch);
2955 d = _mm_max_ps(d,_mm_setzero_ps());
2956 d2 = _mm_mul_ps(d,d);
2957 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
2959 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
2961 /* Evaluate switch function */
2962 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
2963 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv21,_mm_mul_ps(velec,dsw)) );
2964 cutoff_mask = _mm_cmplt_ps(rsq21,rcutoff2);
2968 fscal = _mm_and_ps(fscal,cutoff_mask);
2970 fscal = _mm_andnot_ps(dummy_mask,fscal);
2972 /* Calculate temporary vectorial force */
2973 tx = _mm_mul_ps(fscal,dx21);
2974 ty = _mm_mul_ps(fscal,dy21);
2975 tz = _mm_mul_ps(fscal,dz21);
2977 /* Update vectorial force */
2978 fix2 = _mm_add_ps(fix2,tx);
2979 fiy2 = _mm_add_ps(fiy2,ty);
2980 fiz2 = _mm_add_ps(fiz2,tz);
2982 fjx1 = _mm_add_ps(fjx1,tx);
2983 fjy1 = _mm_add_ps(fjy1,ty);
2984 fjz1 = _mm_add_ps(fjz1,tz);
2988 /**************************
2989 * CALCULATE INTERACTIONS *
2990 **************************/
2992 if (gmx_mm_any_lt(rsq22,rcutoff2))
2995 r22 = _mm_mul_ps(rsq22,rinv22);
2996 r22 = _mm_andnot_ps(dummy_mask,r22);
2998 /* EWALD ELECTROSTATICS */
3000 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
3001 ewrt = _mm_mul_ps(r22,ewtabscale);
3002 ewitab = _mm_cvttps_epi32(ewrt);
3003 eweps = _mm_sub_ps(ewrt,_mm_round_ps(ewrt, _MM_FROUND_FLOOR));
3004 ewitab = _mm_slli_epi32(ewitab,2);
3005 ewtabF = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,0) );
3006 ewtabD = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,1) );
3007 ewtabV = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,2) );
3008 ewtabFn = _mm_load_ps( ewtab + gmx_mm_extract_epi32(ewitab,3) );
3009 _MM_TRANSPOSE4_PS(ewtabF,ewtabD,ewtabV,ewtabFn);
3010 felec = _mm_add_ps(ewtabF,_mm_mul_ps(eweps,ewtabD));
3011 velec = _mm_sub_ps(ewtabV,_mm_mul_ps(_mm_mul_ps(ewtabhalfspace,eweps),_mm_add_ps(ewtabF,felec)));
3012 velec = _mm_mul_ps(qq22,_mm_sub_ps(rinv22,velec));
3013 felec = _mm_mul_ps(_mm_mul_ps(qq22,rinv22),_mm_sub_ps(rinvsq22,felec));
3015 d = _mm_sub_ps(r22,rswitch);
3016 d = _mm_max_ps(d,_mm_setzero_ps());
3017 d2 = _mm_mul_ps(d,d);
3018 sw = _mm_add_ps(one,_mm_mul_ps(d2,_mm_mul_ps(d,_mm_add_ps(swV3,_mm_mul_ps(d,_mm_add_ps(swV4,_mm_mul_ps(d,swV5)))))));
3020 dsw = _mm_mul_ps(d2,_mm_add_ps(swF2,_mm_mul_ps(d,_mm_add_ps(swF3,_mm_mul_ps(d,swF4)))));
3022 /* Evaluate switch function */
3023 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
3024 felec = _mm_sub_ps( _mm_mul_ps(felec,sw) , _mm_mul_ps(rinv22,_mm_mul_ps(velec,dsw)) );
3025 cutoff_mask = _mm_cmplt_ps(rsq22,rcutoff2);
3029 fscal = _mm_and_ps(fscal,cutoff_mask);
3031 fscal = _mm_andnot_ps(dummy_mask,fscal);
3033 /* Calculate temporary vectorial force */
3034 tx = _mm_mul_ps(fscal,dx22);
3035 ty = _mm_mul_ps(fscal,dy22);
3036 tz = _mm_mul_ps(fscal,dz22);
3038 /* Update vectorial force */
3039 fix2 = _mm_add_ps(fix2,tx);
3040 fiy2 = _mm_add_ps(fiy2,ty);
3041 fiz2 = _mm_add_ps(fiz2,tz);
3043 fjx2 = _mm_add_ps(fjx2,tx);
3044 fjy2 = _mm_add_ps(fjy2,ty);
3045 fjz2 = _mm_add_ps(fjz2,tz);
3049 fjptrA = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
3050 fjptrB = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
3051 fjptrC = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
3052 fjptrD = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
3054 gmx_mm_decrement_3rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,
3055 fjx0,fjy0,fjz0,fjx1,fjy1,fjz1,fjx2,fjy2,fjz2);
3057 /* Inner loop uses 567 flops */
3060 /* End of innermost loop */
3062 gmx_mm_update_iforce_3atom_swizzle_ps(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,
3063 f+i_coord_offset,fshift+i_shift_offset);
3065 /* Increment number of inner iterations */
3066 inneriter += j_index_end - j_index_start;
3068 /* Outer loop uses 18 flops */
3071 /* Increment number of outer iterations */
3074 /* Update outer/inner flops */
3076 inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_W3W3_F,outeriter*18 + inneriter*567);