?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Remove nb-parameters from t_forcerec
[gromacs.git] / src / gromacs / gmxlib / nonbonded / nb_kernel_c / nb_kernel_ElecRFCut_VdwLJSw_GeomW4W4_c.c
blobbfdd9fbf2506b825eb602f2eca6a0860e8f77fa1
1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 2012,2013,2014.2015,2017, 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 c kernel generator.
37 */
38 #include "gmxpre.h"
39
40 #include "config.h"
41
42 #include <math.h>
43
44 #include "../nb_kernel.h"
45 #include "gromacs/gmxlib/nrnb.h"
46
47 /*
48 * Gromacs nonbonded kernel: nb_kernel_ElecRFCut_VdwLJSw_GeomW4W4_VF_c
49 * Electrostatics interaction: ReactionField
50 * VdW interaction: LennardJones
51 * Geometry: Water4-Water4
52 * Calculate force/pot: PotentialAndForce
53 */
54 void
55 nb_kernel_ElecRFCut_VdwLJSw_GeomW4W4_VF_c
56 (t_nblist * gmx_restrict nlist,
57 rvec * gmx_restrict xx,
58 rvec * gmx_restrict ff,
59 struct t_forcerec * gmx_restrict fr,
60 t_mdatoms * gmx_restrict mdatoms,
61 nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
62 t_nrnb * gmx_restrict nrnb)
63 {
64 int i_shift_offset,i_coord_offset,j_coord_offset;
65 int j_index_start,j_index_end;
66 int nri,inr,ggid,iidx,jidx,jnr,outeriter,inneriter;
67 real shX,shY,shZ,tx,ty,tz,fscal,rcutoff,rcutoff2;
68 int *iinr,*jindex,*jjnr,*shiftidx,*gid;
69 real *shiftvec,*fshift,*x,*f;
70 int vdwioffset0;
71 real ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
72 int vdwioffset1;
73 real ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
74 int vdwioffset2;
75 real ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
76 int vdwioffset3;
77 real ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
78 int vdwjidx0;
79 real jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
80 int vdwjidx1;
81 real jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
82 int vdwjidx2;
83 real jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
84 int vdwjidx3;
85 real jx3,jy3,jz3,fjx3,fjy3,fjz3,jq3,isaj3;
86 real dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00,cexp1_00,cexp2_00;
87 real dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11,cexp1_11,cexp2_11;
88 real dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12,cexp1_12,cexp2_12;
89 real dx13,dy13,dz13,rsq13,rinv13,rinvsq13,r13,qq13,c6_13,c12_13,cexp1_13,cexp2_13;
90 real dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21,cexp1_21,cexp2_21;
91 real dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22,cexp1_22,cexp2_22;
92 real dx23,dy23,dz23,rsq23,rinv23,rinvsq23,r23,qq23,c6_23,c12_23,cexp1_23,cexp2_23;
93 real dx31,dy31,dz31,rsq31,rinv31,rinvsq31,r31,qq31,c6_31,c12_31,cexp1_31,cexp2_31;
94 real dx32,dy32,dz32,rsq32,rinv32,rinvsq32,r32,qq32,c6_32,c12_32,cexp1_32,cexp2_32;
95 real dx33,dy33,dz33,rsq33,rinv33,rinvsq33,r33,qq33,c6_33,c12_33,cexp1_33,cexp2_33;
96 real velec,felec,velecsum,facel,crf,krf,krf2;
97 real *charge;
98 int nvdwtype;
99 real rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,br,vvdwexp,sh_vdw_invrcut6;
100 int *vdwtype;
101 real *vdwparam;
102 real rswitch,swV3,swV4,swV5,swF2,swF3,swF4,d,d2,sw,dsw;
103
104 x = xx[0];
105 f = ff[0];
106
107 nri = nlist->nri;
108 iinr = nlist->iinr;
109 jindex = nlist->jindex;
110 jjnr = nlist->jjnr;
111 shiftidx = nlist->shift;
112 gid = nlist->gid;
113 shiftvec = fr->shift_vec[0];
114 fshift = fr->fshift[0];
115 facel = fr->ic->epsfac;
116 charge = mdatoms->chargeA;
117 krf = fr->ic->k_rf;
118 krf2 = krf*2.0;
119 crf = fr->ic->c_rf;
120 nvdwtype = fr->ntype;
121 vdwparam = fr->nbfp;
122 vdwtype = mdatoms->typeA;
123
124 /* Setup water-specific parameters */
125 inr = nlist->iinr[0];
126 iq1 = facel*charge[inr+1];
127 iq2 = facel*charge[inr+2];
128 iq3 = facel*charge[inr+3];
129 vdwioffset0 = 2*nvdwtype*vdwtype[inr+0];
130
131 jq1 = charge[inr+1];
132 jq2 = charge[inr+2];
133 jq3 = charge[inr+3];
134 vdwjidx0 = 2*vdwtype[inr+0];
135 c6_00 = vdwparam[vdwioffset0+vdwjidx0];
136 c12_00 = vdwparam[vdwioffset0+vdwjidx0+1];
137 qq11 = iq1*jq1;
138 qq12 = iq1*jq2;
139 qq13 = iq1*jq3;
140 qq21 = iq2*jq1;
141 qq22 = iq2*jq2;
142 qq23 = iq2*jq3;
143 qq31 = iq3*jq1;
144 qq32 = iq3*jq2;
145 qq33 = iq3*jq3;
146
147 /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
148 rcutoff = fr->ic->rcoulomb;
149 rcutoff2 = rcutoff*rcutoff;
150
151 rswitch = fr->ic->rvdw_switch;
152 /* Setup switch parameters */
153 d = rcutoff-rswitch;
154 swV3 = -10.0/(d*d*d);
155 swV4 = 15.0/(d*d*d*d);
156 swV5 = -6.0/(d*d*d*d*d);
157 swF2 = -30.0/(d*d*d);
158 swF3 = 60.0/(d*d*d*d);
159 swF4 = -30.0/(d*d*d*d*d);
160
161 outeriter = 0;
162 inneriter = 0;
163
164 /* Start outer loop over neighborlists */
165 for(iidx=0; iidx<nri; iidx++)
166 {
167 /* Load shift vector for this list */
168 i_shift_offset = DIM*shiftidx[iidx];
169 shX = shiftvec[i_shift_offset+XX];
170 shY = shiftvec[i_shift_offset+YY];
171 shZ = shiftvec[i_shift_offset+ZZ];
172
173 /* Load limits for loop over neighbors */
174 j_index_start = jindex[iidx];
175 j_index_end = jindex[iidx+1];
176
177 /* Get outer coordinate index */
178 inr = iinr[iidx];
179 i_coord_offset = DIM*inr;
180
181 /* Load i particle coords and add shift vector */
182 ix0 = shX + x[i_coord_offset+DIM*0+XX];
183 iy0 = shY + x[i_coord_offset+DIM*0+YY];
184 iz0 = shZ + x[i_coord_offset+DIM*0+ZZ];
185 ix1 = shX + x[i_coord_offset+DIM*1+XX];
186 iy1 = shY + x[i_coord_offset+DIM*1+YY];
187 iz1 = shZ + x[i_coord_offset+DIM*1+ZZ];
188 ix2 = shX + x[i_coord_offset+DIM*2+XX];
189 iy2 = shY + x[i_coord_offset+DIM*2+YY];
190 iz2 = shZ + x[i_coord_offset+DIM*2+ZZ];
191 ix3 = shX + x[i_coord_offset+DIM*3+XX];
192 iy3 = shY + x[i_coord_offset+DIM*3+YY];
193 iz3 = shZ + x[i_coord_offset+DIM*3+ZZ];
194
195 fix0 = 0.0;
196 fiy0 = 0.0;
197 fiz0 = 0.0;
198 fix1 = 0.0;
199 fiy1 = 0.0;
200 fiz1 = 0.0;
201 fix2 = 0.0;
202 fiy2 = 0.0;
203 fiz2 = 0.0;
204 fix3 = 0.0;
205 fiy3 = 0.0;
206 fiz3 = 0.0;
207
208 /* Reset potential sums */
209 velecsum = 0.0;
210 vvdwsum = 0.0;
211
212 /* Start inner kernel loop */
213 for(jidx=j_index_start; jidx<j_index_end; jidx++)
214 {
215 /* Get j neighbor index, and coordinate index */
216 jnr = jjnr[jidx];
217 j_coord_offset = DIM*jnr;
218
219 /* load j atom coordinates */
220 jx0 = x[j_coord_offset+DIM*0+XX];
221 jy0 = x[j_coord_offset+DIM*0+YY];
222 jz0 = x[j_coord_offset+DIM*0+ZZ];
223 jx1 = x[j_coord_offset+DIM*1+XX];
224 jy1 = x[j_coord_offset+DIM*1+YY];
225 jz1 = x[j_coord_offset+DIM*1+ZZ];
226 jx2 = x[j_coord_offset+DIM*2+XX];
227 jy2 = x[j_coord_offset+DIM*2+YY];
228 jz2 = x[j_coord_offset+DIM*2+ZZ];
229 jx3 = x[j_coord_offset+DIM*3+XX];
230 jy3 = x[j_coord_offset+DIM*3+YY];
231 jz3 = x[j_coord_offset+DIM*3+ZZ];
232
233 /* Calculate displacement vector */
234 dx00 = ix0 - jx0;
235 dy00 = iy0 - jy0;
236 dz00 = iz0 - jz0;
237 dx11 = ix1 - jx1;
238 dy11 = iy1 - jy1;
239 dz11 = iz1 - jz1;
240 dx12 = ix1 - jx2;
241 dy12 = iy1 - jy2;
242 dz12 = iz1 - jz2;
243 dx13 = ix1 - jx3;
244 dy13 = iy1 - jy3;
245 dz13 = iz1 - jz3;
246 dx21 = ix2 - jx1;
247 dy21 = iy2 - jy1;
248 dz21 = iz2 - jz1;
249 dx22 = ix2 - jx2;
250 dy22 = iy2 - jy2;
251 dz22 = iz2 - jz2;
252 dx23 = ix2 - jx3;
253 dy23 = iy2 - jy3;
254 dz23 = iz2 - jz3;
255 dx31 = ix3 - jx1;
256 dy31 = iy3 - jy1;
257 dz31 = iz3 - jz1;
258 dx32 = ix3 - jx2;
259 dy32 = iy3 - jy2;
260 dz32 = iz3 - jz2;
261 dx33 = ix3 - jx3;
262 dy33 = iy3 - jy3;
263 dz33 = iz3 - jz3;
264
265 /* Calculate squared distance and things based on it */
266 rsq00 = dx00*dx00+dy00*dy00+dz00*dz00;
267 rsq11 = dx11*dx11+dy11*dy11+dz11*dz11;
268 rsq12 = dx12*dx12+dy12*dy12+dz12*dz12;
269 rsq13 = dx13*dx13+dy13*dy13+dz13*dz13;
270 rsq21 = dx21*dx21+dy21*dy21+dz21*dz21;
271 rsq22 = dx22*dx22+dy22*dy22+dz22*dz22;
272 rsq23 = dx23*dx23+dy23*dy23+dz23*dz23;
273 rsq31 = dx31*dx31+dy31*dy31+dz31*dz31;
274 rsq32 = dx32*dx32+dy32*dy32+dz32*dz32;
275 rsq33 = dx33*dx33+dy33*dy33+dz33*dz33;
276
277 rinv00 = 1.0/sqrt(rsq00);
278 rinv11 = 1.0/sqrt(rsq11);
279 rinv12 = 1.0/sqrt(rsq12);
280 rinv13 = 1.0/sqrt(rsq13);
281 rinv21 = 1.0/sqrt(rsq21);
282 rinv22 = 1.0/sqrt(rsq22);
283 rinv23 = 1.0/sqrt(rsq23);
284 rinv31 = 1.0/sqrt(rsq31);
285 rinv32 = 1.0/sqrt(rsq32);
286 rinv33 = 1.0/sqrt(rsq33);
287
288 rinvsq00 = rinv00*rinv00;
289 rinvsq11 = rinv11*rinv11;
290 rinvsq12 = rinv12*rinv12;
291 rinvsq13 = rinv13*rinv13;
292 rinvsq21 = rinv21*rinv21;
293 rinvsq22 = rinv22*rinv22;
294 rinvsq23 = rinv23*rinv23;
295 rinvsq31 = rinv31*rinv31;
296 rinvsq32 = rinv32*rinv32;
297 rinvsq33 = rinv33*rinv33;
298
299 /**************************
300 * CALCULATE INTERACTIONS *
301 **************************/
302
303 if (rsq00<rcutoff2)
304 {
305
306 r00 = rsq00*rinv00;
307
308 /* LENNARD-JONES DISPERSION/REPULSION */
309
310 rinvsix = rinvsq00*rinvsq00*rinvsq00;
311 vvdw6 = c6_00*rinvsix;
312 vvdw12 = c12_00*rinvsix*rinvsix;
313 vvdw = vvdw12*(1.0/12.0) - vvdw6*(1.0/6.0);
314 fvdw = (vvdw12-vvdw6)*rinvsq00;
315
316 d = r00-rswitch;
317 d = (d>0.0) ? d : 0.0;
318 d2 = d*d;
319 sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
320
321 dsw = d2*(swF2+d*(swF3+d*swF4));
322
323 /* Evaluate switch function */
324 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
325 fvdw = fvdw*sw - rinv00*vvdw*dsw;
326 vvdw *= sw;
327
328 /* Update potential sums from outer loop */
329 vvdwsum += vvdw;
330
331 fscal = fvdw;
332
333 /* Calculate temporary vectorial force */
334 tx = fscal*dx00;
335 ty = fscal*dy00;
336 tz = fscal*dz00;
337
338 /* Update vectorial force */
339 fix0 += tx;
340 fiy0 += ty;
341 fiz0 += tz;
342 f[j_coord_offset+DIM*0+XX] -= tx;
343 f[j_coord_offset+DIM*0+YY] -= ty;
344 f[j_coord_offset+DIM*0+ZZ] -= tz;
345
346 }
347
348 /**************************
349 * CALCULATE INTERACTIONS *
350 **************************/
351
352 if (rsq11<rcutoff2)
353 {
354
355 /* REACTION-FIELD ELECTROSTATICS */
356 velec = qq11*(rinv11+krf*rsq11-crf);
357 felec = qq11*(rinv11*rinvsq11-krf2);
358
359 /* Update potential sums from outer loop */
360 velecsum += velec;
361
362 fscal = felec;
363
364 /* Calculate temporary vectorial force */
365 tx = fscal*dx11;
366 ty = fscal*dy11;
367 tz = fscal*dz11;
368
369 /* Update vectorial force */
370 fix1 += tx;
371 fiy1 += ty;
372 fiz1 += tz;
373 f[j_coord_offset+DIM*1+XX] -= tx;
374 f[j_coord_offset+DIM*1+YY] -= ty;
375 f[j_coord_offset+DIM*1+ZZ] -= tz;
376
377 }
378
379 /**************************
380 * CALCULATE INTERACTIONS *
381 **************************/
382
383 if (rsq12<rcutoff2)
384 {
385
386 /* REACTION-FIELD ELECTROSTATICS */
387 velec = qq12*(rinv12+krf*rsq12-crf);
388 felec = qq12*(rinv12*rinvsq12-krf2);
389
390 /* Update potential sums from outer loop */
391 velecsum += velec;
392
393 fscal = felec;
394
395 /* Calculate temporary vectorial force */
396 tx = fscal*dx12;
397 ty = fscal*dy12;
398 tz = fscal*dz12;
399
400 /* Update vectorial force */
401 fix1 += tx;
402 fiy1 += ty;
403 fiz1 += tz;
404 f[j_coord_offset+DIM*2+XX] -= tx;
405 f[j_coord_offset+DIM*2+YY] -= ty;
406 f[j_coord_offset+DIM*2+ZZ] -= tz;
407
408 }
409
410 /**************************
411 * CALCULATE INTERACTIONS *
412 **************************/
413
414 if (rsq13<rcutoff2)
415 {
416
417 /* REACTION-FIELD ELECTROSTATICS */
418 velec = qq13*(rinv13+krf*rsq13-crf);
419 felec = qq13*(rinv13*rinvsq13-krf2);
420
421 /* Update potential sums from outer loop */
422 velecsum += velec;
423
424 fscal = felec;
425
426 /* Calculate temporary vectorial force */
427 tx = fscal*dx13;
428 ty = fscal*dy13;
429 tz = fscal*dz13;
430
431 /* Update vectorial force */
432 fix1 += tx;
433 fiy1 += ty;
434 fiz1 += tz;
435 f[j_coord_offset+DIM*3+XX] -= tx;
436 f[j_coord_offset+DIM*3+YY] -= ty;
437 f[j_coord_offset+DIM*3+ZZ] -= tz;
438
439 }
440
441 /**************************
442 * CALCULATE INTERACTIONS *
443 **************************/
444
445 if (rsq21<rcutoff2)
446 {
447
448 /* REACTION-FIELD ELECTROSTATICS */
449 velec = qq21*(rinv21+krf*rsq21-crf);
450 felec = qq21*(rinv21*rinvsq21-krf2);
451
452 /* Update potential sums from outer loop */
453 velecsum += velec;
454
455 fscal = felec;
456
457 /* Calculate temporary vectorial force */
458 tx = fscal*dx21;
459 ty = fscal*dy21;
460 tz = fscal*dz21;
461
462 /* Update vectorial force */
463 fix2 += tx;
464 fiy2 += ty;
465 fiz2 += tz;
466 f[j_coord_offset+DIM*1+XX] -= tx;
467 f[j_coord_offset+DIM*1+YY] -= ty;
468 f[j_coord_offset+DIM*1+ZZ] -= tz;
469
470 }
471
472 /**************************
473 * CALCULATE INTERACTIONS *
474 **************************/
475
476 if (rsq22<rcutoff2)
477 {
478
479 /* REACTION-FIELD ELECTROSTATICS */
480 velec = qq22*(rinv22+krf*rsq22-crf);
481 felec = qq22*(rinv22*rinvsq22-krf2);
482
483 /* Update potential sums from outer loop */
484 velecsum += velec;
485
486 fscal = felec;
487
488 /* Calculate temporary vectorial force */
489 tx = fscal*dx22;
490 ty = fscal*dy22;
491 tz = fscal*dz22;
492
493 /* Update vectorial force */
494 fix2 += tx;
495 fiy2 += ty;
496 fiz2 += tz;
497 f[j_coord_offset+DIM*2+XX] -= tx;
498 f[j_coord_offset+DIM*2+YY] -= ty;
499 f[j_coord_offset+DIM*2+ZZ] -= tz;
500
501 }
502
503 /**************************
504 * CALCULATE INTERACTIONS *
505 **************************/
506
507 if (rsq23<rcutoff2)
508 {
509
510 /* REACTION-FIELD ELECTROSTATICS */
511 velec = qq23*(rinv23+krf*rsq23-crf);
512 felec = qq23*(rinv23*rinvsq23-krf2);
513
514 /* Update potential sums from outer loop */
515 velecsum += velec;
516
517 fscal = felec;
518
519 /* Calculate temporary vectorial force */
520 tx = fscal*dx23;
521 ty = fscal*dy23;
522 tz = fscal*dz23;
523
524 /* Update vectorial force */
525 fix2 += tx;
526 fiy2 += ty;
527 fiz2 += tz;
528 f[j_coord_offset+DIM*3+XX] -= tx;
529 f[j_coord_offset+DIM*3+YY] -= ty;
530 f[j_coord_offset+DIM*3+ZZ] -= tz;
531
532 }
533
534 /**************************
535 * CALCULATE INTERACTIONS *
536 **************************/
537
538 if (rsq31<rcutoff2)
539 {
540
541 /* REACTION-FIELD ELECTROSTATICS */
542 velec = qq31*(rinv31+krf*rsq31-crf);
543 felec = qq31*(rinv31*rinvsq31-krf2);
544
545 /* Update potential sums from outer loop */
546 velecsum += velec;
547
548 fscal = felec;
549
550 /* Calculate temporary vectorial force */
551 tx = fscal*dx31;
552 ty = fscal*dy31;
553 tz = fscal*dz31;
554
555 /* Update vectorial force */
556 fix3 += tx;
557 fiy3 += ty;
558 fiz3 += tz;
559 f[j_coord_offset+DIM*1+XX] -= tx;
560 f[j_coord_offset+DIM*1+YY] -= ty;
561 f[j_coord_offset+DIM*1+ZZ] -= tz;
562
563 }
564
565 /**************************
566 * CALCULATE INTERACTIONS *
567 **************************/
568
569 if (rsq32<rcutoff2)
570 {
571
572 /* REACTION-FIELD ELECTROSTATICS */
573 velec = qq32*(rinv32+krf*rsq32-crf);
574 felec = qq32*(rinv32*rinvsq32-krf2);
575
576 /* Update potential sums from outer loop */
577 velecsum += velec;
578
579 fscal = felec;
580
581 /* Calculate temporary vectorial force */
582 tx = fscal*dx32;
583 ty = fscal*dy32;
584 tz = fscal*dz32;
585
586 /* Update vectorial force */
587 fix3 += tx;
588 fiy3 += ty;
589 fiz3 += tz;
590 f[j_coord_offset+DIM*2+XX] -= tx;
591 f[j_coord_offset+DIM*2+YY] -= ty;
592 f[j_coord_offset+DIM*2+ZZ] -= tz;
593
594 }
595
596 /**************************
597 * CALCULATE INTERACTIONS *
598 **************************/
599
600 if (rsq33<rcutoff2)
601 {
602
603 /* REACTION-FIELD ELECTROSTATICS */
604 velec = qq33*(rinv33+krf*rsq33-crf);
605 felec = qq33*(rinv33*rinvsq33-krf2);
606
607 /* Update potential sums from outer loop */
608 velecsum += velec;
609
610 fscal = felec;
611
612 /* Calculate temporary vectorial force */
613 tx = fscal*dx33;
614 ty = fscal*dy33;
615 tz = fscal*dz33;
616
617 /* Update vectorial force */
618 fix3 += tx;
619 fiy3 += ty;
620 fiz3 += tz;
621 f[j_coord_offset+DIM*3+XX] -= tx;
622 f[j_coord_offset+DIM*3+YY] -= ty;
623 f[j_coord_offset+DIM*3+ZZ] -= tz;
624
625 }
626
627 /* Inner loop uses 332 flops */
628 }
629 /* End of innermost loop */
630
631 tx = ty = tz = 0;
632 f[i_coord_offset+DIM*0+XX] += fix0;
633 f[i_coord_offset+DIM*0+YY] += fiy0;
634 f[i_coord_offset+DIM*0+ZZ] += fiz0;
635 tx += fix0;
636 ty += fiy0;
637 tz += fiz0;
638 f[i_coord_offset+DIM*1+XX] += fix1;
639 f[i_coord_offset+DIM*1+YY] += fiy1;
640 f[i_coord_offset+DIM*1+ZZ] += fiz1;
641 tx += fix1;
642 ty += fiy1;
643 tz += fiz1;
644 f[i_coord_offset+DIM*2+XX] += fix2;
645 f[i_coord_offset+DIM*2+YY] += fiy2;
646 f[i_coord_offset+DIM*2+ZZ] += fiz2;
647 tx += fix2;
648 ty += fiy2;
649 tz += fiz2;
650 f[i_coord_offset+DIM*3+XX] += fix3;
651 f[i_coord_offset+DIM*3+YY] += fiy3;
652 f[i_coord_offset+DIM*3+ZZ] += fiz3;
653 tx += fix3;
654 ty += fiy3;
655 tz += fiz3;
656 fshift[i_shift_offset+XX] += tx;
657 fshift[i_shift_offset+YY] += ty;
658 fshift[i_shift_offset+ZZ] += tz;
659
660 ggid = gid[iidx];
661 /* Update potential energies */
662 kernel_data->energygrp_elec[ggid] += velecsum;
663 kernel_data->energygrp_vdw[ggid] += vvdwsum;
664
665 /* Increment number of inner iterations */
666 inneriter += j_index_end - j_index_start;
667
668 /* Outer loop uses 41 flops */
669 }
670
671 /* Increment number of outer iterations */
672 outeriter += nri;
673
674 /* Update outer/inner flops */
675
676 inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4W4_VF,outeriter*41 + inneriter*332);
677 }
678 /*
679 * Gromacs nonbonded kernel: nb_kernel_ElecRFCut_VdwLJSw_GeomW4W4_F_c
680 * Electrostatics interaction: ReactionField
681 * VdW interaction: LennardJones
682 * Geometry: Water4-Water4
683 * Calculate force/pot: Force
684 */
685 void
686 nb_kernel_ElecRFCut_VdwLJSw_GeomW4W4_F_c
687 (t_nblist * gmx_restrict nlist,
688 rvec * gmx_restrict xx,
689 rvec * gmx_restrict ff,
690 struct t_forcerec * gmx_restrict fr,
691 t_mdatoms * gmx_restrict mdatoms,
692 nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
693 t_nrnb * gmx_restrict nrnb)
694 {
695 int i_shift_offset,i_coord_offset,j_coord_offset;
696 int j_index_start,j_index_end;
697 int nri,inr,ggid,iidx,jidx,jnr,outeriter,inneriter;
698 real shX,shY,shZ,tx,ty,tz,fscal,rcutoff,rcutoff2;
699 int *iinr,*jindex,*jjnr,*shiftidx,*gid;
700 real *shiftvec,*fshift,*x,*f;
701 int vdwioffset0;
702 real ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
703 int vdwioffset1;
704 real ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
705 int vdwioffset2;
706 real ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
707 int vdwioffset3;
708 real ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
709 int vdwjidx0;
710 real jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
711 int vdwjidx1;
712 real jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
713 int vdwjidx2;
714 real jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
715 int vdwjidx3;
716 real jx3,jy3,jz3,fjx3,fjy3,fjz3,jq3,isaj3;
717 real dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00,cexp1_00,cexp2_00;
718 real dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11,cexp1_11,cexp2_11;
719 real dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12,cexp1_12,cexp2_12;
720 real dx13,dy13,dz13,rsq13,rinv13,rinvsq13,r13,qq13,c6_13,c12_13,cexp1_13,cexp2_13;
721 real dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21,cexp1_21,cexp2_21;
722 real dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22,cexp1_22,cexp2_22;
723 real dx23,dy23,dz23,rsq23,rinv23,rinvsq23,r23,qq23,c6_23,c12_23,cexp1_23,cexp2_23;
724 real dx31,dy31,dz31,rsq31,rinv31,rinvsq31,r31,qq31,c6_31,c12_31,cexp1_31,cexp2_31;
725 real dx32,dy32,dz32,rsq32,rinv32,rinvsq32,r32,qq32,c6_32,c12_32,cexp1_32,cexp2_32;
726 real dx33,dy33,dz33,rsq33,rinv33,rinvsq33,r33,qq33,c6_33,c12_33,cexp1_33,cexp2_33;
727 real velec,felec,velecsum,facel,crf,krf,krf2;
728 real *charge;
729 int nvdwtype;
730 real rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,br,vvdwexp,sh_vdw_invrcut6;
731 int *vdwtype;
732 real *vdwparam;
733 real rswitch,swV3,swV4,swV5,swF2,swF3,swF4,d,d2,sw,dsw;
734
735 x = xx[0];
736 f = ff[0];
737
738 nri = nlist->nri;
739 iinr = nlist->iinr;
740 jindex = nlist->jindex;
741 jjnr = nlist->jjnr;
742 shiftidx = nlist->shift;
743 gid = nlist->gid;
744 shiftvec = fr->shift_vec[0];
745 fshift = fr->fshift[0];
746 facel = fr->ic->epsfac;
747 charge = mdatoms->chargeA;
748 krf = fr->ic->k_rf;
749 krf2 = krf*2.0;
750 crf = fr->ic->c_rf;
751 nvdwtype = fr->ntype;
752 vdwparam = fr->nbfp;
753 vdwtype = mdatoms->typeA;
754
755 /* Setup water-specific parameters */
756 inr = nlist->iinr[0];
757 iq1 = facel*charge[inr+1];
758 iq2 = facel*charge[inr+2];
759 iq3 = facel*charge[inr+3];
760 vdwioffset0 = 2*nvdwtype*vdwtype[inr+0];
761
762 jq1 = charge[inr+1];
763 jq2 = charge[inr+2];
764 jq3 = charge[inr+3];
765 vdwjidx0 = 2*vdwtype[inr+0];
766 c6_00 = vdwparam[vdwioffset0+vdwjidx0];
767 c12_00 = vdwparam[vdwioffset0+vdwjidx0+1];
768 qq11 = iq1*jq1;
769 qq12 = iq1*jq2;
770 qq13 = iq1*jq3;
771 qq21 = iq2*jq1;
772 qq22 = iq2*jq2;
773 qq23 = iq2*jq3;
774 qq31 = iq3*jq1;
775 qq32 = iq3*jq2;
776 qq33 = iq3*jq3;
777
778 /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
779 rcutoff = fr->ic->rcoulomb;
780 rcutoff2 = rcutoff*rcutoff;
781
782 rswitch = fr->ic->rvdw_switch;
783 /* Setup switch parameters */
784 d = rcutoff-rswitch;
785 swV3 = -10.0/(d*d*d);
786 swV4 = 15.0/(d*d*d*d);
787 swV5 = -6.0/(d*d*d*d*d);
788 swF2 = -30.0/(d*d*d);
789 swF3 = 60.0/(d*d*d*d);
790 swF4 = -30.0/(d*d*d*d*d);
791
792 outeriter = 0;
793 inneriter = 0;
794
795 /* Start outer loop over neighborlists */
796 for(iidx=0; iidx<nri; iidx++)
797 {
798 /* Load shift vector for this list */
799 i_shift_offset = DIM*shiftidx[iidx];
800 shX = shiftvec[i_shift_offset+XX];
801 shY = shiftvec[i_shift_offset+YY];
802 shZ = shiftvec[i_shift_offset+ZZ];
803
804 /* Load limits for loop over neighbors */
805 j_index_start = jindex[iidx];
806 j_index_end = jindex[iidx+1];
807
808 /* Get outer coordinate index */
809 inr = iinr[iidx];
810 i_coord_offset = DIM*inr;
811
812 /* Load i particle coords and add shift vector */
813 ix0 = shX + x[i_coord_offset+DIM*0+XX];
814 iy0 = shY + x[i_coord_offset+DIM*0+YY];
815 iz0 = shZ + x[i_coord_offset+DIM*0+ZZ];
816 ix1 = shX + x[i_coord_offset+DIM*1+XX];
817 iy1 = shY + x[i_coord_offset+DIM*1+YY];
818 iz1 = shZ + x[i_coord_offset+DIM*1+ZZ];
819 ix2 = shX + x[i_coord_offset+DIM*2+XX];
820 iy2 = shY + x[i_coord_offset+DIM*2+YY];
821 iz2 = shZ + x[i_coord_offset+DIM*2+ZZ];
822 ix3 = shX + x[i_coord_offset+DIM*3+XX];
823 iy3 = shY + x[i_coord_offset+DIM*3+YY];
824 iz3 = shZ + x[i_coord_offset+DIM*3+ZZ];
825
826 fix0 = 0.0;
827 fiy0 = 0.0;
828 fiz0 = 0.0;
829 fix1 = 0.0;
830 fiy1 = 0.0;
831 fiz1 = 0.0;
832 fix2 = 0.0;
833 fiy2 = 0.0;
834 fiz2 = 0.0;
835 fix3 = 0.0;
836 fiy3 = 0.0;
837 fiz3 = 0.0;
838
839 /* Start inner kernel loop */
840 for(jidx=j_index_start; jidx<j_index_end; jidx++)
841 {
842 /* Get j neighbor index, and coordinate index */
843 jnr = jjnr[jidx];
844 j_coord_offset = DIM*jnr;
845
846 /* load j atom coordinates */
847 jx0 = x[j_coord_offset+DIM*0+XX];
848 jy0 = x[j_coord_offset+DIM*0+YY];
849 jz0 = x[j_coord_offset+DIM*0+ZZ];
850 jx1 = x[j_coord_offset+DIM*1+XX];
851 jy1 = x[j_coord_offset+DIM*1+YY];
852 jz1 = x[j_coord_offset+DIM*1+ZZ];
853 jx2 = x[j_coord_offset+DIM*2+XX];
854 jy2 = x[j_coord_offset+DIM*2+YY];
855 jz2 = x[j_coord_offset+DIM*2+ZZ];
856 jx3 = x[j_coord_offset+DIM*3+XX];
857 jy3 = x[j_coord_offset+DIM*3+YY];
858 jz3 = x[j_coord_offset+DIM*3+ZZ];
859
860 /* Calculate displacement vector */
861 dx00 = ix0 - jx0;
862 dy00 = iy0 - jy0;
863 dz00 = iz0 - jz0;
864 dx11 = ix1 - jx1;
865 dy11 = iy1 - jy1;
866 dz11 = iz1 - jz1;
867 dx12 = ix1 - jx2;
868 dy12 = iy1 - jy2;
869 dz12 = iz1 - jz2;
870 dx13 = ix1 - jx3;
871 dy13 = iy1 - jy3;
872 dz13 = iz1 - jz3;
873 dx21 = ix2 - jx1;
874 dy21 = iy2 - jy1;
875 dz21 = iz2 - jz1;
876 dx22 = ix2 - jx2;
877 dy22 = iy2 - jy2;
878 dz22 = iz2 - jz2;
879 dx23 = ix2 - jx3;
880 dy23 = iy2 - jy3;
881 dz23 = iz2 - jz3;
882 dx31 = ix3 - jx1;
883 dy31 = iy3 - jy1;
884 dz31 = iz3 - jz1;
885 dx32 = ix3 - jx2;
886 dy32 = iy3 - jy2;
887 dz32 = iz3 - jz2;
888 dx33 = ix3 - jx3;
889 dy33 = iy3 - jy3;
890 dz33 = iz3 - jz3;
891
892 /* Calculate squared distance and things based on it */
893 rsq00 = dx00*dx00+dy00*dy00+dz00*dz00;
894 rsq11 = dx11*dx11+dy11*dy11+dz11*dz11;
895 rsq12 = dx12*dx12+dy12*dy12+dz12*dz12;
896 rsq13 = dx13*dx13+dy13*dy13+dz13*dz13;
897 rsq21 = dx21*dx21+dy21*dy21+dz21*dz21;
898 rsq22 = dx22*dx22+dy22*dy22+dz22*dz22;
899 rsq23 = dx23*dx23+dy23*dy23+dz23*dz23;
900 rsq31 = dx31*dx31+dy31*dy31+dz31*dz31;
901 rsq32 = dx32*dx32+dy32*dy32+dz32*dz32;
902 rsq33 = dx33*dx33+dy33*dy33+dz33*dz33;
903
904 rinv00 = 1.0/sqrt(rsq00);
905 rinv11 = 1.0/sqrt(rsq11);
906 rinv12 = 1.0/sqrt(rsq12);
907 rinv13 = 1.0/sqrt(rsq13);
908 rinv21 = 1.0/sqrt(rsq21);
909 rinv22 = 1.0/sqrt(rsq22);
910 rinv23 = 1.0/sqrt(rsq23);
911 rinv31 = 1.0/sqrt(rsq31);
912 rinv32 = 1.0/sqrt(rsq32);
913 rinv33 = 1.0/sqrt(rsq33);
914
915 rinvsq00 = rinv00*rinv00;
916 rinvsq11 = rinv11*rinv11;
917 rinvsq12 = rinv12*rinv12;
918 rinvsq13 = rinv13*rinv13;
919 rinvsq21 = rinv21*rinv21;
920 rinvsq22 = rinv22*rinv22;
921 rinvsq23 = rinv23*rinv23;
922 rinvsq31 = rinv31*rinv31;
923 rinvsq32 = rinv32*rinv32;
924 rinvsq33 = rinv33*rinv33;
925
926 /**************************
927 * CALCULATE INTERACTIONS *
928 **************************/
929
930 if (rsq00<rcutoff2)
931 {
932
933 r00 = rsq00*rinv00;
934
935 /* LENNARD-JONES DISPERSION/REPULSION */
936
937 rinvsix = rinvsq00*rinvsq00*rinvsq00;
938 vvdw6 = c6_00*rinvsix;
939 vvdw12 = c12_00*rinvsix*rinvsix;
940 vvdw = vvdw12*(1.0/12.0) - vvdw6*(1.0/6.0);
941 fvdw = (vvdw12-vvdw6)*rinvsq00;
942
943 d = r00-rswitch;
944 d = (d>0.0) ? d : 0.0;
945 d2 = d*d;
946 sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
947
948 dsw = d2*(swF2+d*(swF3+d*swF4));
949
950 /* Evaluate switch function */
951 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
952 fvdw = fvdw*sw - rinv00*vvdw*dsw;
953
954 fscal = fvdw;
955
956 /* Calculate temporary vectorial force */
957 tx = fscal*dx00;
958 ty = fscal*dy00;
959 tz = fscal*dz00;
960
961 /* Update vectorial force */
962 fix0 += tx;
963 fiy0 += ty;
964 fiz0 += tz;
965 f[j_coord_offset+DIM*0+XX] -= tx;
966 f[j_coord_offset+DIM*0+YY] -= ty;
967 f[j_coord_offset+DIM*0+ZZ] -= tz;
968
969 }
970
971 /**************************
972 * CALCULATE INTERACTIONS *
973 **************************/
974
975 if (rsq11<rcutoff2)
976 {
977
978 /* REACTION-FIELD ELECTROSTATICS */
979 felec = qq11*(rinv11*rinvsq11-krf2);
980
981 fscal = felec;
982
983 /* Calculate temporary vectorial force */
984 tx = fscal*dx11;
985 ty = fscal*dy11;
986 tz = fscal*dz11;
987
988 /* Update vectorial force */
989 fix1 += tx;
990 fiy1 += ty;
991 fiz1 += tz;
992 f[j_coord_offset+DIM*1+XX] -= tx;
993 f[j_coord_offset+DIM*1+YY] -= ty;
994 f[j_coord_offset+DIM*1+ZZ] -= tz;
995
996 }
997
998 /**************************
999 * CALCULATE INTERACTIONS *
1000 **************************/
1001
1002 if (rsq12<rcutoff2)
1003 {
1004
1005 /* REACTION-FIELD ELECTROSTATICS */
1006 felec = qq12*(rinv12*rinvsq12-krf2);
1007
1008 fscal = felec;
1009
1010 /* Calculate temporary vectorial force */
1011 tx = fscal*dx12;
1012 ty = fscal*dy12;
1013 tz = fscal*dz12;
1014
1015 /* Update vectorial force */
1016 fix1 += tx;
1017 fiy1 += ty;
1018 fiz1 += tz;
1019 f[j_coord_offset+DIM*2+XX] -= tx;
1020 f[j_coord_offset+DIM*2+YY] -= ty;
1021 f[j_coord_offset+DIM*2+ZZ] -= tz;
1022
1023 }
1024
1025 /**************************
1026 * CALCULATE INTERACTIONS *
1027 **************************/
1028
1029 if (rsq13<rcutoff2)
1030 {
1031
1032 /* REACTION-FIELD ELECTROSTATICS */
1033 felec = qq13*(rinv13*rinvsq13-krf2);
1034
1035 fscal = felec;
1036
1037 /* Calculate temporary vectorial force */
1038 tx = fscal*dx13;
1039 ty = fscal*dy13;
1040 tz = fscal*dz13;
1041
1042 /* Update vectorial force */
1043 fix1 += tx;
1044 fiy1 += ty;
1045 fiz1 += tz;
1046 f[j_coord_offset+DIM*3+XX] -= tx;
1047 f[j_coord_offset+DIM*3+YY] -= ty;
1048 f[j_coord_offset+DIM*3+ZZ] -= tz;
1049
1050 }
1051
1052 /**************************
1053 * CALCULATE INTERACTIONS *
1054 **************************/
1055
1056 if (rsq21<rcutoff2)
1057 {
1058
1059 /* REACTION-FIELD ELECTROSTATICS */
1060 felec = qq21*(rinv21*rinvsq21-krf2);
1061
1062 fscal = felec;
1063
1064 /* Calculate temporary vectorial force */
1065 tx = fscal*dx21;
1066 ty = fscal*dy21;
1067 tz = fscal*dz21;
1068
1069 /* Update vectorial force */
1070 fix2 += tx;
1071 fiy2 += ty;
1072 fiz2 += tz;
1073 f[j_coord_offset+DIM*1+XX] -= tx;
1074 f[j_coord_offset+DIM*1+YY] -= ty;
1075 f[j_coord_offset+DIM*1+ZZ] -= tz;
1076
1077 }
1078
1079 /**************************
1080 * CALCULATE INTERACTIONS *
1081 **************************/
1082
1083 if (rsq22<rcutoff2)
1084 {
1085
1086 /* REACTION-FIELD ELECTROSTATICS */
1087 felec = qq22*(rinv22*rinvsq22-krf2);
1088
1089 fscal = felec;
1090
1091 /* Calculate temporary vectorial force */
1092 tx = fscal*dx22;
1093 ty = fscal*dy22;
1094 tz = fscal*dz22;
1095
1096 /* Update vectorial force */
1097 fix2 += tx;
1098 fiy2 += ty;
1099 fiz2 += tz;
1100 f[j_coord_offset+DIM*2+XX] -= tx;
1101 f[j_coord_offset+DIM*2+YY] -= ty;
1102 f[j_coord_offset+DIM*2+ZZ] -= tz;
1103
1104 }
1105
1106 /**************************
1107 * CALCULATE INTERACTIONS *
1108 **************************/
1109
1110 if (rsq23<rcutoff2)
1111 {
1112
1113 /* REACTION-FIELD ELECTROSTATICS */
1114 felec = qq23*(rinv23*rinvsq23-krf2);
1115
1116 fscal = felec;
1117
1118 /* Calculate temporary vectorial force */
1119 tx = fscal*dx23;
1120 ty = fscal*dy23;
1121 tz = fscal*dz23;
1122
1123 /* Update vectorial force */
1124 fix2 += tx;
1125 fiy2 += ty;
1126 fiz2 += tz;
1127 f[j_coord_offset+DIM*3+XX] -= tx;
1128 f[j_coord_offset+DIM*3+YY] -= ty;
1129 f[j_coord_offset+DIM*3+ZZ] -= tz;
1130
1131 }
1132
1133 /**************************
1134 * CALCULATE INTERACTIONS *
1135 **************************/
1136
1137 if (rsq31<rcutoff2)
1138 {
1139
1140 /* REACTION-FIELD ELECTROSTATICS */
1141 felec = qq31*(rinv31*rinvsq31-krf2);
1142
1143 fscal = felec;
1144
1145 /* Calculate temporary vectorial force */
1146 tx = fscal*dx31;
1147 ty = fscal*dy31;
1148 tz = fscal*dz31;
1149
1150 /* Update vectorial force */
1151 fix3 += tx;
1152 fiy3 += ty;
1153 fiz3 += tz;
1154 f[j_coord_offset+DIM*1+XX] -= tx;
1155 f[j_coord_offset+DIM*1+YY] -= ty;
1156 f[j_coord_offset+DIM*1+ZZ] -= tz;
1157
1158 }
1159
1160 /**************************
1161 * CALCULATE INTERACTIONS *
1162 **************************/
1163
1164 if (rsq32<rcutoff2)
1165 {
1166
1167 /* REACTION-FIELD ELECTROSTATICS */
1168 felec = qq32*(rinv32*rinvsq32-krf2);
1169
1170 fscal = felec;
1171
1172 /* Calculate temporary vectorial force */
1173 tx = fscal*dx32;
1174 ty = fscal*dy32;
1175 tz = fscal*dz32;
1176
1177 /* Update vectorial force */
1178 fix3 += tx;
1179 fiy3 += ty;
1180 fiz3 += tz;
1181 f[j_coord_offset+DIM*2+XX] -= tx;
1182 f[j_coord_offset+DIM*2+YY] -= ty;
1183 f[j_coord_offset+DIM*2+ZZ] -= tz;
1184
1185 }
1186
1187 /**************************
1188 * CALCULATE INTERACTIONS *
1189 **************************/
1190
1191 if (rsq33<rcutoff2)
1192 {
1193
1194 /* REACTION-FIELD ELECTROSTATICS */
1195 felec = qq33*(rinv33*rinvsq33-krf2);
1196
1197 fscal = felec;
1198
1199 /* Calculate temporary vectorial force */
1200 tx = fscal*dx33;
1201 ty = fscal*dy33;
1202 tz = fscal*dz33;
1203
1204 /* Update vectorial force */
1205 fix3 += tx;
1206 fiy3 += ty;
1207 fiz3 += tz;
1208 f[j_coord_offset+DIM*3+XX] -= tx;
1209 f[j_coord_offset+DIM*3+YY] -= ty;
1210 f[j_coord_offset+DIM*3+ZZ] -= tz;
1211
1212 }
1213
1214 /* Inner loop uses 285 flops */
1215 }
1216 /* End of innermost loop */
1217
1218 tx = ty = tz = 0;
1219 f[i_coord_offset+DIM*0+XX] += fix0;
1220 f[i_coord_offset+DIM*0+YY] += fiy0;
1221 f[i_coord_offset+DIM*0+ZZ] += fiz0;
1222 tx += fix0;
1223 ty += fiy0;
1224 tz += fiz0;
1225 f[i_coord_offset+DIM*1+XX] += fix1;
1226 f[i_coord_offset+DIM*1+YY] += fiy1;
1227 f[i_coord_offset+DIM*1+ZZ] += fiz1;
1228 tx += fix1;
1229 ty += fiy1;
1230 tz += fiz1;
1231 f[i_coord_offset+DIM*2+XX] += fix2;
1232 f[i_coord_offset+DIM*2+YY] += fiy2;
1233 f[i_coord_offset+DIM*2+ZZ] += fiz2;
1234 tx += fix2;
1235 ty += fiy2;
1236 tz += fiz2;
1237 f[i_coord_offset+DIM*3+XX] += fix3;
1238 f[i_coord_offset+DIM*3+YY] += fiy3;
1239 f[i_coord_offset+DIM*3+ZZ] += fiz3;
1240 tx += fix3;
1241 ty += fiy3;
1242 tz += fiz3;
1243 fshift[i_shift_offset+XX] += tx;
1244 fshift[i_shift_offset+YY] += ty;
1245 fshift[i_shift_offset+ZZ] += tz;
1246
1247 /* Increment number of inner iterations */
1248 inneriter += j_index_end - j_index_start;
1249
1250 /* Outer loop uses 39 flops */
1251 }
1252
1253 /* Increment number of outer iterations */
1254 outeriter += nri;
1255
1256 /* Update outer/inner flops */
1257
1258 inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4W4_F,outeriter*39 + inneriter*285);
1259 }
The ultimate molecular dynamics simulation package