?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Remove all unnecessary HAVE_CONFIG_H
[gromacs.git] / src / gromacs / gmxlib / nonbonded / nb_kernel_sse4_1_double / nb_kernel_ElecEw_VdwCSTab_GeomW4P1_sse4_1_double.c
blobaf684cd223a19001cd8c1b4352cfa2236043b166
1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
8 *
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.
13 *
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.
18 *
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.
23 *
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.
31 *
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.
34 */
35 /*
36 * Note: this file was generated by the GROMACS sse4_1_double kernel generator.
37 */
38 #include "config.h"
39
40 #include <math.h>
41
42 #include "../nb_kernel.h"
43 #include "types/simple.h"
44 #include "gromacs/math/vec.h"
45 #include "nrnb.h"
46
47 #include "gromacs/simd/math_x86_sse4_1_double.h"
48 #include "kernelutil_x86_sse4_1_double.h"
49
50 /*
51 * Gromacs nonbonded kernel: nb_kernel_ElecEw_VdwCSTab_GeomW4P1_VF_sse4_1_double
52 * Electrostatics interaction: Ewald
53 * VdW interaction: CubicSplineTable
54 * Geometry: Water4-Particle
55 * Calculate force/pot: PotentialAndForce
56 */
57 void
58 nb_kernel_ElecEw_VdwCSTab_GeomW4P1_VF_sse4_1_double
59 (t_nblist * gmx_restrict nlist,
60 rvec * gmx_restrict xx,
61 rvec * gmx_restrict ff,
62 t_forcerec * gmx_restrict fr,
63 t_mdatoms * gmx_restrict mdatoms,
64 nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
65 t_nrnb * gmx_restrict nrnb)
66 {
67 /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
68 * just 0 for non-waters.
69 * Suffixes A,B refer to j loop unrolling done with SSE double precision, e.g. for the two different
70 * jnr indices corresponding to data put in the four positions in the SIMD register.
71 */
72 int i_shift_offset,i_coord_offset,outeriter,inneriter;
73 int j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
74 int jnrA,jnrB;
75 int j_coord_offsetA,j_coord_offsetB;
76 int *iinr,*jindex,*jjnr,*shiftidx,*gid;
77 real rcutoff_scalar;
78 real *shiftvec,*fshift,*x,*f;
79 __m128d tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
80 int vdwioffset0;
81 __m128d ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
82 int vdwioffset1;
83 __m128d ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
84 int vdwioffset2;
85 __m128d ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
86 int vdwioffset3;
87 __m128d ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
88 int vdwjidx0A,vdwjidx0B;
89 __m128d jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
90 __m128d dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
91 __m128d dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
92 __m128d dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
93 __m128d dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
94 __m128d velec,felec,velecsum,facel,crf,krf,krf2;
95 real *charge;
96 int nvdwtype;
97 __m128d rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
98 int *vdwtype;
99 real *vdwparam;
100 __m128d one_sixth = _mm_set1_pd(1.0/6.0);
101 __m128d one_twelfth = _mm_set1_pd(1.0/12.0);
102 __m128i vfitab;
103 __m128i ifour = _mm_set1_epi32(4);
104 __m128d rt,vfeps,vftabscale,Y,F,G,H,Heps,Fp,VV,FF;
105 real *vftab;
106 __m128i ewitab;
107 __m128d ewtabscale,eweps,sh_ewald,ewrt,ewtabhalfspace,ewtabF,ewtabFn,ewtabD,ewtabV;
108 real *ewtab;
109 __m128d dummy_mask,cutoff_mask;
110 __m128d signbit = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
111 __m128d one = _mm_set1_pd(1.0);
112 __m128d two = _mm_set1_pd(2.0);
113 x = xx[0];
114 f = ff[0];
115
116 nri = nlist->nri;
117 iinr = nlist->iinr;
118 jindex = nlist->jindex;
119 jjnr = nlist->jjnr;
120 shiftidx = nlist->shift;
121 gid = nlist->gid;
122 shiftvec = fr->shift_vec[0];
123 fshift = fr->fshift[0];
124 facel = _mm_set1_pd(fr->epsfac);
125 charge = mdatoms->chargeA;
126 nvdwtype = fr->ntype;
127 vdwparam = fr->nbfp;
128 vdwtype = mdatoms->typeA;
129
130 vftab = kernel_data->table_vdw->data;
131 vftabscale = _mm_set1_pd(kernel_data->table_vdw->scale);
132
133 sh_ewald = _mm_set1_pd(fr->ic->sh_ewald);
134 ewtab = fr->ic->tabq_coul_FDV0;
135 ewtabscale = _mm_set1_pd(fr->ic->tabq_scale);
136 ewtabhalfspace = _mm_set1_pd(0.5/fr->ic->tabq_scale);
137
138 /* Setup water-specific parameters */
139 inr = nlist->iinr[0];
140 iq1 = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+1]));
141 iq2 = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+2]));
142 iq3 = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+3]));
143 vdwioffset0 = 2*nvdwtype*vdwtype[inr+0];
144
145 /* Avoid stupid compiler warnings */
146 jnrA = jnrB = 0;
147 j_coord_offsetA = 0;
148 j_coord_offsetB = 0;
149
150 outeriter = 0;
151 inneriter = 0;
152
153 /* Start outer loop over neighborlists */
154 for(iidx=0; iidx<nri; iidx++)
155 {
156 /* Load shift vector for this list */
157 i_shift_offset = DIM*shiftidx[iidx];
158
159 /* Load limits for loop over neighbors */
160 j_index_start = jindex[iidx];
161 j_index_end = jindex[iidx+1];
162
163 /* Get outer coordinate index */
164 inr = iinr[iidx];
165 i_coord_offset = DIM*inr;
166
167 /* Load i particle coords and add shift vector */
168 gmx_mm_load_shift_and_4rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,
169 &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
170
171 fix0 = _mm_setzero_pd();
172 fiy0 = _mm_setzero_pd();
173 fiz0 = _mm_setzero_pd();
174 fix1 = _mm_setzero_pd();
175 fiy1 = _mm_setzero_pd();
176 fiz1 = _mm_setzero_pd();
177 fix2 = _mm_setzero_pd();
178 fiy2 = _mm_setzero_pd();
179 fiz2 = _mm_setzero_pd();
180 fix3 = _mm_setzero_pd();
181 fiy3 = _mm_setzero_pd();
182 fiz3 = _mm_setzero_pd();
183
184 /* Reset potential sums */
185 velecsum = _mm_setzero_pd();
186 vvdwsum = _mm_setzero_pd();
187
188 /* Start inner kernel loop */
189 for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
190 {
191
192 /* Get j neighbor index, and coordinate index */
193 jnrA = jjnr[jidx];
194 jnrB = jjnr[jidx+1];
195 j_coord_offsetA = DIM*jnrA;
196 j_coord_offsetB = DIM*jnrB;
197
198 /* load j atom coordinates */
199 gmx_mm_load_1rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
200 &jx0,&jy0,&jz0);
201
202 /* Calculate displacement vector */
203 dx00 = _mm_sub_pd(ix0,jx0);
204 dy00 = _mm_sub_pd(iy0,jy0);
205 dz00 = _mm_sub_pd(iz0,jz0);
206 dx10 = _mm_sub_pd(ix1,jx0);
207 dy10 = _mm_sub_pd(iy1,jy0);
208 dz10 = _mm_sub_pd(iz1,jz0);
209 dx20 = _mm_sub_pd(ix2,jx0);
210 dy20 = _mm_sub_pd(iy2,jy0);
211 dz20 = _mm_sub_pd(iz2,jz0);
212 dx30 = _mm_sub_pd(ix3,jx0);
213 dy30 = _mm_sub_pd(iy3,jy0);
214 dz30 = _mm_sub_pd(iz3,jz0);
215
216 /* Calculate squared distance and things based on it */
217 rsq00 = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
218 rsq10 = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
219 rsq20 = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
220 rsq30 = gmx_mm_calc_rsq_pd(dx30,dy30,dz30);
221
222 rinv00 = gmx_mm_invsqrt_pd(rsq00);
223 rinv10 = gmx_mm_invsqrt_pd(rsq10);
224 rinv20 = gmx_mm_invsqrt_pd(rsq20);
225 rinv30 = gmx_mm_invsqrt_pd(rsq30);
226
227 rinvsq10 = _mm_mul_pd(rinv10,rinv10);
228 rinvsq20 = _mm_mul_pd(rinv20,rinv20);
229 rinvsq30 = _mm_mul_pd(rinv30,rinv30);
230
231 /* Load parameters for j particles */
232 jq0 = gmx_mm_load_2real_swizzle_pd(charge+jnrA+0,charge+jnrB+0);
233 vdwjidx0A = 2*vdwtype[jnrA+0];
234 vdwjidx0B = 2*vdwtype[jnrB+0];
235
236 fjx0 = _mm_setzero_pd();
237 fjy0 = _mm_setzero_pd();
238 fjz0 = _mm_setzero_pd();
239
240 /**************************
241 * CALCULATE INTERACTIONS *
242 **************************/
243
244 r00 = _mm_mul_pd(rsq00,rinv00);
245
246 /* Compute parameters for interactions between i and j atoms */
247 gmx_mm_load_2pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,
248 vdwparam+vdwioffset0+vdwjidx0B,&c6_00,&c12_00);
249
250 /* Calculate table index by multiplying r with table scale and truncate to integer */
251 rt = _mm_mul_pd(r00,vftabscale);
252 vfitab = _mm_cvttpd_epi32(rt);
253 vfeps = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
254 vfitab = _mm_slli_epi32(vfitab,3);
255
256 /* CUBIC SPLINE TABLE DISPERSION */
257 Y = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) );
258 F = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,1) );
259 GMX_MM_TRANSPOSE2_PD(Y,F);
260 G = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) +2);
261 H = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,1) +2);
262 GMX_MM_TRANSPOSE2_PD(G,H);
263 Heps = _mm_mul_pd(vfeps,H);
264 Fp = _mm_add_pd(F,_mm_mul_pd(vfeps,_mm_add_pd(G,Heps)));
265 VV = _mm_add_pd(Y,_mm_mul_pd(vfeps,Fp));
266 vvdw6 = _mm_mul_pd(c6_00,VV);
267 FF = _mm_add_pd(Fp,_mm_mul_pd(vfeps,_mm_add_pd(G,_mm_add_pd(Heps,Heps))));
268 fvdw6 = _mm_mul_pd(c6_00,FF);
269
270 /* CUBIC SPLINE TABLE REPULSION */
271 vfitab = _mm_add_epi32(vfitab,ifour);
272 Y = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) );
273 F = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,1) );
274 GMX_MM_TRANSPOSE2_PD(Y,F);
275 G = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) +2);
276 H = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,1) +2);
277 GMX_MM_TRANSPOSE2_PD(G,H);
278 Heps = _mm_mul_pd(vfeps,H);
279 Fp = _mm_add_pd(F,_mm_mul_pd(vfeps,_mm_add_pd(G,Heps)));
280 VV = _mm_add_pd(Y,_mm_mul_pd(vfeps,Fp));
281 vvdw12 = _mm_mul_pd(c12_00,VV);
282 FF = _mm_add_pd(Fp,_mm_mul_pd(vfeps,_mm_add_pd(G,_mm_add_pd(Heps,Heps))));
283 fvdw12 = _mm_mul_pd(c12_00,FF);
284 vvdw = _mm_add_pd(vvdw12,vvdw6);
285 fvdw = _mm_xor_pd(signbit,_mm_mul_pd(_mm_add_pd(fvdw6,fvdw12),_mm_mul_pd(vftabscale,rinv00)));
286
287 /* Update potential sum for this i atom from the interaction with this j atom. */
288 vvdwsum = _mm_add_pd(vvdwsum,vvdw);
289
290 fscal = fvdw;
291
292 /* Calculate temporary vectorial force */
293 tx = _mm_mul_pd(fscal,dx00);
294 ty = _mm_mul_pd(fscal,dy00);
295 tz = _mm_mul_pd(fscal,dz00);
296
297 /* Update vectorial force */
298 fix0 = _mm_add_pd(fix0,tx);
299 fiy0 = _mm_add_pd(fiy0,ty);
300 fiz0 = _mm_add_pd(fiz0,tz);
301
302 fjx0 = _mm_add_pd(fjx0,tx);
303 fjy0 = _mm_add_pd(fjy0,ty);
304 fjz0 = _mm_add_pd(fjz0,tz);
305
306 /**************************
307 * CALCULATE INTERACTIONS *
308 **************************/
309
310 r10 = _mm_mul_pd(rsq10,rinv10);
311
312 /* Compute parameters for interactions between i and j atoms */
313 qq10 = _mm_mul_pd(iq1,jq0);
314
315 /* EWALD ELECTROSTATICS */
316
317 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
318 ewrt = _mm_mul_pd(r10,ewtabscale);
319 ewitab = _mm_cvttpd_epi32(ewrt);
320 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
321 ewitab = _mm_slli_epi32(ewitab,2);
322 ewtabF = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
323 ewtabD = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
324 GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
325 ewtabV = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
326 ewtabFn = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
327 GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
328 felec = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
329 velec = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
330 velec = _mm_mul_pd(qq10,_mm_sub_pd(rinv10,velec));
331 felec = _mm_mul_pd(_mm_mul_pd(qq10,rinv10),_mm_sub_pd(rinvsq10,felec));
332
333 /* Update potential sum for this i atom from the interaction with this j atom. */
334 velecsum = _mm_add_pd(velecsum,velec);
335
336 fscal = felec;
337
338 /* Calculate temporary vectorial force */
339 tx = _mm_mul_pd(fscal,dx10);
340 ty = _mm_mul_pd(fscal,dy10);
341 tz = _mm_mul_pd(fscal,dz10);
342
343 /* Update vectorial force */
344 fix1 = _mm_add_pd(fix1,tx);
345 fiy1 = _mm_add_pd(fiy1,ty);
346 fiz1 = _mm_add_pd(fiz1,tz);
347
348 fjx0 = _mm_add_pd(fjx0,tx);
349 fjy0 = _mm_add_pd(fjy0,ty);
350 fjz0 = _mm_add_pd(fjz0,tz);
351
352 /**************************
353 * CALCULATE INTERACTIONS *
354 **************************/
355
356 r20 = _mm_mul_pd(rsq20,rinv20);
357
358 /* Compute parameters for interactions between i and j atoms */
359 qq20 = _mm_mul_pd(iq2,jq0);
360
361 /* EWALD ELECTROSTATICS */
362
363 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
364 ewrt = _mm_mul_pd(r20,ewtabscale);
365 ewitab = _mm_cvttpd_epi32(ewrt);
366 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
367 ewitab = _mm_slli_epi32(ewitab,2);
368 ewtabF = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
369 ewtabD = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
370 GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
371 ewtabV = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
372 ewtabFn = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
373 GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
374 felec = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
375 velec = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
376 velec = _mm_mul_pd(qq20,_mm_sub_pd(rinv20,velec));
377 felec = _mm_mul_pd(_mm_mul_pd(qq20,rinv20),_mm_sub_pd(rinvsq20,felec));
378
379 /* Update potential sum for this i atom from the interaction with this j atom. */
380 velecsum = _mm_add_pd(velecsum,velec);
381
382 fscal = felec;
383
384 /* Calculate temporary vectorial force */
385 tx = _mm_mul_pd(fscal,dx20);
386 ty = _mm_mul_pd(fscal,dy20);
387 tz = _mm_mul_pd(fscal,dz20);
388
389 /* Update vectorial force */
390 fix2 = _mm_add_pd(fix2,tx);
391 fiy2 = _mm_add_pd(fiy2,ty);
392 fiz2 = _mm_add_pd(fiz2,tz);
393
394 fjx0 = _mm_add_pd(fjx0,tx);
395 fjy0 = _mm_add_pd(fjy0,ty);
396 fjz0 = _mm_add_pd(fjz0,tz);
397
398 /**************************
399 * CALCULATE INTERACTIONS *
400 **************************/
401
402 r30 = _mm_mul_pd(rsq30,rinv30);
403
404 /* Compute parameters for interactions between i and j atoms */
405 qq30 = _mm_mul_pd(iq3,jq0);
406
407 /* EWALD ELECTROSTATICS */
408
409 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
410 ewrt = _mm_mul_pd(r30,ewtabscale);
411 ewitab = _mm_cvttpd_epi32(ewrt);
412 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
413 ewitab = _mm_slli_epi32(ewitab,2);
414 ewtabF = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
415 ewtabD = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,1) );
416 GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
417 ewtabV = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
418 ewtabFn = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,1) +2);
419 GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
420 felec = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
421 velec = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
422 velec = _mm_mul_pd(qq30,_mm_sub_pd(rinv30,velec));
423 felec = _mm_mul_pd(_mm_mul_pd(qq30,rinv30),_mm_sub_pd(rinvsq30,felec));
424
425 /* Update potential sum for this i atom from the interaction with this j atom. */
426 velecsum = _mm_add_pd(velecsum,velec);
427
428 fscal = felec;
429
430 /* Calculate temporary vectorial force */
431 tx = _mm_mul_pd(fscal,dx30);
432 ty = _mm_mul_pd(fscal,dy30);
433 tz = _mm_mul_pd(fscal,dz30);
434
435 /* Update vectorial force */
436 fix3 = _mm_add_pd(fix3,tx);
437 fiy3 = _mm_add_pd(fiy3,ty);
438 fiz3 = _mm_add_pd(fiz3,tz);
439
440 fjx0 = _mm_add_pd(fjx0,tx);
441 fjy0 = _mm_add_pd(fjy0,ty);
442 fjz0 = _mm_add_pd(fjz0,tz);
443
444 gmx_mm_decrement_1rvec_2ptr_swizzle_pd(f+j_coord_offsetA,f+j_coord_offsetB,fjx0,fjy0,fjz0);
445
446 /* Inner loop uses 182 flops */
447 }
448
449 if(jidx<j_index_end)
450 {
451
452 jnrA = jjnr[jidx];
453 j_coord_offsetA = DIM*jnrA;
454
455 /* load j atom coordinates */
456 gmx_mm_load_1rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
457 &jx0,&jy0,&jz0);
458
459 /* Calculate displacement vector */
460 dx00 = _mm_sub_pd(ix0,jx0);
461 dy00 = _mm_sub_pd(iy0,jy0);
462 dz00 = _mm_sub_pd(iz0,jz0);
463 dx10 = _mm_sub_pd(ix1,jx0);
464 dy10 = _mm_sub_pd(iy1,jy0);
465 dz10 = _mm_sub_pd(iz1,jz0);
466 dx20 = _mm_sub_pd(ix2,jx0);
467 dy20 = _mm_sub_pd(iy2,jy0);
468 dz20 = _mm_sub_pd(iz2,jz0);
469 dx30 = _mm_sub_pd(ix3,jx0);
470 dy30 = _mm_sub_pd(iy3,jy0);
471 dz30 = _mm_sub_pd(iz3,jz0);
472
473 /* Calculate squared distance and things based on it */
474 rsq00 = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
475 rsq10 = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
476 rsq20 = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
477 rsq30 = gmx_mm_calc_rsq_pd(dx30,dy30,dz30);
478
479 rinv00 = gmx_mm_invsqrt_pd(rsq00);
480 rinv10 = gmx_mm_invsqrt_pd(rsq10);
481 rinv20 = gmx_mm_invsqrt_pd(rsq20);
482 rinv30 = gmx_mm_invsqrt_pd(rsq30);
483
484 rinvsq10 = _mm_mul_pd(rinv10,rinv10);
485 rinvsq20 = _mm_mul_pd(rinv20,rinv20);
486 rinvsq30 = _mm_mul_pd(rinv30,rinv30);
487
488 /* Load parameters for j particles */
489 jq0 = _mm_load_sd(charge+jnrA+0);
490 vdwjidx0A = 2*vdwtype[jnrA+0];
491
492 fjx0 = _mm_setzero_pd();
493 fjy0 = _mm_setzero_pd();
494 fjz0 = _mm_setzero_pd();
495
496 /**************************
497 * CALCULATE INTERACTIONS *
498 **************************/
499
500 r00 = _mm_mul_pd(rsq00,rinv00);
501
502 /* Compute parameters for interactions between i and j atoms */
503 gmx_mm_load_1pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,&c6_00,&c12_00);
504
505 /* Calculate table index by multiplying r with table scale and truncate to integer */
506 rt = _mm_mul_pd(r00,vftabscale);
507 vfitab = _mm_cvttpd_epi32(rt);
508 vfeps = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
509 vfitab = _mm_slli_epi32(vfitab,3);
510
511 /* CUBIC SPLINE TABLE DISPERSION */
512 Y = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) );
513 F = _mm_setzero_pd();
514 GMX_MM_TRANSPOSE2_PD(Y,F);
515 G = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) +2);
516 H = _mm_setzero_pd();
517 GMX_MM_TRANSPOSE2_PD(G,H);
518 Heps = _mm_mul_pd(vfeps,H);
519 Fp = _mm_add_pd(F,_mm_mul_pd(vfeps,_mm_add_pd(G,Heps)));
520 VV = _mm_add_pd(Y,_mm_mul_pd(vfeps,Fp));
521 vvdw6 = _mm_mul_pd(c6_00,VV);
522 FF = _mm_add_pd(Fp,_mm_mul_pd(vfeps,_mm_add_pd(G,_mm_add_pd(Heps,Heps))));
523 fvdw6 = _mm_mul_pd(c6_00,FF);
524
525 /* CUBIC SPLINE TABLE REPULSION */
526 vfitab = _mm_add_epi32(vfitab,ifour);
527 Y = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) );
528 F = _mm_setzero_pd();
529 GMX_MM_TRANSPOSE2_PD(Y,F);
530 G = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) +2);
531 H = _mm_setzero_pd();
532 GMX_MM_TRANSPOSE2_PD(G,H);
533 Heps = _mm_mul_pd(vfeps,H);
534 Fp = _mm_add_pd(F,_mm_mul_pd(vfeps,_mm_add_pd(G,Heps)));
535 VV = _mm_add_pd(Y,_mm_mul_pd(vfeps,Fp));
536 vvdw12 = _mm_mul_pd(c12_00,VV);
537 FF = _mm_add_pd(Fp,_mm_mul_pd(vfeps,_mm_add_pd(G,_mm_add_pd(Heps,Heps))));
538 fvdw12 = _mm_mul_pd(c12_00,FF);
539 vvdw = _mm_add_pd(vvdw12,vvdw6);
540 fvdw = _mm_xor_pd(signbit,_mm_mul_pd(_mm_add_pd(fvdw6,fvdw12),_mm_mul_pd(vftabscale,rinv00)));
541
542 /* Update potential sum for this i atom from the interaction with this j atom. */
543 vvdw = _mm_unpacklo_pd(vvdw,_mm_setzero_pd());
544 vvdwsum = _mm_add_pd(vvdwsum,vvdw);
545
546 fscal = fvdw;
547
548 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
549
550 /* Calculate temporary vectorial force */
551 tx = _mm_mul_pd(fscal,dx00);
552 ty = _mm_mul_pd(fscal,dy00);
553 tz = _mm_mul_pd(fscal,dz00);
554
555 /* Update vectorial force */
556 fix0 = _mm_add_pd(fix0,tx);
557 fiy0 = _mm_add_pd(fiy0,ty);
558 fiz0 = _mm_add_pd(fiz0,tz);
559
560 fjx0 = _mm_add_pd(fjx0,tx);
561 fjy0 = _mm_add_pd(fjy0,ty);
562 fjz0 = _mm_add_pd(fjz0,tz);
563
564 /**************************
565 * CALCULATE INTERACTIONS *
566 **************************/
567
568 r10 = _mm_mul_pd(rsq10,rinv10);
569
570 /* Compute parameters for interactions between i and j atoms */
571 qq10 = _mm_mul_pd(iq1,jq0);
572
573 /* EWALD ELECTROSTATICS */
574
575 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
576 ewrt = _mm_mul_pd(r10,ewtabscale);
577 ewitab = _mm_cvttpd_epi32(ewrt);
578 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
579 ewitab = _mm_slli_epi32(ewitab,2);
580 ewtabF = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
581 ewtabD = _mm_setzero_pd();
582 GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
583 ewtabV = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
584 ewtabFn = _mm_setzero_pd();
585 GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
586 felec = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
587 velec = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
588 velec = _mm_mul_pd(qq10,_mm_sub_pd(rinv10,velec));
589 felec = _mm_mul_pd(_mm_mul_pd(qq10,rinv10),_mm_sub_pd(rinvsq10,felec));
590
591 /* Update potential sum for this i atom from the interaction with this j atom. */
592 velec = _mm_unpacklo_pd(velec,_mm_setzero_pd());
593 velecsum = _mm_add_pd(velecsum,velec);
594
595 fscal = felec;
596
597 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
598
599 /* Calculate temporary vectorial force */
600 tx = _mm_mul_pd(fscal,dx10);
601 ty = _mm_mul_pd(fscal,dy10);
602 tz = _mm_mul_pd(fscal,dz10);
603
604 /* Update vectorial force */
605 fix1 = _mm_add_pd(fix1,tx);
606 fiy1 = _mm_add_pd(fiy1,ty);
607 fiz1 = _mm_add_pd(fiz1,tz);
608
609 fjx0 = _mm_add_pd(fjx0,tx);
610 fjy0 = _mm_add_pd(fjy0,ty);
611 fjz0 = _mm_add_pd(fjz0,tz);
612
613 /**************************
614 * CALCULATE INTERACTIONS *
615 **************************/
616
617 r20 = _mm_mul_pd(rsq20,rinv20);
618
619 /* Compute parameters for interactions between i and j atoms */
620 qq20 = _mm_mul_pd(iq2,jq0);
621
622 /* EWALD ELECTROSTATICS */
623
624 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
625 ewrt = _mm_mul_pd(r20,ewtabscale);
626 ewitab = _mm_cvttpd_epi32(ewrt);
627 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
628 ewitab = _mm_slli_epi32(ewitab,2);
629 ewtabF = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
630 ewtabD = _mm_setzero_pd();
631 GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
632 ewtabV = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
633 ewtabFn = _mm_setzero_pd();
634 GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
635 felec = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
636 velec = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
637 velec = _mm_mul_pd(qq20,_mm_sub_pd(rinv20,velec));
638 felec = _mm_mul_pd(_mm_mul_pd(qq20,rinv20),_mm_sub_pd(rinvsq20,felec));
639
640 /* Update potential sum for this i atom from the interaction with this j atom. */
641 velec = _mm_unpacklo_pd(velec,_mm_setzero_pd());
642 velecsum = _mm_add_pd(velecsum,velec);
643
644 fscal = felec;
645
646 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
647
648 /* Calculate temporary vectorial force */
649 tx = _mm_mul_pd(fscal,dx20);
650 ty = _mm_mul_pd(fscal,dy20);
651 tz = _mm_mul_pd(fscal,dz20);
652
653 /* Update vectorial force */
654 fix2 = _mm_add_pd(fix2,tx);
655 fiy2 = _mm_add_pd(fiy2,ty);
656 fiz2 = _mm_add_pd(fiz2,tz);
657
658 fjx0 = _mm_add_pd(fjx0,tx);
659 fjy0 = _mm_add_pd(fjy0,ty);
660 fjz0 = _mm_add_pd(fjz0,tz);
661
662 /**************************
663 * CALCULATE INTERACTIONS *
664 **************************/
665
666 r30 = _mm_mul_pd(rsq30,rinv30);
667
668 /* Compute parameters for interactions between i and j atoms */
669 qq30 = _mm_mul_pd(iq3,jq0);
670
671 /* EWALD ELECTROSTATICS */
672
673 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
674 ewrt = _mm_mul_pd(r30,ewtabscale);
675 ewitab = _mm_cvttpd_epi32(ewrt);
676 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
677 ewitab = _mm_slli_epi32(ewitab,2);
678 ewtabF = _mm_load_pd( ewtab + gmx_mm_extract_epi32(ewitab,0) );
679 ewtabD = _mm_setzero_pd();
680 GMX_MM_TRANSPOSE2_PD(ewtabF,ewtabD);
681 ewtabV = _mm_load_sd( ewtab + gmx_mm_extract_epi32(ewitab,0) +2);
682 ewtabFn = _mm_setzero_pd();
683 GMX_MM_TRANSPOSE2_PD(ewtabV,ewtabFn);
684 felec = _mm_add_pd(ewtabF,_mm_mul_pd(eweps,ewtabD));
685 velec = _mm_sub_pd(ewtabV,_mm_mul_pd(_mm_mul_pd(ewtabhalfspace,eweps),_mm_add_pd(ewtabF,felec)));
686 velec = _mm_mul_pd(qq30,_mm_sub_pd(rinv30,velec));
687 felec = _mm_mul_pd(_mm_mul_pd(qq30,rinv30),_mm_sub_pd(rinvsq30,felec));
688
689 /* Update potential sum for this i atom from the interaction with this j atom. */
690 velec = _mm_unpacklo_pd(velec,_mm_setzero_pd());
691 velecsum = _mm_add_pd(velecsum,velec);
692
693 fscal = felec;
694
695 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
696
697 /* Calculate temporary vectorial force */
698 tx = _mm_mul_pd(fscal,dx30);
699 ty = _mm_mul_pd(fscal,dy30);
700 tz = _mm_mul_pd(fscal,dz30);
701
702 /* Update vectorial force */
703 fix3 = _mm_add_pd(fix3,tx);
704 fiy3 = _mm_add_pd(fiy3,ty);
705 fiz3 = _mm_add_pd(fiz3,tz);
706
707 fjx0 = _mm_add_pd(fjx0,tx);
708 fjy0 = _mm_add_pd(fjy0,ty);
709 fjz0 = _mm_add_pd(fjz0,tz);
710
711 gmx_mm_decrement_1rvec_1ptr_swizzle_pd(f+j_coord_offsetA,fjx0,fjy0,fjz0);
712
713 /* Inner loop uses 182 flops */
714 }
715
716 /* End of innermost loop */
717
718 gmx_mm_update_iforce_4atom_swizzle_pd(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
719 f+i_coord_offset,fshift+i_shift_offset);
720
721 ggid = gid[iidx];
722 /* Update potential energies */
723 gmx_mm_update_1pot_pd(velecsum,kernel_data->energygrp_elec+ggid);
724 gmx_mm_update_1pot_pd(vvdwsum,kernel_data->energygrp_vdw+ggid);
725
726 /* Increment number of inner iterations */
727 inneriter += j_index_end - j_index_start;
728
729 /* Outer loop uses 26 flops */
730 }
731
732 /* Increment number of outer iterations */
733 outeriter += nri;
734
735 /* Update outer/inner flops */
736
737 inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4_VF,outeriter*26 + inneriter*182);
738 }
739 /*
740 * Gromacs nonbonded kernel: nb_kernel_ElecEw_VdwCSTab_GeomW4P1_F_sse4_1_double
741 * Electrostatics interaction: Ewald
742 * VdW interaction: CubicSplineTable
743 * Geometry: Water4-Particle
744 * Calculate force/pot: Force
745 */
746 void
747 nb_kernel_ElecEw_VdwCSTab_GeomW4P1_F_sse4_1_double
748 (t_nblist * gmx_restrict nlist,
749 rvec * gmx_restrict xx,
750 rvec * gmx_restrict ff,
751 t_forcerec * gmx_restrict fr,
752 t_mdatoms * gmx_restrict mdatoms,
753 nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
754 t_nrnb * gmx_restrict nrnb)
755 {
756 /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or
757 * just 0 for non-waters.
758 * Suffixes A,B refer to j loop unrolling done with SSE double precision, e.g. for the two different
759 * jnr indices corresponding to data put in the four positions in the SIMD register.
760 */
761 int i_shift_offset,i_coord_offset,outeriter,inneriter;
762 int j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
763 int jnrA,jnrB;
764 int j_coord_offsetA,j_coord_offsetB;
765 int *iinr,*jindex,*jjnr,*shiftidx,*gid;
766 real rcutoff_scalar;
767 real *shiftvec,*fshift,*x,*f;
768 __m128d tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
769 int vdwioffset0;
770 __m128d ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
771 int vdwioffset1;
772 __m128d ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
773 int vdwioffset2;
774 __m128d ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
775 int vdwioffset3;
776 __m128d ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
777 int vdwjidx0A,vdwjidx0B;
778 __m128d jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
779 __m128d dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
780 __m128d dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
781 __m128d dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
782 __m128d dx30,dy30,dz30,rsq30,rinv30,rinvsq30,r30,qq30,c6_30,c12_30;
783 __m128d velec,felec,velecsum,facel,crf,krf,krf2;
784 real *charge;
785 int nvdwtype;
786 __m128d rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
787 int *vdwtype;
788 real *vdwparam;
789 __m128d one_sixth = _mm_set1_pd(1.0/6.0);
790 __m128d one_twelfth = _mm_set1_pd(1.0/12.0);
791 __m128i vfitab;
792 __m128i ifour = _mm_set1_epi32(4);
793 __m128d rt,vfeps,vftabscale,Y,F,G,H,Heps,Fp,VV,FF;
794 real *vftab;
795 __m128i ewitab;
796 __m128d ewtabscale,eweps,sh_ewald,ewrt,ewtabhalfspace,ewtabF,ewtabFn,ewtabD,ewtabV;
797 real *ewtab;
798 __m128d dummy_mask,cutoff_mask;
799 __m128d signbit = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
800 __m128d one = _mm_set1_pd(1.0);
801 __m128d two = _mm_set1_pd(2.0);
802 x = xx[0];
803 f = ff[0];
804
805 nri = nlist->nri;
806 iinr = nlist->iinr;
807 jindex = nlist->jindex;
808 jjnr = nlist->jjnr;
809 shiftidx = nlist->shift;
810 gid = nlist->gid;
811 shiftvec = fr->shift_vec[0];
812 fshift = fr->fshift[0];
813 facel = _mm_set1_pd(fr->epsfac);
814 charge = mdatoms->chargeA;
815 nvdwtype = fr->ntype;
816 vdwparam = fr->nbfp;
817 vdwtype = mdatoms->typeA;
818
819 vftab = kernel_data->table_vdw->data;
820 vftabscale = _mm_set1_pd(kernel_data->table_vdw->scale);
821
822 sh_ewald = _mm_set1_pd(fr->ic->sh_ewald);
823 ewtab = fr->ic->tabq_coul_F;
824 ewtabscale = _mm_set1_pd(fr->ic->tabq_scale);
825 ewtabhalfspace = _mm_set1_pd(0.5/fr->ic->tabq_scale);
826
827 /* Setup water-specific parameters */
828 inr = nlist->iinr[0];
829 iq1 = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+1]));
830 iq2 = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+2]));
831 iq3 = _mm_mul_pd(facel,_mm_set1_pd(charge[inr+3]));
832 vdwioffset0 = 2*nvdwtype*vdwtype[inr+0];
833
834 /* Avoid stupid compiler warnings */
835 jnrA = jnrB = 0;
836 j_coord_offsetA = 0;
837 j_coord_offsetB = 0;
838
839 outeriter = 0;
840 inneriter = 0;
841
842 /* Start outer loop over neighborlists */
843 for(iidx=0; iidx<nri; iidx++)
844 {
845 /* Load shift vector for this list */
846 i_shift_offset = DIM*shiftidx[iidx];
847
848 /* Load limits for loop over neighbors */
849 j_index_start = jindex[iidx];
850 j_index_end = jindex[iidx+1];
851
852 /* Get outer coordinate index */
853 inr = iinr[iidx];
854 i_coord_offset = DIM*inr;
855
856 /* Load i particle coords and add shift vector */
857 gmx_mm_load_shift_and_4rvec_broadcast_pd(shiftvec+i_shift_offset,x+i_coord_offset,
858 &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2,&ix3,&iy3,&iz3);
859
860 fix0 = _mm_setzero_pd();
861 fiy0 = _mm_setzero_pd();
862 fiz0 = _mm_setzero_pd();
863 fix1 = _mm_setzero_pd();
864 fiy1 = _mm_setzero_pd();
865 fiz1 = _mm_setzero_pd();
866 fix2 = _mm_setzero_pd();
867 fiy2 = _mm_setzero_pd();
868 fiz2 = _mm_setzero_pd();
869 fix3 = _mm_setzero_pd();
870 fiy3 = _mm_setzero_pd();
871 fiz3 = _mm_setzero_pd();
872
873 /* Start inner kernel loop */
874 for(jidx=j_index_start; jidx<j_index_end-1; jidx+=2)
875 {
876
877 /* Get j neighbor index, and coordinate index */
878 jnrA = jjnr[jidx];
879 jnrB = jjnr[jidx+1];
880 j_coord_offsetA = DIM*jnrA;
881 j_coord_offsetB = DIM*jnrB;
882
883 /* load j atom coordinates */
884 gmx_mm_load_1rvec_2ptr_swizzle_pd(x+j_coord_offsetA,x+j_coord_offsetB,
885 &jx0,&jy0,&jz0);
886
887 /* Calculate displacement vector */
888 dx00 = _mm_sub_pd(ix0,jx0);
889 dy00 = _mm_sub_pd(iy0,jy0);
890 dz00 = _mm_sub_pd(iz0,jz0);
891 dx10 = _mm_sub_pd(ix1,jx0);
892 dy10 = _mm_sub_pd(iy1,jy0);
893 dz10 = _mm_sub_pd(iz1,jz0);
894 dx20 = _mm_sub_pd(ix2,jx0);
895 dy20 = _mm_sub_pd(iy2,jy0);
896 dz20 = _mm_sub_pd(iz2,jz0);
897 dx30 = _mm_sub_pd(ix3,jx0);
898 dy30 = _mm_sub_pd(iy3,jy0);
899 dz30 = _mm_sub_pd(iz3,jz0);
900
901 /* Calculate squared distance and things based on it */
902 rsq00 = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
903 rsq10 = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
904 rsq20 = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
905 rsq30 = gmx_mm_calc_rsq_pd(dx30,dy30,dz30);
906
907 rinv00 = gmx_mm_invsqrt_pd(rsq00);
908 rinv10 = gmx_mm_invsqrt_pd(rsq10);
909 rinv20 = gmx_mm_invsqrt_pd(rsq20);
910 rinv30 = gmx_mm_invsqrt_pd(rsq30);
911
912 rinvsq10 = _mm_mul_pd(rinv10,rinv10);
913 rinvsq20 = _mm_mul_pd(rinv20,rinv20);
914 rinvsq30 = _mm_mul_pd(rinv30,rinv30);
915
916 /* Load parameters for j particles */
917 jq0 = gmx_mm_load_2real_swizzle_pd(charge+jnrA+0,charge+jnrB+0);
918 vdwjidx0A = 2*vdwtype[jnrA+0];
919 vdwjidx0B = 2*vdwtype[jnrB+0];
920
921 fjx0 = _mm_setzero_pd();
922 fjy0 = _mm_setzero_pd();
923 fjz0 = _mm_setzero_pd();
924
925 /**************************
926 * CALCULATE INTERACTIONS *
927 **************************/
928
929 r00 = _mm_mul_pd(rsq00,rinv00);
930
931 /* Compute parameters for interactions between i and j atoms */
932 gmx_mm_load_2pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,
933 vdwparam+vdwioffset0+vdwjidx0B,&c6_00,&c12_00);
934
935 /* Calculate table index by multiplying r with table scale and truncate to integer */
936 rt = _mm_mul_pd(r00,vftabscale);
937 vfitab = _mm_cvttpd_epi32(rt);
938 vfeps = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
939 vfitab = _mm_slli_epi32(vfitab,3);
940
941 /* CUBIC SPLINE TABLE DISPERSION */
942 Y = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) );
943 F = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,1) );
944 GMX_MM_TRANSPOSE2_PD(Y,F);
945 G = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) +2);
946 H = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,1) +2);
947 GMX_MM_TRANSPOSE2_PD(G,H);
948 Heps = _mm_mul_pd(vfeps,H);
949 Fp = _mm_add_pd(F,_mm_mul_pd(vfeps,_mm_add_pd(G,Heps)));
950 FF = _mm_add_pd(Fp,_mm_mul_pd(vfeps,_mm_add_pd(G,_mm_add_pd(Heps,Heps))));
951 fvdw6 = _mm_mul_pd(c6_00,FF);
952
953 /* CUBIC SPLINE TABLE REPULSION */
954 vfitab = _mm_add_epi32(vfitab,ifour);
955 Y = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) );
956 F = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,1) );
957 GMX_MM_TRANSPOSE2_PD(Y,F);
958 G = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) +2);
959 H = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,1) +2);
960 GMX_MM_TRANSPOSE2_PD(G,H);
961 Heps = _mm_mul_pd(vfeps,H);
962 Fp = _mm_add_pd(F,_mm_mul_pd(vfeps,_mm_add_pd(G,Heps)));
963 FF = _mm_add_pd(Fp,_mm_mul_pd(vfeps,_mm_add_pd(G,_mm_add_pd(Heps,Heps))));
964 fvdw12 = _mm_mul_pd(c12_00,FF);
965 fvdw = _mm_xor_pd(signbit,_mm_mul_pd(_mm_add_pd(fvdw6,fvdw12),_mm_mul_pd(vftabscale,rinv00)));
966
967 fscal = fvdw;
968
969 /* Calculate temporary vectorial force */
970 tx = _mm_mul_pd(fscal,dx00);
971 ty = _mm_mul_pd(fscal,dy00);
972 tz = _mm_mul_pd(fscal,dz00);
973
974 /* Update vectorial force */
975 fix0 = _mm_add_pd(fix0,tx);
976 fiy0 = _mm_add_pd(fiy0,ty);
977 fiz0 = _mm_add_pd(fiz0,tz);
978
979 fjx0 = _mm_add_pd(fjx0,tx);
980 fjy0 = _mm_add_pd(fjy0,ty);
981 fjz0 = _mm_add_pd(fjz0,tz);
982
983 /**************************
984 * CALCULATE INTERACTIONS *
985 **************************/
986
987 r10 = _mm_mul_pd(rsq10,rinv10);
988
989 /* Compute parameters for interactions between i and j atoms */
990 qq10 = _mm_mul_pd(iq1,jq0);
991
992 /* EWALD ELECTROSTATICS */
993
994 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
995 ewrt = _mm_mul_pd(r10,ewtabscale);
996 ewitab = _mm_cvttpd_epi32(ewrt);
997 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
998 gmx_mm_load_2pair_swizzle_pd(ewtab+gmx_mm_extract_epi32(ewitab,0),ewtab+gmx_mm_extract_epi32(ewitab,1),
999 &ewtabF,&ewtabFn);
1000 felec = _mm_add_pd(_mm_mul_pd( _mm_sub_pd(one,eweps),ewtabF),_mm_mul_pd(eweps,ewtabFn));
1001 felec = _mm_mul_pd(_mm_mul_pd(qq10,rinv10),_mm_sub_pd(rinvsq10,felec));
1002
1003 fscal = felec;
1004
1005 /* Calculate temporary vectorial force */
1006 tx = _mm_mul_pd(fscal,dx10);
1007 ty = _mm_mul_pd(fscal,dy10);
1008 tz = _mm_mul_pd(fscal,dz10);
1009
1010 /* Update vectorial force */
1011 fix1 = _mm_add_pd(fix1,tx);
1012 fiy1 = _mm_add_pd(fiy1,ty);
1013 fiz1 = _mm_add_pd(fiz1,tz);
1014
1015 fjx0 = _mm_add_pd(fjx0,tx);
1016 fjy0 = _mm_add_pd(fjy0,ty);
1017 fjz0 = _mm_add_pd(fjz0,tz);
1018
1019 /**************************
1020 * CALCULATE INTERACTIONS *
1021 **************************/
1022
1023 r20 = _mm_mul_pd(rsq20,rinv20);
1024
1025 /* Compute parameters for interactions between i and j atoms */
1026 qq20 = _mm_mul_pd(iq2,jq0);
1027
1028 /* EWALD ELECTROSTATICS */
1029
1030 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1031 ewrt = _mm_mul_pd(r20,ewtabscale);
1032 ewitab = _mm_cvttpd_epi32(ewrt);
1033 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1034 gmx_mm_load_2pair_swizzle_pd(ewtab+gmx_mm_extract_epi32(ewitab,0),ewtab+gmx_mm_extract_epi32(ewitab,1),
1035 &ewtabF,&ewtabFn);
1036 felec = _mm_add_pd(_mm_mul_pd( _mm_sub_pd(one,eweps),ewtabF),_mm_mul_pd(eweps,ewtabFn));
1037 felec = _mm_mul_pd(_mm_mul_pd(qq20,rinv20),_mm_sub_pd(rinvsq20,felec));
1038
1039 fscal = felec;
1040
1041 /* Calculate temporary vectorial force */
1042 tx = _mm_mul_pd(fscal,dx20);
1043 ty = _mm_mul_pd(fscal,dy20);
1044 tz = _mm_mul_pd(fscal,dz20);
1045
1046 /* Update vectorial force */
1047 fix2 = _mm_add_pd(fix2,tx);
1048 fiy2 = _mm_add_pd(fiy2,ty);
1049 fiz2 = _mm_add_pd(fiz2,tz);
1050
1051 fjx0 = _mm_add_pd(fjx0,tx);
1052 fjy0 = _mm_add_pd(fjy0,ty);
1053 fjz0 = _mm_add_pd(fjz0,tz);
1054
1055 /**************************
1056 * CALCULATE INTERACTIONS *
1057 **************************/
1058
1059 r30 = _mm_mul_pd(rsq30,rinv30);
1060
1061 /* Compute parameters for interactions between i and j atoms */
1062 qq30 = _mm_mul_pd(iq3,jq0);
1063
1064 /* EWALD ELECTROSTATICS */
1065
1066 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1067 ewrt = _mm_mul_pd(r30,ewtabscale);
1068 ewitab = _mm_cvttpd_epi32(ewrt);
1069 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1070 gmx_mm_load_2pair_swizzle_pd(ewtab+gmx_mm_extract_epi32(ewitab,0),ewtab+gmx_mm_extract_epi32(ewitab,1),
1071 &ewtabF,&ewtabFn);
1072 felec = _mm_add_pd(_mm_mul_pd( _mm_sub_pd(one,eweps),ewtabF),_mm_mul_pd(eweps,ewtabFn));
1073 felec = _mm_mul_pd(_mm_mul_pd(qq30,rinv30),_mm_sub_pd(rinvsq30,felec));
1074
1075 fscal = felec;
1076
1077 /* Calculate temporary vectorial force */
1078 tx = _mm_mul_pd(fscal,dx30);
1079 ty = _mm_mul_pd(fscal,dy30);
1080 tz = _mm_mul_pd(fscal,dz30);
1081
1082 /* Update vectorial force */
1083 fix3 = _mm_add_pd(fix3,tx);
1084 fiy3 = _mm_add_pd(fiy3,ty);
1085 fiz3 = _mm_add_pd(fiz3,tz);
1086
1087 fjx0 = _mm_add_pd(fjx0,tx);
1088 fjy0 = _mm_add_pd(fjy0,ty);
1089 fjz0 = _mm_add_pd(fjz0,tz);
1090
1091 gmx_mm_decrement_1rvec_2ptr_swizzle_pd(f+j_coord_offsetA,f+j_coord_offsetB,fjx0,fjy0,fjz0);
1092
1093 /* Inner loop uses 159 flops */
1094 }
1095
1096 if(jidx<j_index_end)
1097 {
1098
1099 jnrA = jjnr[jidx];
1100 j_coord_offsetA = DIM*jnrA;
1101
1102 /* load j atom coordinates */
1103 gmx_mm_load_1rvec_1ptr_swizzle_pd(x+j_coord_offsetA,
1104 &jx0,&jy0,&jz0);
1105
1106 /* Calculate displacement vector */
1107 dx00 = _mm_sub_pd(ix0,jx0);
1108 dy00 = _mm_sub_pd(iy0,jy0);
1109 dz00 = _mm_sub_pd(iz0,jz0);
1110 dx10 = _mm_sub_pd(ix1,jx0);
1111 dy10 = _mm_sub_pd(iy1,jy0);
1112 dz10 = _mm_sub_pd(iz1,jz0);
1113 dx20 = _mm_sub_pd(ix2,jx0);
1114 dy20 = _mm_sub_pd(iy2,jy0);
1115 dz20 = _mm_sub_pd(iz2,jz0);
1116 dx30 = _mm_sub_pd(ix3,jx0);
1117 dy30 = _mm_sub_pd(iy3,jy0);
1118 dz30 = _mm_sub_pd(iz3,jz0);
1119
1120 /* Calculate squared distance and things based on it */
1121 rsq00 = gmx_mm_calc_rsq_pd(dx00,dy00,dz00);
1122 rsq10 = gmx_mm_calc_rsq_pd(dx10,dy10,dz10);
1123 rsq20 = gmx_mm_calc_rsq_pd(dx20,dy20,dz20);
1124 rsq30 = gmx_mm_calc_rsq_pd(dx30,dy30,dz30);
1125
1126 rinv00 = gmx_mm_invsqrt_pd(rsq00);
1127 rinv10 = gmx_mm_invsqrt_pd(rsq10);
1128 rinv20 = gmx_mm_invsqrt_pd(rsq20);
1129 rinv30 = gmx_mm_invsqrt_pd(rsq30);
1130
1131 rinvsq10 = _mm_mul_pd(rinv10,rinv10);
1132 rinvsq20 = _mm_mul_pd(rinv20,rinv20);
1133 rinvsq30 = _mm_mul_pd(rinv30,rinv30);
1134
1135 /* Load parameters for j particles */
1136 jq0 = _mm_load_sd(charge+jnrA+0);
1137 vdwjidx0A = 2*vdwtype[jnrA+0];
1138
1139 fjx0 = _mm_setzero_pd();
1140 fjy0 = _mm_setzero_pd();
1141 fjz0 = _mm_setzero_pd();
1142
1143 /**************************
1144 * CALCULATE INTERACTIONS *
1145 **************************/
1146
1147 r00 = _mm_mul_pd(rsq00,rinv00);
1148
1149 /* Compute parameters for interactions between i and j atoms */
1150 gmx_mm_load_1pair_swizzle_pd(vdwparam+vdwioffset0+vdwjidx0A,&c6_00,&c12_00);
1151
1152 /* Calculate table index by multiplying r with table scale and truncate to integer */
1153 rt = _mm_mul_pd(r00,vftabscale);
1154 vfitab = _mm_cvttpd_epi32(rt);
1155 vfeps = _mm_sub_pd(rt,_mm_round_pd(rt, _MM_FROUND_FLOOR));
1156 vfitab = _mm_slli_epi32(vfitab,3);
1157
1158 /* CUBIC SPLINE TABLE DISPERSION */
1159 Y = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) );
1160 F = _mm_setzero_pd();
1161 GMX_MM_TRANSPOSE2_PD(Y,F);
1162 G = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) +2);
1163 H = _mm_setzero_pd();
1164 GMX_MM_TRANSPOSE2_PD(G,H);
1165 Heps = _mm_mul_pd(vfeps,H);
1166 Fp = _mm_add_pd(F,_mm_mul_pd(vfeps,_mm_add_pd(G,Heps)));
1167 FF = _mm_add_pd(Fp,_mm_mul_pd(vfeps,_mm_add_pd(G,_mm_add_pd(Heps,Heps))));
1168 fvdw6 = _mm_mul_pd(c6_00,FF);
1169
1170 /* CUBIC SPLINE TABLE REPULSION */
1171 vfitab = _mm_add_epi32(vfitab,ifour);
1172 Y = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) );
1173 F = _mm_setzero_pd();
1174 GMX_MM_TRANSPOSE2_PD(Y,F);
1175 G = _mm_load_pd( vftab + gmx_mm_extract_epi32(vfitab,0) +2);
1176 H = _mm_setzero_pd();
1177 GMX_MM_TRANSPOSE2_PD(G,H);
1178 Heps = _mm_mul_pd(vfeps,H);
1179 Fp = _mm_add_pd(F,_mm_mul_pd(vfeps,_mm_add_pd(G,Heps)));
1180 FF = _mm_add_pd(Fp,_mm_mul_pd(vfeps,_mm_add_pd(G,_mm_add_pd(Heps,Heps))));
1181 fvdw12 = _mm_mul_pd(c12_00,FF);
1182 fvdw = _mm_xor_pd(signbit,_mm_mul_pd(_mm_add_pd(fvdw6,fvdw12),_mm_mul_pd(vftabscale,rinv00)));
1183
1184 fscal = fvdw;
1185
1186 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1187
1188 /* Calculate temporary vectorial force */
1189 tx = _mm_mul_pd(fscal,dx00);
1190 ty = _mm_mul_pd(fscal,dy00);
1191 tz = _mm_mul_pd(fscal,dz00);
1192
1193 /* Update vectorial force */
1194 fix0 = _mm_add_pd(fix0,tx);
1195 fiy0 = _mm_add_pd(fiy0,ty);
1196 fiz0 = _mm_add_pd(fiz0,tz);
1197
1198 fjx0 = _mm_add_pd(fjx0,tx);
1199 fjy0 = _mm_add_pd(fjy0,ty);
1200 fjz0 = _mm_add_pd(fjz0,tz);
1201
1202 /**************************
1203 * CALCULATE INTERACTIONS *
1204 **************************/
1205
1206 r10 = _mm_mul_pd(rsq10,rinv10);
1207
1208 /* Compute parameters for interactions between i and j atoms */
1209 qq10 = _mm_mul_pd(iq1,jq0);
1210
1211 /* EWALD ELECTROSTATICS */
1212
1213 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1214 ewrt = _mm_mul_pd(r10,ewtabscale);
1215 ewitab = _mm_cvttpd_epi32(ewrt);
1216 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1217 gmx_mm_load_1pair_swizzle_pd(ewtab+gmx_mm_extract_epi32(ewitab,0),&ewtabF,&ewtabFn);
1218 felec = _mm_add_pd(_mm_mul_pd( _mm_sub_pd(one,eweps),ewtabF),_mm_mul_pd(eweps,ewtabFn));
1219 felec = _mm_mul_pd(_mm_mul_pd(qq10,rinv10),_mm_sub_pd(rinvsq10,felec));
1220
1221 fscal = felec;
1222
1223 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1224
1225 /* Calculate temporary vectorial force */
1226 tx = _mm_mul_pd(fscal,dx10);
1227 ty = _mm_mul_pd(fscal,dy10);
1228 tz = _mm_mul_pd(fscal,dz10);
1229
1230 /* Update vectorial force */
1231 fix1 = _mm_add_pd(fix1,tx);
1232 fiy1 = _mm_add_pd(fiy1,ty);
1233 fiz1 = _mm_add_pd(fiz1,tz);
1234
1235 fjx0 = _mm_add_pd(fjx0,tx);
1236 fjy0 = _mm_add_pd(fjy0,ty);
1237 fjz0 = _mm_add_pd(fjz0,tz);
1238
1239 /**************************
1240 * CALCULATE INTERACTIONS *
1241 **************************/
1242
1243 r20 = _mm_mul_pd(rsq20,rinv20);
1244
1245 /* Compute parameters for interactions between i and j atoms */
1246 qq20 = _mm_mul_pd(iq2,jq0);
1247
1248 /* EWALD ELECTROSTATICS */
1249
1250 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1251 ewrt = _mm_mul_pd(r20,ewtabscale);
1252 ewitab = _mm_cvttpd_epi32(ewrt);
1253 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1254 gmx_mm_load_1pair_swizzle_pd(ewtab+gmx_mm_extract_epi32(ewitab,0),&ewtabF,&ewtabFn);
1255 felec = _mm_add_pd(_mm_mul_pd( _mm_sub_pd(one,eweps),ewtabF),_mm_mul_pd(eweps,ewtabFn));
1256 felec = _mm_mul_pd(_mm_mul_pd(qq20,rinv20),_mm_sub_pd(rinvsq20,felec));
1257
1258 fscal = felec;
1259
1260 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1261
1262 /* Calculate temporary vectorial force */
1263 tx = _mm_mul_pd(fscal,dx20);
1264 ty = _mm_mul_pd(fscal,dy20);
1265 tz = _mm_mul_pd(fscal,dz20);
1266
1267 /* Update vectorial force */
1268 fix2 = _mm_add_pd(fix2,tx);
1269 fiy2 = _mm_add_pd(fiy2,ty);
1270 fiz2 = _mm_add_pd(fiz2,tz);
1271
1272 fjx0 = _mm_add_pd(fjx0,tx);
1273 fjy0 = _mm_add_pd(fjy0,ty);
1274 fjz0 = _mm_add_pd(fjz0,tz);
1275
1276 /**************************
1277 * CALCULATE INTERACTIONS *
1278 **************************/
1279
1280 r30 = _mm_mul_pd(rsq30,rinv30);
1281
1282 /* Compute parameters for interactions between i and j atoms */
1283 qq30 = _mm_mul_pd(iq3,jq0);
1284
1285 /* EWALD ELECTROSTATICS */
1286
1287 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1288 ewrt = _mm_mul_pd(r30,ewtabscale);
1289 ewitab = _mm_cvttpd_epi32(ewrt);
1290 eweps = _mm_sub_pd(ewrt,_mm_round_pd(ewrt, _MM_FROUND_FLOOR));
1291 gmx_mm_load_1pair_swizzle_pd(ewtab+gmx_mm_extract_epi32(ewitab,0),&ewtabF,&ewtabFn);
1292 felec = _mm_add_pd(_mm_mul_pd( _mm_sub_pd(one,eweps),ewtabF),_mm_mul_pd(eweps,ewtabFn));
1293 felec = _mm_mul_pd(_mm_mul_pd(qq30,rinv30),_mm_sub_pd(rinvsq30,felec));
1294
1295 fscal = felec;
1296
1297 fscal = _mm_unpacklo_pd(fscal,_mm_setzero_pd());
1298
1299 /* Calculate temporary vectorial force */
1300 tx = _mm_mul_pd(fscal,dx30);
1301 ty = _mm_mul_pd(fscal,dy30);
1302 tz = _mm_mul_pd(fscal,dz30);
1303
1304 /* Update vectorial force */
1305 fix3 = _mm_add_pd(fix3,tx);
1306 fiy3 = _mm_add_pd(fiy3,ty);
1307 fiz3 = _mm_add_pd(fiz3,tz);
1308
1309 fjx0 = _mm_add_pd(fjx0,tx);
1310 fjy0 = _mm_add_pd(fjy0,ty);
1311 fjz0 = _mm_add_pd(fjz0,tz);
1312
1313 gmx_mm_decrement_1rvec_1ptr_swizzle_pd(f+j_coord_offsetA,fjx0,fjy0,fjz0);
1314
1315 /* Inner loop uses 159 flops */
1316 }
1317
1318 /* End of innermost loop */
1319
1320 gmx_mm_update_iforce_4atom_swizzle_pd(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,fix3,fiy3,fiz3,
1321 f+i_coord_offset,fshift+i_shift_offset);
1322
1323 /* Increment number of inner iterations */
1324 inneriter += j_index_end - j_index_start;
1325
1326 /* Outer loop uses 24 flops */
1327 }
1328
1329 /* Increment number of outer iterations */
1330 outeriter += nri;
1331
1332 /* Update outer/inner flops */
1333
1334 inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4_F,outeriter*24 + inneriter*159);
1335 }
The ultimate molecular dynamics simulation package