c8d647153c0c79241443e8bf09c42b6f0645b3c4
[alexxy/gromacs.git] / src / gmxlib / nonbonded / nb_kernel_sse2_single / nb_kernel_ElecCoul_VdwLJ_GeomW3P1_sse2_single.c
1 /*
2  * Note: this file was generated by the Gromacs sse2_single kernel generator.
3  *
4  *                This source code is part of
5  *
6  *                 G   R   O   M   A   C   S
7  *
8  * Copyright (c) 2001-2012, The GROMACS Development Team
9  *
10  * Gromacs is a library for molecular simulation and trajectory analysis,
11  * written by Erik Lindahl, David van der Spoel, Berk Hess, and others - for
12  * a full list of developers and information, check out http://www.gromacs.org
13  *
14  * This program is free software; you can redistribute it and/or modify it under
15  * the terms of the GNU Lesser General Public License as published by the Free
16  * Software Foundation; either version 2 of the License, or (at your option) any
17  * later version.
18  *
19  * To help fund GROMACS development, we humbly ask that you cite
20  * the papers people have written on it - you can find them on the website.
21  */
22 #ifdef HAVE_CONFIG_H
23 #include <config.h>
24 #endif
25
26 #include <math.h>
27
28 #include "../nb_kernel.h"
29 #include "types/simple.h"
30 #include "vec.h"
31 #include "nrnb.h"
32
33 #include "gmx_math_x86_sse2_single.h"
34 #include "kernelutil_x86_sse2_single.h"
35
36 /*
37  * Gromacs nonbonded kernel:   nb_kernel_ElecCoul_VdwLJ_GeomW3P1_VF_sse2_single
38  * Electrostatics interaction: Coulomb
39  * VdW interaction:            LennardJones
40  * Geometry:                   Water3-Particle
41  * Calculate force/pot:        PotentialAndForce
42  */
43 void
44 nb_kernel_ElecCoul_VdwLJ_GeomW3P1_VF_sse2_single
45                     (t_nblist * gmx_restrict                nlist,
46                      rvec * gmx_restrict                    xx,
47                      rvec * gmx_restrict                    ff,
48                      t_forcerec * gmx_restrict              fr,
49                      t_mdatoms * gmx_restrict               mdatoms,
50                      nb_kernel_data_t * gmx_restrict        kernel_data,
51                      t_nrnb * gmx_restrict                  nrnb)
52 {
53     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or 
54      * just 0 for non-waters.
55      * Suffixes A,B,C,D refer to j loop unrolling done with SSE, e.g. for the four different
56      * jnr indices corresponding to data put in the four positions in the SIMD register.
57      */
58     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
59     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
60     int              jnrA,jnrB,jnrC,jnrD;
61     int              jnrlistA,jnrlistB,jnrlistC,jnrlistD;
62     int              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
63     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
64     real             rcutoff_scalar;
65     real             *shiftvec,*fshift,*x,*f;
66     real             *fjptrA,*fjptrB,*fjptrC,*fjptrD;
67     real             scratch[4*DIM];
68     __m128           tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
69     int              vdwioffset0;
70     __m128           ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
71     int              vdwioffset1;
72     __m128           ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
73     int              vdwioffset2;
74     __m128           ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
75     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
76     __m128           jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
77     __m128           dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
78     __m128           dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
79     __m128           dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
80     __m128           velec,felec,velecsum,facel,crf,krf,krf2;
81     real             *charge;
82     int              nvdwtype;
83     __m128           rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
84     int              *vdwtype;
85     real             *vdwparam;
86     __m128           one_sixth   = _mm_set1_ps(1.0/6.0);
87     __m128           one_twelfth = _mm_set1_ps(1.0/12.0);
88     __m128           dummy_mask,cutoff_mask;
89     __m128           signbit = _mm_castsi128_ps( _mm_set1_epi32(0x80000000) );
90     __m128           one     = _mm_set1_ps(1.0);
91     __m128           two     = _mm_set1_ps(2.0);
92     x                = xx[0];
93     f                = ff[0];
94
95     nri              = nlist->nri;
96     iinr             = nlist->iinr;
97     jindex           = nlist->jindex;
98     jjnr             = nlist->jjnr;
99     shiftidx         = nlist->shift;
100     gid              = nlist->gid;
101     shiftvec         = fr->shift_vec[0];
102     fshift           = fr->fshift[0];
103     facel            = _mm_set1_ps(fr->epsfac);
104     charge           = mdatoms->chargeA;
105     nvdwtype         = fr->ntype;
106     vdwparam         = fr->nbfp;
107     vdwtype          = mdatoms->typeA;
108
109     /* Setup water-specific parameters */
110     inr              = nlist->iinr[0];
111     iq0              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+0]));
112     iq1              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+1]));
113     iq2              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+2]));
114     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
115
116     /* Avoid stupid compiler warnings */
117     jnrA = jnrB = jnrC = jnrD = 0;
118     j_coord_offsetA = 0;
119     j_coord_offsetB = 0;
120     j_coord_offsetC = 0;
121     j_coord_offsetD = 0;
122
123     outeriter        = 0;
124     inneriter        = 0;
125
126     for(iidx=0;iidx<4*DIM;iidx++)
127     {
128         scratch[iidx] = 0.0;
129     }  
130
131     /* Start outer loop over neighborlists */
132     for(iidx=0; iidx<nri; iidx++)
133     {
134         /* Load shift vector for this list */
135         i_shift_offset   = DIM*shiftidx[iidx];
136
137         /* Load limits for loop over neighbors */
138         j_index_start    = jindex[iidx];
139         j_index_end      = jindex[iidx+1];
140
141         /* Get outer coordinate index */
142         inr              = iinr[iidx];
143         i_coord_offset   = DIM*inr;
144
145         /* Load i particle coords and add shift vector */
146         gmx_mm_load_shift_and_3rvec_broadcast_ps(shiftvec+i_shift_offset,x+i_coord_offset,
147                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2);
148         
149         fix0             = _mm_setzero_ps();
150         fiy0             = _mm_setzero_ps();
151         fiz0             = _mm_setzero_ps();
152         fix1             = _mm_setzero_ps();
153         fiy1             = _mm_setzero_ps();
154         fiz1             = _mm_setzero_ps();
155         fix2             = _mm_setzero_ps();
156         fiy2             = _mm_setzero_ps();
157         fiz2             = _mm_setzero_ps();
158
159         /* Reset potential sums */
160         velecsum         = _mm_setzero_ps();
161         vvdwsum          = _mm_setzero_ps();
162
163         /* Start inner kernel loop */
164         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
165         {
166
167             /* Get j neighbor index, and coordinate index */
168             jnrA             = jjnr[jidx];
169             jnrB             = jjnr[jidx+1];
170             jnrC             = jjnr[jidx+2];
171             jnrD             = jjnr[jidx+3];
172             j_coord_offsetA  = DIM*jnrA;
173             j_coord_offsetB  = DIM*jnrB;
174             j_coord_offsetC  = DIM*jnrC;
175             j_coord_offsetD  = DIM*jnrD;
176
177             /* load j atom coordinates */
178             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
179                                               x+j_coord_offsetC,x+j_coord_offsetD,
180                                               &jx0,&jy0,&jz0);
181
182             /* Calculate displacement vector */
183             dx00             = _mm_sub_ps(ix0,jx0);
184             dy00             = _mm_sub_ps(iy0,jy0);
185             dz00             = _mm_sub_ps(iz0,jz0);
186             dx10             = _mm_sub_ps(ix1,jx0);
187             dy10             = _mm_sub_ps(iy1,jy0);
188             dz10             = _mm_sub_ps(iz1,jz0);
189             dx20             = _mm_sub_ps(ix2,jx0);
190             dy20             = _mm_sub_ps(iy2,jy0);
191             dz20             = _mm_sub_ps(iz2,jz0);
192
193             /* Calculate squared distance and things based on it */
194             rsq00            = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
195             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
196             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
197
198             rinv00           = gmx_mm_invsqrt_ps(rsq00);
199             rinv10           = gmx_mm_invsqrt_ps(rsq10);
200             rinv20           = gmx_mm_invsqrt_ps(rsq20);
201
202             rinvsq00         = _mm_mul_ps(rinv00,rinv00);
203             rinvsq10         = _mm_mul_ps(rinv10,rinv10);
204             rinvsq20         = _mm_mul_ps(rinv20,rinv20);
205
206             /* Load parameters for j particles */
207             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
208                                                               charge+jnrC+0,charge+jnrD+0);
209             vdwjidx0A        = 2*vdwtype[jnrA+0];
210             vdwjidx0B        = 2*vdwtype[jnrB+0];
211             vdwjidx0C        = 2*vdwtype[jnrC+0];
212             vdwjidx0D        = 2*vdwtype[jnrD+0];
213
214             fjx0             = _mm_setzero_ps();
215             fjy0             = _mm_setzero_ps();
216             fjz0             = _mm_setzero_ps();
217
218             /**************************
219              * CALCULATE INTERACTIONS *
220              **************************/
221
222             /* Compute parameters for interactions between i and j atoms */
223             qq00             = _mm_mul_ps(iq0,jq0);
224             gmx_mm_load_4pair_swizzle_ps(vdwparam+vdwioffset0+vdwjidx0A,
225                                          vdwparam+vdwioffset0+vdwjidx0B,
226                                          vdwparam+vdwioffset0+vdwjidx0C,
227                                          vdwparam+vdwioffset0+vdwjidx0D,
228                                          &c6_00,&c12_00);
229
230             /* COULOMB ELECTROSTATICS */
231             velec            = _mm_mul_ps(qq00,rinv00);
232             felec            = _mm_mul_ps(velec,rinvsq00);
233
234             /* LENNARD-JONES DISPERSION/REPULSION */
235
236             rinvsix          = _mm_mul_ps(_mm_mul_ps(rinvsq00,rinvsq00),rinvsq00);
237             vvdw6            = _mm_mul_ps(c6_00,rinvsix);
238             vvdw12           = _mm_mul_ps(c12_00,_mm_mul_ps(rinvsix,rinvsix));
239             vvdw             = _mm_sub_ps( _mm_mul_ps(vvdw12,one_twelfth) , _mm_mul_ps(vvdw6,one_sixth) );
240             fvdw             = _mm_mul_ps(_mm_sub_ps(vvdw12,vvdw6),rinvsq00);
241
242             /* Update potential sum for this i atom from the interaction with this j atom. */
243             velecsum         = _mm_add_ps(velecsum,velec);
244             vvdwsum          = _mm_add_ps(vvdwsum,vvdw);
245
246             fscal            = _mm_add_ps(felec,fvdw);
247
248             /* Calculate temporary vectorial force */
249             tx               = _mm_mul_ps(fscal,dx00);
250             ty               = _mm_mul_ps(fscal,dy00);
251             tz               = _mm_mul_ps(fscal,dz00);
252
253             /* Update vectorial force */
254             fix0             = _mm_add_ps(fix0,tx);
255             fiy0             = _mm_add_ps(fiy0,ty);
256             fiz0             = _mm_add_ps(fiz0,tz);
257
258             fjx0             = _mm_add_ps(fjx0,tx);
259             fjy0             = _mm_add_ps(fjy0,ty);
260             fjz0             = _mm_add_ps(fjz0,tz);
261             
262             /**************************
263              * CALCULATE INTERACTIONS *
264              **************************/
265
266             /* Compute parameters for interactions between i and j atoms */
267             qq10             = _mm_mul_ps(iq1,jq0);
268
269             /* COULOMB ELECTROSTATICS */
270             velec            = _mm_mul_ps(qq10,rinv10);
271             felec            = _mm_mul_ps(velec,rinvsq10);
272
273             /* Update potential sum for this i atom from the interaction with this j atom. */
274             velecsum         = _mm_add_ps(velecsum,velec);
275
276             fscal            = felec;
277
278             /* Calculate temporary vectorial force */
279             tx               = _mm_mul_ps(fscal,dx10);
280             ty               = _mm_mul_ps(fscal,dy10);
281             tz               = _mm_mul_ps(fscal,dz10);
282
283             /* Update vectorial force */
284             fix1             = _mm_add_ps(fix1,tx);
285             fiy1             = _mm_add_ps(fiy1,ty);
286             fiz1             = _mm_add_ps(fiz1,tz);
287
288             fjx0             = _mm_add_ps(fjx0,tx);
289             fjy0             = _mm_add_ps(fjy0,ty);
290             fjz0             = _mm_add_ps(fjz0,tz);
291             
292             /**************************
293              * CALCULATE INTERACTIONS *
294              **************************/
295
296             /* Compute parameters for interactions between i and j atoms */
297             qq20             = _mm_mul_ps(iq2,jq0);
298
299             /* COULOMB ELECTROSTATICS */
300             velec            = _mm_mul_ps(qq20,rinv20);
301             felec            = _mm_mul_ps(velec,rinvsq20);
302
303             /* Update potential sum for this i atom from the interaction with this j atom. */
304             velecsum         = _mm_add_ps(velecsum,velec);
305
306             fscal            = felec;
307
308             /* Calculate temporary vectorial force */
309             tx               = _mm_mul_ps(fscal,dx20);
310             ty               = _mm_mul_ps(fscal,dy20);
311             tz               = _mm_mul_ps(fscal,dz20);
312
313             /* Update vectorial force */
314             fix2             = _mm_add_ps(fix2,tx);
315             fiy2             = _mm_add_ps(fiy2,ty);
316             fiz2             = _mm_add_ps(fiz2,tz);
317
318             fjx0             = _mm_add_ps(fjx0,tx);
319             fjy0             = _mm_add_ps(fjy0,ty);
320             fjz0             = _mm_add_ps(fjz0,tz);
321             
322             fjptrA             = f+j_coord_offsetA;
323             fjptrB             = f+j_coord_offsetB;
324             fjptrC             = f+j_coord_offsetC;
325             fjptrD             = f+j_coord_offsetD;
326
327             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,fjx0,fjy0,fjz0);
328
329             /* Inner loop uses 96 flops */
330         }
331
332         if(jidx<j_index_end)
333         {
334
335             /* Get j neighbor index, and coordinate index */
336             jnrlistA         = jjnr[jidx];
337             jnrlistB         = jjnr[jidx+1];
338             jnrlistC         = jjnr[jidx+2];
339             jnrlistD         = jjnr[jidx+3];
340             /* Sign of each element will be negative for non-real atoms.
341              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
342              * so use it as val = _mm_andnot_ps(mask,val) to clear dummy entries.
343              */
344             dummy_mask = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
345             jnrA       = (jnrlistA>=0) ? jnrlistA : 0;
346             jnrB       = (jnrlistB>=0) ? jnrlistB : 0;
347             jnrC       = (jnrlistC>=0) ? jnrlistC : 0;
348             jnrD       = (jnrlistD>=0) ? jnrlistD : 0;
349             j_coord_offsetA  = DIM*jnrA;
350             j_coord_offsetB  = DIM*jnrB;
351             j_coord_offsetC  = DIM*jnrC;
352             j_coord_offsetD  = DIM*jnrD;
353
354             /* load j atom coordinates */
355             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
356                                               x+j_coord_offsetC,x+j_coord_offsetD,
357                                               &jx0,&jy0,&jz0);
358
359             /* Calculate displacement vector */
360             dx00             = _mm_sub_ps(ix0,jx0);
361             dy00             = _mm_sub_ps(iy0,jy0);
362             dz00             = _mm_sub_ps(iz0,jz0);
363             dx10             = _mm_sub_ps(ix1,jx0);
364             dy10             = _mm_sub_ps(iy1,jy0);
365             dz10             = _mm_sub_ps(iz1,jz0);
366             dx20             = _mm_sub_ps(ix2,jx0);
367             dy20             = _mm_sub_ps(iy2,jy0);
368             dz20             = _mm_sub_ps(iz2,jz0);
369
370             /* Calculate squared distance and things based on it */
371             rsq00            = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
372             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
373             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
374
375             rinv00           = gmx_mm_invsqrt_ps(rsq00);
376             rinv10           = gmx_mm_invsqrt_ps(rsq10);
377             rinv20           = gmx_mm_invsqrt_ps(rsq20);
378
379             rinvsq00         = _mm_mul_ps(rinv00,rinv00);
380             rinvsq10         = _mm_mul_ps(rinv10,rinv10);
381             rinvsq20         = _mm_mul_ps(rinv20,rinv20);
382
383             /* Load parameters for j particles */
384             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
385                                                               charge+jnrC+0,charge+jnrD+0);
386             vdwjidx0A        = 2*vdwtype[jnrA+0];
387             vdwjidx0B        = 2*vdwtype[jnrB+0];
388             vdwjidx0C        = 2*vdwtype[jnrC+0];
389             vdwjidx0D        = 2*vdwtype[jnrD+0];
390
391             fjx0             = _mm_setzero_ps();
392             fjy0             = _mm_setzero_ps();
393             fjz0             = _mm_setzero_ps();
394
395             /**************************
396              * CALCULATE INTERACTIONS *
397              **************************/
398
399             /* Compute parameters for interactions between i and j atoms */
400             qq00             = _mm_mul_ps(iq0,jq0);
401             gmx_mm_load_4pair_swizzle_ps(vdwparam+vdwioffset0+vdwjidx0A,
402                                          vdwparam+vdwioffset0+vdwjidx0B,
403                                          vdwparam+vdwioffset0+vdwjidx0C,
404                                          vdwparam+vdwioffset0+vdwjidx0D,
405                                          &c6_00,&c12_00);
406
407             /* COULOMB ELECTROSTATICS */
408             velec            = _mm_mul_ps(qq00,rinv00);
409             felec            = _mm_mul_ps(velec,rinvsq00);
410
411             /* LENNARD-JONES DISPERSION/REPULSION */
412
413             rinvsix          = _mm_mul_ps(_mm_mul_ps(rinvsq00,rinvsq00),rinvsq00);
414             vvdw6            = _mm_mul_ps(c6_00,rinvsix);
415             vvdw12           = _mm_mul_ps(c12_00,_mm_mul_ps(rinvsix,rinvsix));
416             vvdw             = _mm_sub_ps( _mm_mul_ps(vvdw12,one_twelfth) , _mm_mul_ps(vvdw6,one_sixth) );
417             fvdw             = _mm_mul_ps(_mm_sub_ps(vvdw12,vvdw6),rinvsq00);
418
419             /* Update potential sum for this i atom from the interaction with this j atom. */
420             velec            = _mm_andnot_ps(dummy_mask,velec);
421             velecsum         = _mm_add_ps(velecsum,velec);
422             vvdw             = _mm_andnot_ps(dummy_mask,vvdw);
423             vvdwsum          = _mm_add_ps(vvdwsum,vvdw);
424
425             fscal            = _mm_add_ps(felec,fvdw);
426
427             fscal            = _mm_andnot_ps(dummy_mask,fscal);
428
429             /* Calculate temporary vectorial force */
430             tx               = _mm_mul_ps(fscal,dx00);
431             ty               = _mm_mul_ps(fscal,dy00);
432             tz               = _mm_mul_ps(fscal,dz00);
433
434             /* Update vectorial force */
435             fix0             = _mm_add_ps(fix0,tx);
436             fiy0             = _mm_add_ps(fiy0,ty);
437             fiz0             = _mm_add_ps(fiz0,tz);
438
439             fjx0             = _mm_add_ps(fjx0,tx);
440             fjy0             = _mm_add_ps(fjy0,ty);
441             fjz0             = _mm_add_ps(fjz0,tz);
442             
443             /**************************
444              * CALCULATE INTERACTIONS *
445              **************************/
446
447             /* Compute parameters for interactions between i and j atoms */
448             qq10             = _mm_mul_ps(iq1,jq0);
449
450             /* COULOMB ELECTROSTATICS */
451             velec            = _mm_mul_ps(qq10,rinv10);
452             felec            = _mm_mul_ps(velec,rinvsq10);
453
454             /* Update potential sum for this i atom from the interaction with this j atom. */
455             velec            = _mm_andnot_ps(dummy_mask,velec);
456             velecsum         = _mm_add_ps(velecsum,velec);
457
458             fscal            = felec;
459
460             fscal            = _mm_andnot_ps(dummy_mask,fscal);
461
462             /* Calculate temporary vectorial force */
463             tx               = _mm_mul_ps(fscal,dx10);
464             ty               = _mm_mul_ps(fscal,dy10);
465             tz               = _mm_mul_ps(fscal,dz10);
466
467             /* Update vectorial force */
468             fix1             = _mm_add_ps(fix1,tx);
469             fiy1             = _mm_add_ps(fiy1,ty);
470             fiz1             = _mm_add_ps(fiz1,tz);
471
472             fjx0             = _mm_add_ps(fjx0,tx);
473             fjy0             = _mm_add_ps(fjy0,ty);
474             fjz0             = _mm_add_ps(fjz0,tz);
475             
476             /**************************
477              * CALCULATE INTERACTIONS *
478              **************************/
479
480             /* Compute parameters for interactions between i and j atoms */
481             qq20             = _mm_mul_ps(iq2,jq0);
482
483             /* COULOMB ELECTROSTATICS */
484             velec            = _mm_mul_ps(qq20,rinv20);
485             felec            = _mm_mul_ps(velec,rinvsq20);
486
487             /* Update potential sum for this i atom from the interaction with this j atom. */
488             velec            = _mm_andnot_ps(dummy_mask,velec);
489             velecsum         = _mm_add_ps(velecsum,velec);
490
491             fscal            = felec;
492
493             fscal            = _mm_andnot_ps(dummy_mask,fscal);
494
495             /* Calculate temporary vectorial force */
496             tx               = _mm_mul_ps(fscal,dx20);
497             ty               = _mm_mul_ps(fscal,dy20);
498             tz               = _mm_mul_ps(fscal,dz20);
499
500             /* Update vectorial force */
501             fix2             = _mm_add_ps(fix2,tx);
502             fiy2             = _mm_add_ps(fiy2,ty);
503             fiz2             = _mm_add_ps(fiz2,tz);
504
505             fjx0             = _mm_add_ps(fjx0,tx);
506             fjy0             = _mm_add_ps(fjy0,ty);
507             fjz0             = _mm_add_ps(fjz0,tz);
508             
509             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
510             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
511             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
512             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
513
514             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,fjx0,fjy0,fjz0);
515
516             /* Inner loop uses 96 flops */
517         }
518
519         /* End of innermost loop */
520
521         gmx_mm_update_iforce_3atom_swizzle_ps(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,
522                                               f+i_coord_offset,fshift+i_shift_offset);
523
524         ggid                        = gid[iidx];
525         /* Update potential energies */
526         gmx_mm_update_1pot_ps(velecsum,kernel_data->energygrp_elec+ggid);
527         gmx_mm_update_1pot_ps(vvdwsum,kernel_data->energygrp_vdw+ggid);
528
529         /* Increment number of inner iterations */
530         inneriter                  += j_index_end - j_index_start;
531
532         /* Outer loop uses 20 flops */
533     }
534
535     /* Increment number of outer iterations */
536     outeriter        += nri;
537
538     /* Update outer/inner flops */
539
540     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W3_VF,outeriter*20 + inneriter*96);
541 }
542 /*
543  * Gromacs nonbonded kernel:   nb_kernel_ElecCoul_VdwLJ_GeomW3P1_F_sse2_single
544  * Electrostatics interaction: Coulomb
545  * VdW interaction:            LennardJones
546  * Geometry:                   Water3-Particle
547  * Calculate force/pot:        Force
548  */
549 void
550 nb_kernel_ElecCoul_VdwLJ_GeomW3P1_F_sse2_single
551                     (t_nblist * gmx_restrict                nlist,
552                      rvec * gmx_restrict                    xx,
553                      rvec * gmx_restrict                    ff,
554                      t_forcerec * gmx_restrict              fr,
555                      t_mdatoms * gmx_restrict               mdatoms,
556                      nb_kernel_data_t * gmx_restrict        kernel_data,
557                      t_nrnb * gmx_restrict                  nrnb)
558 {
559     /* Suffixes 0,1,2,3 refer to particle indices for waters in the inner or outer loop, or 
560      * just 0 for non-waters.
561      * Suffixes A,B,C,D refer to j loop unrolling done with SSE, e.g. for the four different
562      * jnr indices corresponding to data put in the four positions in the SIMD register.
563      */
564     int              i_shift_offset,i_coord_offset,outeriter,inneriter;
565     int              j_index_start,j_index_end,jidx,nri,inr,ggid,iidx;
566     int              jnrA,jnrB,jnrC,jnrD;
567     int              jnrlistA,jnrlistB,jnrlistC,jnrlistD;
568     int              j_coord_offsetA,j_coord_offsetB,j_coord_offsetC,j_coord_offsetD;
569     int              *iinr,*jindex,*jjnr,*shiftidx,*gid;
570     real             rcutoff_scalar;
571     real             *shiftvec,*fshift,*x,*f;
572     real             *fjptrA,*fjptrB,*fjptrC,*fjptrD;
573     real             scratch[4*DIM];
574     __m128           tx,ty,tz,fscal,rcutoff,rcutoff2,jidxall;
575     int              vdwioffset0;
576     __m128           ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
577     int              vdwioffset1;
578     __m128           ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
579     int              vdwioffset2;
580     __m128           ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
581     int              vdwjidx0A,vdwjidx0B,vdwjidx0C,vdwjidx0D;
582     __m128           jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
583     __m128           dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00;
584     __m128           dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10;
585     __m128           dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20;
586     __m128           velec,felec,velecsum,facel,crf,krf,krf2;
587     real             *charge;
588     int              nvdwtype;
589     __m128           rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,sh_vdw_invrcut6;
590     int              *vdwtype;
591     real             *vdwparam;
592     __m128           one_sixth   = _mm_set1_ps(1.0/6.0);
593     __m128           one_twelfth = _mm_set1_ps(1.0/12.0);
594     __m128           dummy_mask,cutoff_mask;
595     __m128           signbit = _mm_castsi128_ps( _mm_set1_epi32(0x80000000) );
596     __m128           one     = _mm_set1_ps(1.0);
597     __m128           two     = _mm_set1_ps(2.0);
598     x                = xx[0];
599     f                = ff[0];
600
601     nri              = nlist->nri;
602     iinr             = nlist->iinr;
603     jindex           = nlist->jindex;
604     jjnr             = nlist->jjnr;
605     shiftidx         = nlist->shift;
606     gid              = nlist->gid;
607     shiftvec         = fr->shift_vec[0];
608     fshift           = fr->fshift[0];
609     facel            = _mm_set1_ps(fr->epsfac);
610     charge           = mdatoms->chargeA;
611     nvdwtype         = fr->ntype;
612     vdwparam         = fr->nbfp;
613     vdwtype          = mdatoms->typeA;
614
615     /* Setup water-specific parameters */
616     inr              = nlist->iinr[0];
617     iq0              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+0]));
618     iq1              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+1]));
619     iq2              = _mm_mul_ps(facel,_mm_set1_ps(charge[inr+2]));
620     vdwioffset0      = 2*nvdwtype*vdwtype[inr+0];
621
622     /* Avoid stupid compiler warnings */
623     jnrA = jnrB = jnrC = jnrD = 0;
624     j_coord_offsetA = 0;
625     j_coord_offsetB = 0;
626     j_coord_offsetC = 0;
627     j_coord_offsetD = 0;
628
629     outeriter        = 0;
630     inneriter        = 0;
631
632     for(iidx=0;iidx<4*DIM;iidx++)
633     {
634         scratch[iidx] = 0.0;
635     }  
636
637     /* Start outer loop over neighborlists */
638     for(iidx=0; iidx<nri; iidx++)
639     {
640         /* Load shift vector for this list */
641         i_shift_offset   = DIM*shiftidx[iidx];
642
643         /* Load limits for loop over neighbors */
644         j_index_start    = jindex[iidx];
645         j_index_end      = jindex[iidx+1];
646
647         /* Get outer coordinate index */
648         inr              = iinr[iidx];
649         i_coord_offset   = DIM*inr;
650
651         /* Load i particle coords and add shift vector */
652         gmx_mm_load_shift_and_3rvec_broadcast_ps(shiftvec+i_shift_offset,x+i_coord_offset,
653                                                  &ix0,&iy0,&iz0,&ix1,&iy1,&iz1,&ix2,&iy2,&iz2);
654         
655         fix0             = _mm_setzero_ps();
656         fiy0             = _mm_setzero_ps();
657         fiz0             = _mm_setzero_ps();
658         fix1             = _mm_setzero_ps();
659         fiy1             = _mm_setzero_ps();
660         fiz1             = _mm_setzero_ps();
661         fix2             = _mm_setzero_ps();
662         fiy2             = _mm_setzero_ps();
663         fiz2             = _mm_setzero_ps();
664
665         /* Start inner kernel loop */
666         for(jidx=j_index_start; jidx<j_index_end && jjnr[jidx+3]>=0; jidx+=4)
667         {
668
669             /* Get j neighbor index, and coordinate index */
670             jnrA             = jjnr[jidx];
671             jnrB             = jjnr[jidx+1];
672             jnrC             = jjnr[jidx+2];
673             jnrD             = jjnr[jidx+3];
674             j_coord_offsetA  = DIM*jnrA;
675             j_coord_offsetB  = DIM*jnrB;
676             j_coord_offsetC  = DIM*jnrC;
677             j_coord_offsetD  = DIM*jnrD;
678
679             /* load j atom coordinates */
680             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
681                                               x+j_coord_offsetC,x+j_coord_offsetD,
682                                               &jx0,&jy0,&jz0);
683
684             /* Calculate displacement vector */
685             dx00             = _mm_sub_ps(ix0,jx0);
686             dy00             = _mm_sub_ps(iy0,jy0);
687             dz00             = _mm_sub_ps(iz0,jz0);
688             dx10             = _mm_sub_ps(ix1,jx0);
689             dy10             = _mm_sub_ps(iy1,jy0);
690             dz10             = _mm_sub_ps(iz1,jz0);
691             dx20             = _mm_sub_ps(ix2,jx0);
692             dy20             = _mm_sub_ps(iy2,jy0);
693             dz20             = _mm_sub_ps(iz2,jz0);
694
695             /* Calculate squared distance and things based on it */
696             rsq00            = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
697             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
698             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
699
700             rinv00           = gmx_mm_invsqrt_ps(rsq00);
701             rinv10           = gmx_mm_invsqrt_ps(rsq10);
702             rinv20           = gmx_mm_invsqrt_ps(rsq20);
703
704             rinvsq00         = _mm_mul_ps(rinv00,rinv00);
705             rinvsq10         = _mm_mul_ps(rinv10,rinv10);
706             rinvsq20         = _mm_mul_ps(rinv20,rinv20);
707
708             /* Load parameters for j particles */
709             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
710                                                               charge+jnrC+0,charge+jnrD+0);
711             vdwjidx0A        = 2*vdwtype[jnrA+0];
712             vdwjidx0B        = 2*vdwtype[jnrB+0];
713             vdwjidx0C        = 2*vdwtype[jnrC+0];
714             vdwjidx0D        = 2*vdwtype[jnrD+0];
715
716             fjx0             = _mm_setzero_ps();
717             fjy0             = _mm_setzero_ps();
718             fjz0             = _mm_setzero_ps();
719
720             /**************************
721              * CALCULATE INTERACTIONS *
722              **************************/
723
724             /* Compute parameters for interactions between i and j atoms */
725             qq00             = _mm_mul_ps(iq0,jq0);
726             gmx_mm_load_4pair_swizzle_ps(vdwparam+vdwioffset0+vdwjidx0A,
727                                          vdwparam+vdwioffset0+vdwjidx0B,
728                                          vdwparam+vdwioffset0+vdwjidx0C,
729                                          vdwparam+vdwioffset0+vdwjidx0D,
730                                          &c6_00,&c12_00);
731
732             /* COULOMB ELECTROSTATICS */
733             velec            = _mm_mul_ps(qq00,rinv00);
734             felec            = _mm_mul_ps(velec,rinvsq00);
735
736             /* LENNARD-JONES DISPERSION/REPULSION */
737
738             rinvsix          = _mm_mul_ps(_mm_mul_ps(rinvsq00,rinvsq00),rinvsq00);
739             fvdw             = _mm_mul_ps(_mm_sub_ps(_mm_mul_ps(c12_00,rinvsix),c6_00),_mm_mul_ps(rinvsix,rinvsq00));
740
741             fscal            = _mm_add_ps(felec,fvdw);
742
743             /* Calculate temporary vectorial force */
744             tx               = _mm_mul_ps(fscal,dx00);
745             ty               = _mm_mul_ps(fscal,dy00);
746             tz               = _mm_mul_ps(fscal,dz00);
747
748             /* Update vectorial force */
749             fix0             = _mm_add_ps(fix0,tx);
750             fiy0             = _mm_add_ps(fiy0,ty);
751             fiz0             = _mm_add_ps(fiz0,tz);
752
753             fjx0             = _mm_add_ps(fjx0,tx);
754             fjy0             = _mm_add_ps(fjy0,ty);
755             fjz0             = _mm_add_ps(fjz0,tz);
756             
757             /**************************
758              * CALCULATE INTERACTIONS *
759              **************************/
760
761             /* Compute parameters for interactions between i and j atoms */
762             qq10             = _mm_mul_ps(iq1,jq0);
763
764             /* COULOMB ELECTROSTATICS */
765             velec            = _mm_mul_ps(qq10,rinv10);
766             felec            = _mm_mul_ps(velec,rinvsq10);
767
768             fscal            = felec;
769
770             /* Calculate temporary vectorial force */
771             tx               = _mm_mul_ps(fscal,dx10);
772             ty               = _mm_mul_ps(fscal,dy10);
773             tz               = _mm_mul_ps(fscal,dz10);
774
775             /* Update vectorial force */
776             fix1             = _mm_add_ps(fix1,tx);
777             fiy1             = _mm_add_ps(fiy1,ty);
778             fiz1             = _mm_add_ps(fiz1,tz);
779
780             fjx0             = _mm_add_ps(fjx0,tx);
781             fjy0             = _mm_add_ps(fjy0,ty);
782             fjz0             = _mm_add_ps(fjz0,tz);
783             
784             /**************************
785              * CALCULATE INTERACTIONS *
786              **************************/
787
788             /* Compute parameters for interactions between i and j atoms */
789             qq20             = _mm_mul_ps(iq2,jq0);
790
791             /* COULOMB ELECTROSTATICS */
792             velec            = _mm_mul_ps(qq20,rinv20);
793             felec            = _mm_mul_ps(velec,rinvsq20);
794
795             fscal            = felec;
796
797             /* Calculate temporary vectorial force */
798             tx               = _mm_mul_ps(fscal,dx20);
799             ty               = _mm_mul_ps(fscal,dy20);
800             tz               = _mm_mul_ps(fscal,dz20);
801
802             /* Update vectorial force */
803             fix2             = _mm_add_ps(fix2,tx);
804             fiy2             = _mm_add_ps(fiy2,ty);
805             fiz2             = _mm_add_ps(fiz2,tz);
806
807             fjx0             = _mm_add_ps(fjx0,tx);
808             fjy0             = _mm_add_ps(fjy0,ty);
809             fjz0             = _mm_add_ps(fjz0,tz);
810             
811             fjptrA             = f+j_coord_offsetA;
812             fjptrB             = f+j_coord_offsetB;
813             fjptrC             = f+j_coord_offsetC;
814             fjptrD             = f+j_coord_offsetD;
815
816             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,fjx0,fjy0,fjz0);
817
818             /* Inner loop uses 88 flops */
819         }
820
821         if(jidx<j_index_end)
822         {
823
824             /* Get j neighbor index, and coordinate index */
825             jnrlistA         = jjnr[jidx];
826             jnrlistB         = jjnr[jidx+1];
827             jnrlistC         = jjnr[jidx+2];
828             jnrlistD         = jjnr[jidx+3];
829             /* Sign of each element will be negative for non-real atoms.
830              * This mask will be 0xFFFFFFFF for dummy entries and 0x0 for real ones,
831              * so use it as val = _mm_andnot_ps(mask,val) to clear dummy entries.
832              */
833             dummy_mask = gmx_mm_castsi128_ps(_mm_cmplt_epi32(_mm_loadu_si128((const __m128i *)(jjnr+jidx)),_mm_setzero_si128()));
834             jnrA       = (jnrlistA>=0) ? jnrlistA : 0;
835             jnrB       = (jnrlistB>=0) ? jnrlistB : 0;
836             jnrC       = (jnrlistC>=0) ? jnrlistC : 0;
837             jnrD       = (jnrlistD>=0) ? jnrlistD : 0;
838             j_coord_offsetA  = DIM*jnrA;
839             j_coord_offsetB  = DIM*jnrB;
840             j_coord_offsetC  = DIM*jnrC;
841             j_coord_offsetD  = DIM*jnrD;
842
843             /* load j atom coordinates */
844             gmx_mm_load_1rvec_4ptr_swizzle_ps(x+j_coord_offsetA,x+j_coord_offsetB,
845                                               x+j_coord_offsetC,x+j_coord_offsetD,
846                                               &jx0,&jy0,&jz0);
847
848             /* Calculate displacement vector */
849             dx00             = _mm_sub_ps(ix0,jx0);
850             dy00             = _mm_sub_ps(iy0,jy0);
851             dz00             = _mm_sub_ps(iz0,jz0);
852             dx10             = _mm_sub_ps(ix1,jx0);
853             dy10             = _mm_sub_ps(iy1,jy0);
854             dz10             = _mm_sub_ps(iz1,jz0);
855             dx20             = _mm_sub_ps(ix2,jx0);
856             dy20             = _mm_sub_ps(iy2,jy0);
857             dz20             = _mm_sub_ps(iz2,jz0);
858
859             /* Calculate squared distance and things based on it */
860             rsq00            = gmx_mm_calc_rsq_ps(dx00,dy00,dz00);
861             rsq10            = gmx_mm_calc_rsq_ps(dx10,dy10,dz10);
862             rsq20            = gmx_mm_calc_rsq_ps(dx20,dy20,dz20);
863
864             rinv00           = gmx_mm_invsqrt_ps(rsq00);
865             rinv10           = gmx_mm_invsqrt_ps(rsq10);
866             rinv20           = gmx_mm_invsqrt_ps(rsq20);
867
868             rinvsq00         = _mm_mul_ps(rinv00,rinv00);
869             rinvsq10         = _mm_mul_ps(rinv10,rinv10);
870             rinvsq20         = _mm_mul_ps(rinv20,rinv20);
871
872             /* Load parameters for j particles */
873             jq0              = gmx_mm_load_4real_swizzle_ps(charge+jnrA+0,charge+jnrB+0,
874                                                               charge+jnrC+0,charge+jnrD+0);
875             vdwjidx0A        = 2*vdwtype[jnrA+0];
876             vdwjidx0B        = 2*vdwtype[jnrB+0];
877             vdwjidx0C        = 2*vdwtype[jnrC+0];
878             vdwjidx0D        = 2*vdwtype[jnrD+0];
879
880             fjx0             = _mm_setzero_ps();
881             fjy0             = _mm_setzero_ps();
882             fjz0             = _mm_setzero_ps();
883
884             /**************************
885              * CALCULATE INTERACTIONS *
886              **************************/
887
888             /* Compute parameters for interactions between i and j atoms */
889             qq00             = _mm_mul_ps(iq0,jq0);
890             gmx_mm_load_4pair_swizzle_ps(vdwparam+vdwioffset0+vdwjidx0A,
891                                          vdwparam+vdwioffset0+vdwjidx0B,
892                                          vdwparam+vdwioffset0+vdwjidx0C,
893                                          vdwparam+vdwioffset0+vdwjidx0D,
894                                          &c6_00,&c12_00);
895
896             /* COULOMB ELECTROSTATICS */
897             velec            = _mm_mul_ps(qq00,rinv00);
898             felec            = _mm_mul_ps(velec,rinvsq00);
899
900             /* LENNARD-JONES DISPERSION/REPULSION */
901
902             rinvsix          = _mm_mul_ps(_mm_mul_ps(rinvsq00,rinvsq00),rinvsq00);
903             fvdw             = _mm_mul_ps(_mm_sub_ps(_mm_mul_ps(c12_00,rinvsix),c6_00),_mm_mul_ps(rinvsix,rinvsq00));
904
905             fscal            = _mm_add_ps(felec,fvdw);
906
907             fscal            = _mm_andnot_ps(dummy_mask,fscal);
908
909             /* Calculate temporary vectorial force */
910             tx               = _mm_mul_ps(fscal,dx00);
911             ty               = _mm_mul_ps(fscal,dy00);
912             tz               = _mm_mul_ps(fscal,dz00);
913
914             /* Update vectorial force */
915             fix0             = _mm_add_ps(fix0,tx);
916             fiy0             = _mm_add_ps(fiy0,ty);
917             fiz0             = _mm_add_ps(fiz0,tz);
918
919             fjx0             = _mm_add_ps(fjx0,tx);
920             fjy0             = _mm_add_ps(fjy0,ty);
921             fjz0             = _mm_add_ps(fjz0,tz);
922             
923             /**************************
924              * CALCULATE INTERACTIONS *
925              **************************/
926
927             /* Compute parameters for interactions between i and j atoms */
928             qq10             = _mm_mul_ps(iq1,jq0);
929
930             /* COULOMB ELECTROSTATICS */
931             velec            = _mm_mul_ps(qq10,rinv10);
932             felec            = _mm_mul_ps(velec,rinvsq10);
933
934             fscal            = felec;
935
936             fscal            = _mm_andnot_ps(dummy_mask,fscal);
937
938             /* Calculate temporary vectorial force */
939             tx               = _mm_mul_ps(fscal,dx10);
940             ty               = _mm_mul_ps(fscal,dy10);
941             tz               = _mm_mul_ps(fscal,dz10);
942
943             /* Update vectorial force */
944             fix1             = _mm_add_ps(fix1,tx);
945             fiy1             = _mm_add_ps(fiy1,ty);
946             fiz1             = _mm_add_ps(fiz1,tz);
947
948             fjx0             = _mm_add_ps(fjx0,tx);
949             fjy0             = _mm_add_ps(fjy0,ty);
950             fjz0             = _mm_add_ps(fjz0,tz);
951             
952             /**************************
953              * CALCULATE INTERACTIONS *
954              **************************/
955
956             /* Compute parameters for interactions between i and j atoms */
957             qq20             = _mm_mul_ps(iq2,jq0);
958
959             /* COULOMB ELECTROSTATICS */
960             velec            = _mm_mul_ps(qq20,rinv20);
961             felec            = _mm_mul_ps(velec,rinvsq20);
962
963             fscal            = felec;
964
965             fscal            = _mm_andnot_ps(dummy_mask,fscal);
966
967             /* Calculate temporary vectorial force */
968             tx               = _mm_mul_ps(fscal,dx20);
969             ty               = _mm_mul_ps(fscal,dy20);
970             tz               = _mm_mul_ps(fscal,dz20);
971
972             /* Update vectorial force */
973             fix2             = _mm_add_ps(fix2,tx);
974             fiy2             = _mm_add_ps(fiy2,ty);
975             fiz2             = _mm_add_ps(fiz2,tz);
976
977             fjx0             = _mm_add_ps(fjx0,tx);
978             fjy0             = _mm_add_ps(fjy0,ty);
979             fjz0             = _mm_add_ps(fjz0,tz);
980             
981             fjptrA             = (jnrlistA>=0) ? f+j_coord_offsetA : scratch;
982             fjptrB             = (jnrlistB>=0) ? f+j_coord_offsetB : scratch;
983             fjptrC             = (jnrlistC>=0) ? f+j_coord_offsetC : scratch;
984             fjptrD             = (jnrlistD>=0) ? f+j_coord_offsetD : scratch;
985
986             gmx_mm_decrement_1rvec_4ptr_swizzle_ps(fjptrA,fjptrB,fjptrC,fjptrD,fjx0,fjy0,fjz0);
987
988             /* Inner loop uses 88 flops */
989         }
990
991         /* End of innermost loop */
992
993         gmx_mm_update_iforce_3atom_swizzle_ps(fix0,fiy0,fiz0,fix1,fiy1,fiz1,fix2,fiy2,fiz2,
994                                               f+i_coord_offset,fshift+i_shift_offset);
995
996         /* Increment number of inner iterations */
997         inneriter                  += j_index_end - j_index_start;
998
999         /* Outer loop uses 18 flops */
1000     }
1001
1002     /* Increment number of outer iterations */
1003     outeriter        += nri;
1004
1005     /* Update outer/inner flops */
1006
1007     inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W3_F,outeriter*18 + inneriter*88);
1008 }