2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
36 * Note: this file was generated by the GROMACS c kernel generator.
44 #include "../nb_kernel.h"
45 #include "gromacs/legacyheaders/types/simple.h"
46 #include "gromacs/math/vec.h"
47 #include "gromacs/legacyheaders/nrnb.h"
50 * Gromacs nonbonded kernel: nb_kernel_ElecEwSw_VdwLJSw_GeomW3W3_VF_c
51 * Electrostatics interaction: Ewald
52 * VdW interaction: LennardJones
53 * Geometry: Water3-Water3
54 * Calculate force/pot: PotentialAndForce
57 nb_kernel_ElecEwSw_VdwLJSw_GeomW3W3_VF_c
58 (t_nblist * gmx_restrict nlist,
59 rvec * gmx_restrict xx,
60 rvec * gmx_restrict ff,
61 t_forcerec * gmx_restrict fr,
62 t_mdatoms * gmx_restrict mdatoms,
63 nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
64 t_nrnb * gmx_restrict nrnb)
66 int i_shift_offset,i_coord_offset,j_coord_offset;
67 int j_index_start,j_index_end;
68 int nri,inr,ggid,iidx,jidx,jnr,outeriter,inneriter;
69 real shX,shY,shZ,tx,ty,tz,fscal,rcutoff,rcutoff2;
70 int *iinr,*jindex,*jjnr,*shiftidx,*gid;
71 real *shiftvec,*fshift,*x,*f;
73 real ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
75 real ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
77 real ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
79 real jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
81 real jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
83 real jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
84 real dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00,cexp1_00,cexp2_00;
85 real dx01,dy01,dz01,rsq01,rinv01,rinvsq01,r01,qq01,c6_01,c12_01,cexp1_01,cexp2_01;
86 real dx02,dy02,dz02,rsq02,rinv02,rinvsq02,r02,qq02,c6_02,c12_02,cexp1_02,cexp2_02;
87 real dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10,cexp1_10,cexp2_10;
88 real dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11,cexp1_11,cexp2_11;
89 real dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12,cexp1_12,cexp2_12;
90 real dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20,cexp1_20,cexp2_20;
91 real dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21,cexp1_21,cexp2_21;
92 real dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22,cexp1_22,cexp2_22;
93 real velec,felec,velecsum,facel,crf,krf,krf2;
96 real rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,br,vvdwexp,sh_vdw_invrcut6;
100 real ewtabscale,eweps,sh_ewald,ewrt,ewtabhalfspace;
102 real rswitch,swV3,swV4,swV5,swF2,swF3,swF4,d,d2,sw,dsw;
109 jindex = nlist->jindex;
111 shiftidx = nlist->shift;
113 shiftvec = fr->shift_vec[0];
114 fshift = fr->fshift[0];
116 charge = mdatoms->chargeA;
117 nvdwtype = fr->ntype;
119 vdwtype = mdatoms->typeA;
121 sh_ewald = fr->ic->sh_ewald;
122 ewtab = fr->ic->tabq_coul_FDV0;
123 ewtabscale = fr->ic->tabq_scale;
124 ewtabhalfspace = 0.5/ewtabscale;
126 /* Setup water-specific parameters */
127 inr = nlist->iinr[0];
128 iq0 = facel*charge[inr+0];
129 iq1 = facel*charge[inr+1];
130 iq2 = facel*charge[inr+2];
131 vdwioffset0 = 2*nvdwtype*vdwtype[inr+0];
136 vdwjidx0 = 2*vdwtype[inr+0];
138 c6_00 = vdwparam[vdwioffset0+vdwjidx0];
139 c12_00 = vdwparam[vdwioffset0+vdwjidx0+1];
149 /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
150 rcutoff = fr->rcoulomb;
151 rcutoff2 = rcutoff*rcutoff;
153 rswitch = fr->rcoulomb_switch;
154 /* Setup switch parameters */
156 swV3 = -10.0/(d*d*d);
157 swV4 = 15.0/(d*d*d*d);
158 swV5 = -6.0/(d*d*d*d*d);
159 swF2 = -30.0/(d*d*d);
160 swF3 = 60.0/(d*d*d*d);
161 swF4 = -30.0/(d*d*d*d*d);
166 /* Start outer loop over neighborlists */
167 for(iidx=0; iidx<nri; iidx++)
169 /* Load shift vector for this list */
170 i_shift_offset = DIM*shiftidx[iidx];
171 shX = shiftvec[i_shift_offset+XX];
172 shY = shiftvec[i_shift_offset+YY];
173 shZ = shiftvec[i_shift_offset+ZZ];
175 /* Load limits for loop over neighbors */
176 j_index_start = jindex[iidx];
177 j_index_end = jindex[iidx+1];
179 /* Get outer coordinate index */
181 i_coord_offset = DIM*inr;
183 /* Load i particle coords and add shift vector */
184 ix0 = shX + x[i_coord_offset+DIM*0+XX];
185 iy0 = shY + x[i_coord_offset+DIM*0+YY];
186 iz0 = shZ + x[i_coord_offset+DIM*0+ZZ];
187 ix1 = shX + x[i_coord_offset+DIM*1+XX];
188 iy1 = shY + x[i_coord_offset+DIM*1+YY];
189 iz1 = shZ + x[i_coord_offset+DIM*1+ZZ];
190 ix2 = shX + x[i_coord_offset+DIM*2+XX];
191 iy2 = shY + x[i_coord_offset+DIM*2+YY];
192 iz2 = shZ + x[i_coord_offset+DIM*2+ZZ];
204 /* Reset potential sums */
208 /* Start inner kernel loop */
209 for(jidx=j_index_start; jidx<j_index_end; jidx++)
211 /* Get j neighbor index, and coordinate index */
213 j_coord_offset = DIM*jnr;
215 /* load j atom coordinates */
216 jx0 = x[j_coord_offset+DIM*0+XX];
217 jy0 = x[j_coord_offset+DIM*0+YY];
218 jz0 = x[j_coord_offset+DIM*0+ZZ];
219 jx1 = x[j_coord_offset+DIM*1+XX];
220 jy1 = x[j_coord_offset+DIM*1+YY];
221 jz1 = x[j_coord_offset+DIM*1+ZZ];
222 jx2 = x[j_coord_offset+DIM*2+XX];
223 jy2 = x[j_coord_offset+DIM*2+YY];
224 jz2 = x[j_coord_offset+DIM*2+ZZ];
226 /* Calculate displacement vector */
255 /* Calculate squared distance and things based on it */
256 rsq00 = dx00*dx00+dy00*dy00+dz00*dz00;
257 rsq01 = dx01*dx01+dy01*dy01+dz01*dz01;
258 rsq02 = dx02*dx02+dy02*dy02+dz02*dz02;
259 rsq10 = dx10*dx10+dy10*dy10+dz10*dz10;
260 rsq11 = dx11*dx11+dy11*dy11+dz11*dz11;
261 rsq12 = dx12*dx12+dy12*dy12+dz12*dz12;
262 rsq20 = dx20*dx20+dy20*dy20+dz20*dz20;
263 rsq21 = dx21*dx21+dy21*dy21+dz21*dz21;
264 rsq22 = dx22*dx22+dy22*dy22+dz22*dz22;
266 rinv00 = gmx_invsqrt(rsq00);
267 rinv01 = gmx_invsqrt(rsq01);
268 rinv02 = gmx_invsqrt(rsq02);
269 rinv10 = gmx_invsqrt(rsq10);
270 rinv11 = gmx_invsqrt(rsq11);
271 rinv12 = gmx_invsqrt(rsq12);
272 rinv20 = gmx_invsqrt(rsq20);
273 rinv21 = gmx_invsqrt(rsq21);
274 rinv22 = gmx_invsqrt(rsq22);
276 rinvsq00 = rinv00*rinv00;
277 rinvsq01 = rinv01*rinv01;
278 rinvsq02 = rinv02*rinv02;
279 rinvsq10 = rinv10*rinv10;
280 rinvsq11 = rinv11*rinv11;
281 rinvsq12 = rinv12*rinv12;
282 rinvsq20 = rinv20*rinv20;
283 rinvsq21 = rinv21*rinv21;
284 rinvsq22 = rinv22*rinv22;
286 /**************************
287 * CALCULATE INTERACTIONS *
288 **************************/
295 /* EWALD ELECTROSTATICS */
297 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
298 ewrt = r00*ewtabscale;
302 felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
303 velec = qq00*(rinv00-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
304 felec = qq00*rinv00*(rinvsq00-felec);
306 /* LENNARD-JONES DISPERSION/REPULSION */
308 rinvsix = rinvsq00*rinvsq00*rinvsq00;
309 vvdw6 = c6_00*rinvsix;
310 vvdw12 = c12_00*rinvsix*rinvsix;
311 vvdw = vvdw12*(1.0/12.0) - vvdw6*(1.0/6.0);
312 fvdw = (vvdw12-vvdw6)*rinvsq00;
315 d = (d>0.0) ? d : 0.0;
317 sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
319 dsw = d2*(swF2+d*(swF3+d*swF4));
321 /* Evaluate switch function */
322 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
323 felec = felec*sw - rinv00*velec*dsw;
324 fvdw = fvdw*sw - rinv00*vvdw*dsw;
328 /* Update potential sums from outer loop */
334 /* Calculate temporary vectorial force */
339 /* Update vectorial force */
343 f[j_coord_offset+DIM*0+XX] -= tx;
344 f[j_coord_offset+DIM*0+YY] -= ty;
345 f[j_coord_offset+DIM*0+ZZ] -= tz;
349 /**************************
350 * CALCULATE INTERACTIONS *
351 **************************/
358 /* EWALD ELECTROSTATICS */
360 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
361 ewrt = r01*ewtabscale;
365 felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
366 velec = qq01*(rinv01-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
367 felec = qq01*rinv01*(rinvsq01-felec);
370 d = (d>0.0) ? d : 0.0;
372 sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
374 dsw = d2*(swF2+d*(swF3+d*swF4));
376 /* Evaluate switch function */
377 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
378 felec = felec*sw - rinv01*velec*dsw;
381 /* Update potential sums from outer loop */
386 /* Calculate temporary vectorial force */
391 /* Update vectorial force */
395 f[j_coord_offset+DIM*1+XX] -= tx;
396 f[j_coord_offset+DIM*1+YY] -= ty;
397 f[j_coord_offset+DIM*1+ZZ] -= tz;
401 /**************************
402 * CALCULATE INTERACTIONS *
403 **************************/
410 /* EWALD ELECTROSTATICS */
412 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
413 ewrt = r02*ewtabscale;
417 felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
418 velec = qq02*(rinv02-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
419 felec = qq02*rinv02*(rinvsq02-felec);
422 d = (d>0.0) ? d : 0.0;
424 sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
426 dsw = d2*(swF2+d*(swF3+d*swF4));
428 /* Evaluate switch function */
429 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
430 felec = felec*sw - rinv02*velec*dsw;
433 /* Update potential sums from outer loop */
438 /* Calculate temporary vectorial force */
443 /* Update vectorial force */
447 f[j_coord_offset+DIM*2+XX] -= tx;
448 f[j_coord_offset+DIM*2+YY] -= ty;
449 f[j_coord_offset+DIM*2+ZZ] -= tz;
453 /**************************
454 * CALCULATE INTERACTIONS *
455 **************************/
462 /* EWALD ELECTROSTATICS */
464 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
465 ewrt = r10*ewtabscale;
469 felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
470 velec = qq10*(rinv10-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
471 felec = qq10*rinv10*(rinvsq10-felec);
474 d = (d>0.0) ? d : 0.0;
476 sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
478 dsw = d2*(swF2+d*(swF3+d*swF4));
480 /* Evaluate switch function */
481 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
482 felec = felec*sw - rinv10*velec*dsw;
485 /* Update potential sums from outer loop */
490 /* Calculate temporary vectorial force */
495 /* Update vectorial force */
499 f[j_coord_offset+DIM*0+XX] -= tx;
500 f[j_coord_offset+DIM*0+YY] -= ty;
501 f[j_coord_offset+DIM*0+ZZ] -= tz;
505 /**************************
506 * CALCULATE INTERACTIONS *
507 **************************/
514 /* EWALD ELECTROSTATICS */
516 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
517 ewrt = r11*ewtabscale;
521 felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
522 velec = qq11*(rinv11-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
523 felec = qq11*rinv11*(rinvsq11-felec);
526 d = (d>0.0) ? d : 0.0;
528 sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
530 dsw = d2*(swF2+d*(swF3+d*swF4));
532 /* Evaluate switch function */
533 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
534 felec = felec*sw - rinv11*velec*dsw;
537 /* Update potential sums from outer loop */
542 /* Calculate temporary vectorial force */
547 /* Update vectorial force */
551 f[j_coord_offset+DIM*1+XX] -= tx;
552 f[j_coord_offset+DIM*1+YY] -= ty;
553 f[j_coord_offset+DIM*1+ZZ] -= tz;
557 /**************************
558 * CALCULATE INTERACTIONS *
559 **************************/
566 /* EWALD ELECTROSTATICS */
568 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
569 ewrt = r12*ewtabscale;
573 felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
574 velec = qq12*(rinv12-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
575 felec = qq12*rinv12*(rinvsq12-felec);
578 d = (d>0.0) ? d : 0.0;
580 sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
582 dsw = d2*(swF2+d*(swF3+d*swF4));
584 /* Evaluate switch function */
585 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
586 felec = felec*sw - rinv12*velec*dsw;
589 /* Update potential sums from outer loop */
594 /* Calculate temporary vectorial force */
599 /* Update vectorial force */
603 f[j_coord_offset+DIM*2+XX] -= tx;
604 f[j_coord_offset+DIM*2+YY] -= ty;
605 f[j_coord_offset+DIM*2+ZZ] -= tz;
609 /**************************
610 * CALCULATE INTERACTIONS *
611 **************************/
618 /* EWALD ELECTROSTATICS */
620 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
621 ewrt = r20*ewtabscale;
625 felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
626 velec = qq20*(rinv20-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
627 felec = qq20*rinv20*(rinvsq20-felec);
630 d = (d>0.0) ? d : 0.0;
632 sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
634 dsw = d2*(swF2+d*(swF3+d*swF4));
636 /* Evaluate switch function */
637 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
638 felec = felec*sw - rinv20*velec*dsw;
641 /* Update potential sums from outer loop */
646 /* Calculate temporary vectorial force */
651 /* Update vectorial force */
655 f[j_coord_offset+DIM*0+XX] -= tx;
656 f[j_coord_offset+DIM*0+YY] -= ty;
657 f[j_coord_offset+DIM*0+ZZ] -= tz;
661 /**************************
662 * CALCULATE INTERACTIONS *
663 **************************/
670 /* EWALD ELECTROSTATICS */
672 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
673 ewrt = r21*ewtabscale;
677 felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
678 velec = qq21*(rinv21-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
679 felec = qq21*rinv21*(rinvsq21-felec);
682 d = (d>0.0) ? d : 0.0;
684 sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
686 dsw = d2*(swF2+d*(swF3+d*swF4));
688 /* Evaluate switch function */
689 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
690 felec = felec*sw - rinv21*velec*dsw;
693 /* Update potential sums from outer loop */
698 /* Calculate temporary vectorial force */
703 /* Update vectorial force */
707 f[j_coord_offset+DIM*1+XX] -= tx;
708 f[j_coord_offset+DIM*1+YY] -= ty;
709 f[j_coord_offset+DIM*1+ZZ] -= tz;
713 /**************************
714 * CALCULATE INTERACTIONS *
715 **************************/
722 /* EWALD ELECTROSTATICS */
724 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
725 ewrt = r22*ewtabscale;
729 felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
730 velec = qq22*(rinv22-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
731 felec = qq22*rinv22*(rinvsq22-felec);
734 d = (d>0.0) ? d : 0.0;
736 sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
738 dsw = d2*(swF2+d*(swF3+d*swF4));
740 /* Evaluate switch function */
741 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
742 felec = felec*sw - rinv22*velec*dsw;
745 /* Update potential sums from outer loop */
750 /* Calculate temporary vectorial force */
755 /* Update vectorial force */
759 f[j_coord_offset+DIM*2+XX] -= tx;
760 f[j_coord_offset+DIM*2+YY] -= ty;
761 f[j_coord_offset+DIM*2+ZZ] -= tz;
765 /* Inner loop uses 538 flops */
767 /* End of innermost loop */
770 f[i_coord_offset+DIM*0+XX] += fix0;
771 f[i_coord_offset+DIM*0+YY] += fiy0;
772 f[i_coord_offset+DIM*0+ZZ] += fiz0;
776 f[i_coord_offset+DIM*1+XX] += fix1;
777 f[i_coord_offset+DIM*1+YY] += fiy1;
778 f[i_coord_offset+DIM*1+ZZ] += fiz1;
782 f[i_coord_offset+DIM*2+XX] += fix2;
783 f[i_coord_offset+DIM*2+YY] += fiy2;
784 f[i_coord_offset+DIM*2+ZZ] += fiz2;
788 fshift[i_shift_offset+XX] += tx;
789 fshift[i_shift_offset+YY] += ty;
790 fshift[i_shift_offset+ZZ] += tz;
793 /* Update potential energies */
794 kernel_data->energygrp_elec[ggid] += velecsum;
795 kernel_data->energygrp_vdw[ggid] += vvdwsum;
797 /* Increment number of inner iterations */
798 inneriter += j_index_end - j_index_start;
800 /* Outer loop uses 32 flops */
803 /* Increment number of outer iterations */
806 /* Update outer/inner flops */
808 inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W3W3_VF,outeriter*32 + inneriter*538);
811 * Gromacs nonbonded kernel: nb_kernel_ElecEwSw_VdwLJSw_GeomW3W3_F_c
812 * Electrostatics interaction: Ewald
813 * VdW interaction: LennardJones
814 * Geometry: Water3-Water3
815 * Calculate force/pot: Force
818 nb_kernel_ElecEwSw_VdwLJSw_GeomW3W3_F_c
819 (t_nblist * gmx_restrict nlist,
820 rvec * gmx_restrict xx,
821 rvec * gmx_restrict ff,
822 t_forcerec * gmx_restrict fr,
823 t_mdatoms * gmx_restrict mdatoms,
824 nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
825 t_nrnb * gmx_restrict nrnb)
827 int i_shift_offset,i_coord_offset,j_coord_offset;
828 int j_index_start,j_index_end;
829 int nri,inr,ggid,iidx,jidx,jnr,outeriter,inneriter;
830 real shX,shY,shZ,tx,ty,tz,fscal,rcutoff,rcutoff2;
831 int *iinr,*jindex,*jjnr,*shiftidx,*gid;
832 real *shiftvec,*fshift,*x,*f;
834 real ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
836 real ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
838 real ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
840 real jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
842 real jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
844 real jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
845 real dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00,cexp1_00,cexp2_00;
846 real dx01,dy01,dz01,rsq01,rinv01,rinvsq01,r01,qq01,c6_01,c12_01,cexp1_01,cexp2_01;
847 real dx02,dy02,dz02,rsq02,rinv02,rinvsq02,r02,qq02,c6_02,c12_02,cexp1_02,cexp2_02;
848 real dx10,dy10,dz10,rsq10,rinv10,rinvsq10,r10,qq10,c6_10,c12_10,cexp1_10,cexp2_10;
849 real dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11,cexp1_11,cexp2_11;
850 real dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12,cexp1_12,cexp2_12;
851 real dx20,dy20,dz20,rsq20,rinv20,rinvsq20,r20,qq20,c6_20,c12_20,cexp1_20,cexp2_20;
852 real dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21,cexp1_21,cexp2_21;
853 real dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22,cexp1_22,cexp2_22;
854 real velec,felec,velecsum,facel,crf,krf,krf2;
857 real rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,br,vvdwexp,sh_vdw_invrcut6;
861 real ewtabscale,eweps,sh_ewald,ewrt,ewtabhalfspace;
863 real rswitch,swV3,swV4,swV5,swF2,swF3,swF4,d,d2,sw,dsw;
870 jindex = nlist->jindex;
872 shiftidx = nlist->shift;
874 shiftvec = fr->shift_vec[0];
875 fshift = fr->fshift[0];
877 charge = mdatoms->chargeA;
878 nvdwtype = fr->ntype;
880 vdwtype = mdatoms->typeA;
882 sh_ewald = fr->ic->sh_ewald;
883 ewtab = fr->ic->tabq_coul_FDV0;
884 ewtabscale = fr->ic->tabq_scale;
885 ewtabhalfspace = 0.5/ewtabscale;
887 /* Setup water-specific parameters */
888 inr = nlist->iinr[0];
889 iq0 = facel*charge[inr+0];
890 iq1 = facel*charge[inr+1];
891 iq2 = facel*charge[inr+2];
892 vdwioffset0 = 2*nvdwtype*vdwtype[inr+0];
897 vdwjidx0 = 2*vdwtype[inr+0];
899 c6_00 = vdwparam[vdwioffset0+vdwjidx0];
900 c12_00 = vdwparam[vdwioffset0+vdwjidx0+1];
910 /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
911 rcutoff = fr->rcoulomb;
912 rcutoff2 = rcutoff*rcutoff;
914 rswitch = fr->rcoulomb_switch;
915 /* Setup switch parameters */
917 swV3 = -10.0/(d*d*d);
918 swV4 = 15.0/(d*d*d*d);
919 swV5 = -6.0/(d*d*d*d*d);
920 swF2 = -30.0/(d*d*d);
921 swF3 = 60.0/(d*d*d*d);
922 swF4 = -30.0/(d*d*d*d*d);
927 /* Start outer loop over neighborlists */
928 for(iidx=0; iidx<nri; iidx++)
930 /* Load shift vector for this list */
931 i_shift_offset = DIM*shiftidx[iidx];
932 shX = shiftvec[i_shift_offset+XX];
933 shY = shiftvec[i_shift_offset+YY];
934 shZ = shiftvec[i_shift_offset+ZZ];
936 /* Load limits for loop over neighbors */
937 j_index_start = jindex[iidx];
938 j_index_end = jindex[iidx+1];
940 /* Get outer coordinate index */
942 i_coord_offset = DIM*inr;
944 /* Load i particle coords and add shift vector */
945 ix0 = shX + x[i_coord_offset+DIM*0+XX];
946 iy0 = shY + x[i_coord_offset+DIM*0+YY];
947 iz0 = shZ + x[i_coord_offset+DIM*0+ZZ];
948 ix1 = shX + x[i_coord_offset+DIM*1+XX];
949 iy1 = shY + x[i_coord_offset+DIM*1+YY];
950 iz1 = shZ + x[i_coord_offset+DIM*1+ZZ];
951 ix2 = shX + x[i_coord_offset+DIM*2+XX];
952 iy2 = shY + x[i_coord_offset+DIM*2+YY];
953 iz2 = shZ + x[i_coord_offset+DIM*2+ZZ];
965 /* Start inner kernel loop */
966 for(jidx=j_index_start; jidx<j_index_end; jidx++)
968 /* Get j neighbor index, and coordinate index */
970 j_coord_offset = DIM*jnr;
972 /* load j atom coordinates */
973 jx0 = x[j_coord_offset+DIM*0+XX];
974 jy0 = x[j_coord_offset+DIM*0+YY];
975 jz0 = x[j_coord_offset+DIM*0+ZZ];
976 jx1 = x[j_coord_offset+DIM*1+XX];
977 jy1 = x[j_coord_offset+DIM*1+YY];
978 jz1 = x[j_coord_offset+DIM*1+ZZ];
979 jx2 = x[j_coord_offset+DIM*2+XX];
980 jy2 = x[j_coord_offset+DIM*2+YY];
981 jz2 = x[j_coord_offset+DIM*2+ZZ];
983 /* Calculate displacement vector */
1012 /* Calculate squared distance and things based on it */
1013 rsq00 = dx00*dx00+dy00*dy00+dz00*dz00;
1014 rsq01 = dx01*dx01+dy01*dy01+dz01*dz01;
1015 rsq02 = dx02*dx02+dy02*dy02+dz02*dz02;
1016 rsq10 = dx10*dx10+dy10*dy10+dz10*dz10;
1017 rsq11 = dx11*dx11+dy11*dy11+dz11*dz11;
1018 rsq12 = dx12*dx12+dy12*dy12+dz12*dz12;
1019 rsq20 = dx20*dx20+dy20*dy20+dz20*dz20;
1020 rsq21 = dx21*dx21+dy21*dy21+dz21*dz21;
1021 rsq22 = dx22*dx22+dy22*dy22+dz22*dz22;
1023 rinv00 = gmx_invsqrt(rsq00);
1024 rinv01 = gmx_invsqrt(rsq01);
1025 rinv02 = gmx_invsqrt(rsq02);
1026 rinv10 = gmx_invsqrt(rsq10);
1027 rinv11 = gmx_invsqrt(rsq11);
1028 rinv12 = gmx_invsqrt(rsq12);
1029 rinv20 = gmx_invsqrt(rsq20);
1030 rinv21 = gmx_invsqrt(rsq21);
1031 rinv22 = gmx_invsqrt(rsq22);
1033 rinvsq00 = rinv00*rinv00;
1034 rinvsq01 = rinv01*rinv01;
1035 rinvsq02 = rinv02*rinv02;
1036 rinvsq10 = rinv10*rinv10;
1037 rinvsq11 = rinv11*rinv11;
1038 rinvsq12 = rinv12*rinv12;
1039 rinvsq20 = rinv20*rinv20;
1040 rinvsq21 = rinv21*rinv21;
1041 rinvsq22 = rinv22*rinv22;
1043 /**************************
1044 * CALCULATE INTERACTIONS *
1045 **************************/
1052 /* EWALD ELECTROSTATICS */
1054 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1055 ewrt = r00*ewtabscale;
1057 eweps = ewrt-ewitab;
1059 felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
1060 velec = qq00*(rinv00-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
1061 felec = qq00*rinv00*(rinvsq00-felec);
1063 /* LENNARD-JONES DISPERSION/REPULSION */
1065 rinvsix = rinvsq00*rinvsq00*rinvsq00;
1066 vvdw6 = c6_00*rinvsix;
1067 vvdw12 = c12_00*rinvsix*rinvsix;
1068 vvdw = vvdw12*(1.0/12.0) - vvdw6*(1.0/6.0);
1069 fvdw = (vvdw12-vvdw6)*rinvsq00;
1072 d = (d>0.0) ? d : 0.0;
1074 sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
1076 dsw = d2*(swF2+d*(swF3+d*swF4));
1078 /* Evaluate switch function */
1079 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1080 felec = felec*sw - rinv00*velec*dsw;
1081 fvdw = fvdw*sw - rinv00*vvdw*dsw;
1085 /* Calculate temporary vectorial force */
1090 /* Update vectorial force */
1094 f[j_coord_offset+DIM*0+XX] -= tx;
1095 f[j_coord_offset+DIM*0+YY] -= ty;
1096 f[j_coord_offset+DIM*0+ZZ] -= tz;
1100 /**************************
1101 * CALCULATE INTERACTIONS *
1102 **************************/
1109 /* EWALD ELECTROSTATICS */
1111 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1112 ewrt = r01*ewtabscale;
1114 eweps = ewrt-ewitab;
1116 felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
1117 velec = qq01*(rinv01-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
1118 felec = qq01*rinv01*(rinvsq01-felec);
1121 d = (d>0.0) ? d : 0.0;
1123 sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
1125 dsw = d2*(swF2+d*(swF3+d*swF4));
1127 /* Evaluate switch function */
1128 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1129 felec = felec*sw - rinv01*velec*dsw;
1133 /* Calculate temporary vectorial force */
1138 /* Update vectorial force */
1142 f[j_coord_offset+DIM*1+XX] -= tx;
1143 f[j_coord_offset+DIM*1+YY] -= ty;
1144 f[j_coord_offset+DIM*1+ZZ] -= tz;
1148 /**************************
1149 * CALCULATE INTERACTIONS *
1150 **************************/
1157 /* EWALD ELECTROSTATICS */
1159 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1160 ewrt = r02*ewtabscale;
1162 eweps = ewrt-ewitab;
1164 felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
1165 velec = qq02*(rinv02-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
1166 felec = qq02*rinv02*(rinvsq02-felec);
1169 d = (d>0.0) ? d : 0.0;
1171 sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
1173 dsw = d2*(swF2+d*(swF3+d*swF4));
1175 /* Evaluate switch function */
1176 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1177 felec = felec*sw - rinv02*velec*dsw;
1181 /* Calculate temporary vectorial force */
1186 /* Update vectorial force */
1190 f[j_coord_offset+DIM*2+XX] -= tx;
1191 f[j_coord_offset+DIM*2+YY] -= ty;
1192 f[j_coord_offset+DIM*2+ZZ] -= tz;
1196 /**************************
1197 * CALCULATE INTERACTIONS *
1198 **************************/
1205 /* EWALD ELECTROSTATICS */
1207 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1208 ewrt = r10*ewtabscale;
1210 eweps = ewrt-ewitab;
1212 felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
1213 velec = qq10*(rinv10-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
1214 felec = qq10*rinv10*(rinvsq10-felec);
1217 d = (d>0.0) ? d : 0.0;
1219 sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
1221 dsw = d2*(swF2+d*(swF3+d*swF4));
1223 /* Evaluate switch function */
1224 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1225 felec = felec*sw - rinv10*velec*dsw;
1229 /* Calculate temporary vectorial force */
1234 /* Update vectorial force */
1238 f[j_coord_offset+DIM*0+XX] -= tx;
1239 f[j_coord_offset+DIM*0+YY] -= ty;
1240 f[j_coord_offset+DIM*0+ZZ] -= tz;
1244 /**************************
1245 * CALCULATE INTERACTIONS *
1246 **************************/
1253 /* EWALD ELECTROSTATICS */
1255 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1256 ewrt = r11*ewtabscale;
1258 eweps = ewrt-ewitab;
1260 felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
1261 velec = qq11*(rinv11-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
1262 felec = qq11*rinv11*(rinvsq11-felec);
1265 d = (d>0.0) ? d : 0.0;
1267 sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
1269 dsw = d2*(swF2+d*(swF3+d*swF4));
1271 /* Evaluate switch function */
1272 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1273 felec = felec*sw - rinv11*velec*dsw;
1277 /* Calculate temporary vectorial force */
1282 /* Update vectorial force */
1286 f[j_coord_offset+DIM*1+XX] -= tx;
1287 f[j_coord_offset+DIM*1+YY] -= ty;
1288 f[j_coord_offset+DIM*1+ZZ] -= tz;
1292 /**************************
1293 * CALCULATE INTERACTIONS *
1294 **************************/
1301 /* EWALD ELECTROSTATICS */
1303 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1304 ewrt = r12*ewtabscale;
1306 eweps = ewrt-ewitab;
1308 felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
1309 velec = qq12*(rinv12-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
1310 felec = qq12*rinv12*(rinvsq12-felec);
1313 d = (d>0.0) ? d : 0.0;
1315 sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
1317 dsw = d2*(swF2+d*(swF3+d*swF4));
1319 /* Evaluate switch function */
1320 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1321 felec = felec*sw - rinv12*velec*dsw;
1325 /* Calculate temporary vectorial force */
1330 /* Update vectorial force */
1334 f[j_coord_offset+DIM*2+XX] -= tx;
1335 f[j_coord_offset+DIM*2+YY] -= ty;
1336 f[j_coord_offset+DIM*2+ZZ] -= tz;
1340 /**************************
1341 * CALCULATE INTERACTIONS *
1342 **************************/
1349 /* EWALD ELECTROSTATICS */
1351 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1352 ewrt = r20*ewtabscale;
1354 eweps = ewrt-ewitab;
1356 felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
1357 velec = qq20*(rinv20-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
1358 felec = qq20*rinv20*(rinvsq20-felec);
1361 d = (d>0.0) ? d : 0.0;
1363 sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
1365 dsw = d2*(swF2+d*(swF3+d*swF4));
1367 /* Evaluate switch function */
1368 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1369 felec = felec*sw - rinv20*velec*dsw;
1373 /* Calculate temporary vectorial force */
1378 /* Update vectorial force */
1382 f[j_coord_offset+DIM*0+XX] -= tx;
1383 f[j_coord_offset+DIM*0+YY] -= ty;
1384 f[j_coord_offset+DIM*0+ZZ] -= tz;
1388 /**************************
1389 * CALCULATE INTERACTIONS *
1390 **************************/
1397 /* EWALD ELECTROSTATICS */
1399 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1400 ewrt = r21*ewtabscale;
1402 eweps = ewrt-ewitab;
1404 felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
1405 velec = qq21*(rinv21-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
1406 felec = qq21*rinv21*(rinvsq21-felec);
1409 d = (d>0.0) ? d : 0.0;
1411 sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
1413 dsw = d2*(swF2+d*(swF3+d*swF4));
1415 /* Evaluate switch function */
1416 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1417 felec = felec*sw - rinv21*velec*dsw;
1421 /* Calculate temporary vectorial force */
1426 /* Update vectorial force */
1430 f[j_coord_offset+DIM*1+XX] -= tx;
1431 f[j_coord_offset+DIM*1+YY] -= ty;
1432 f[j_coord_offset+DIM*1+ZZ] -= tz;
1436 /**************************
1437 * CALCULATE INTERACTIONS *
1438 **************************/
1445 /* EWALD ELECTROSTATICS */
1447 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1448 ewrt = r22*ewtabscale;
1450 eweps = ewrt-ewitab;
1452 felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
1453 velec = qq22*(rinv22-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
1454 felec = qq22*rinv22*(rinvsq22-felec);
1457 d = (d>0.0) ? d : 0.0;
1459 sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
1461 dsw = d2*(swF2+d*(swF3+d*swF4));
1463 /* Evaluate switch function */
1464 /* fscal'=f'/r=-(v*sw)'/r=-(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=fscal*sw-v*dsw/r */
1465 felec = felec*sw - rinv22*velec*dsw;
1469 /* Calculate temporary vectorial force */
1474 /* Update vectorial force */
1478 f[j_coord_offset+DIM*2+XX] -= tx;
1479 f[j_coord_offset+DIM*2+YY] -= ty;
1480 f[j_coord_offset+DIM*2+ZZ] -= tz;
1484 /* Inner loop uses 518 flops */
1486 /* End of innermost loop */
1489 f[i_coord_offset+DIM*0+XX] += fix0;
1490 f[i_coord_offset+DIM*0+YY] += fiy0;
1491 f[i_coord_offset+DIM*0+ZZ] += fiz0;
1495 f[i_coord_offset+DIM*1+XX] += fix1;
1496 f[i_coord_offset+DIM*1+YY] += fiy1;
1497 f[i_coord_offset+DIM*1+ZZ] += fiz1;
1501 f[i_coord_offset+DIM*2+XX] += fix2;
1502 f[i_coord_offset+DIM*2+YY] += fiy2;
1503 f[i_coord_offset+DIM*2+ZZ] += fiz2;
1507 fshift[i_shift_offset+XX] += tx;
1508 fshift[i_shift_offset+YY] += ty;
1509 fshift[i_shift_offset+ZZ] += tz;
1511 /* Increment number of inner iterations */
1512 inneriter += j_index_end - j_index_start;
1514 /* Outer loop uses 30 flops */
1517 /* Increment number of outer iterations */
1520 /* Update outer/inner flops */
1522 inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W3W3_F,outeriter*30 + inneriter*518);