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.
42 #include "../nb_kernel.h"
43 #include "gromacs/legacyheaders/types/simple.h"
44 #include "gromacs/math/vec.h"
45 #include "gromacs/legacyheaders/nrnb.h"
48 * Gromacs nonbonded kernel: nb_kernel_ElecRFCut_VdwBhamSh_GeomW4W4_VF_c
49 * Electrostatics interaction: ReactionField
50 * VdW interaction: Buckingham
51 * Geometry: Water4-Water4
52 * Calculate force/pot: PotentialAndForce
55 nb_kernel_ElecRFCut_VdwBhamSh_GeomW4W4_VF_c
56 (t_nblist * gmx_restrict nlist,
57 rvec * gmx_restrict xx,
58 rvec * gmx_restrict ff,
59 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)
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;
71 real ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
73 real ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
75 real ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
77 real ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
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;
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;
99 real rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,br,vvdwexp,sh_vdw_invrcut6;
108 jindex = nlist->jindex;
110 shiftidx = nlist->shift;
112 shiftvec = fr->shift_vec[0];
113 fshift = fr->fshift[0];
115 charge = mdatoms->chargeA;
119 nvdwtype = fr->ntype;
121 vdwtype = mdatoms->typeA;
123 /* Setup water-specific parameters */
124 inr = nlist->iinr[0];
125 iq1 = facel*charge[inr+1];
126 iq2 = facel*charge[inr+2];
127 iq3 = facel*charge[inr+3];
128 vdwioffset0 = 3*nvdwtype*vdwtype[inr+0];
133 vdwjidx0 = 3*vdwtype[inr+0];
134 c6_00 = vdwparam[vdwioffset0+vdwjidx0];
135 cexp1_00 = vdwparam[vdwioffset0+vdwjidx0+1];
136 cexp2_00 = vdwparam[vdwioffset0+vdwjidx0+2];
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->rcoulomb;
149 rcutoff2 = rcutoff*rcutoff;
151 sh_vdw_invrcut6 = fr->ic->sh_invrc6;
157 /* Start outer loop over neighborlists */
158 for(iidx=0; iidx<nri; iidx++)
160 /* Load shift vector for this list */
161 i_shift_offset = DIM*shiftidx[iidx];
162 shX = shiftvec[i_shift_offset+XX];
163 shY = shiftvec[i_shift_offset+YY];
164 shZ = shiftvec[i_shift_offset+ZZ];
166 /* Load limits for loop over neighbors */
167 j_index_start = jindex[iidx];
168 j_index_end = jindex[iidx+1];
170 /* Get outer coordinate index */
172 i_coord_offset = DIM*inr;
174 /* Load i particle coords and add shift vector */
175 ix0 = shX + x[i_coord_offset+DIM*0+XX];
176 iy0 = shY + x[i_coord_offset+DIM*0+YY];
177 iz0 = shZ + x[i_coord_offset+DIM*0+ZZ];
178 ix1 = shX + x[i_coord_offset+DIM*1+XX];
179 iy1 = shY + x[i_coord_offset+DIM*1+YY];
180 iz1 = shZ + x[i_coord_offset+DIM*1+ZZ];
181 ix2 = shX + x[i_coord_offset+DIM*2+XX];
182 iy2 = shY + x[i_coord_offset+DIM*2+YY];
183 iz2 = shZ + x[i_coord_offset+DIM*2+ZZ];
184 ix3 = shX + x[i_coord_offset+DIM*3+XX];
185 iy3 = shY + x[i_coord_offset+DIM*3+YY];
186 iz3 = shZ + x[i_coord_offset+DIM*3+ZZ];
201 /* Reset potential sums */
205 /* Start inner kernel loop */
206 for(jidx=j_index_start; jidx<j_index_end; jidx++)
208 /* Get j neighbor index, and coordinate index */
210 j_coord_offset = DIM*jnr;
212 /* load j atom coordinates */
213 jx0 = x[j_coord_offset+DIM*0+XX];
214 jy0 = x[j_coord_offset+DIM*0+YY];
215 jz0 = x[j_coord_offset+DIM*0+ZZ];
216 jx1 = x[j_coord_offset+DIM*1+XX];
217 jy1 = x[j_coord_offset+DIM*1+YY];
218 jz1 = x[j_coord_offset+DIM*1+ZZ];
219 jx2 = x[j_coord_offset+DIM*2+XX];
220 jy2 = x[j_coord_offset+DIM*2+YY];
221 jz2 = x[j_coord_offset+DIM*2+ZZ];
222 jx3 = x[j_coord_offset+DIM*3+XX];
223 jy3 = x[j_coord_offset+DIM*3+YY];
224 jz3 = x[j_coord_offset+DIM*3+ZZ];
226 /* Calculate displacement vector */
258 /* Calculate squared distance and things based on it */
259 rsq00 = dx00*dx00+dy00*dy00+dz00*dz00;
260 rsq11 = dx11*dx11+dy11*dy11+dz11*dz11;
261 rsq12 = dx12*dx12+dy12*dy12+dz12*dz12;
262 rsq13 = dx13*dx13+dy13*dy13+dz13*dz13;
263 rsq21 = dx21*dx21+dy21*dy21+dz21*dz21;
264 rsq22 = dx22*dx22+dy22*dy22+dz22*dz22;
265 rsq23 = dx23*dx23+dy23*dy23+dz23*dz23;
266 rsq31 = dx31*dx31+dy31*dy31+dz31*dz31;
267 rsq32 = dx32*dx32+dy32*dy32+dz32*dz32;
268 rsq33 = dx33*dx33+dy33*dy33+dz33*dz33;
270 rinv00 = gmx_invsqrt(rsq00);
271 rinv11 = gmx_invsqrt(rsq11);
272 rinv12 = gmx_invsqrt(rsq12);
273 rinv13 = gmx_invsqrt(rsq13);
274 rinv21 = gmx_invsqrt(rsq21);
275 rinv22 = gmx_invsqrt(rsq22);
276 rinv23 = gmx_invsqrt(rsq23);
277 rinv31 = gmx_invsqrt(rsq31);
278 rinv32 = gmx_invsqrt(rsq32);
279 rinv33 = gmx_invsqrt(rsq33);
281 rinvsq00 = rinv00*rinv00;
282 rinvsq11 = rinv11*rinv11;
283 rinvsq12 = rinv12*rinv12;
284 rinvsq13 = rinv13*rinv13;
285 rinvsq21 = rinv21*rinv21;
286 rinvsq22 = rinv22*rinv22;
287 rinvsq23 = rinv23*rinv23;
288 rinvsq31 = rinv31*rinv31;
289 rinvsq32 = rinv32*rinv32;
290 rinvsq33 = rinv33*rinv33;
292 /**************************
293 * CALCULATE INTERACTIONS *
294 **************************/
301 /* BUCKINGHAM DISPERSION/REPULSION */
302 rinvsix = rinvsq00*rinvsq00*rinvsq00;
303 vvdw6 = c6_00*rinvsix;
305 vvdwexp = cexp1_00*exp(-br);
306 vvdw = (vvdwexp-cexp1_00*exp(-cexp2_00*rvdw)) - (vvdw6 - c6_00*sh_vdw_invrcut6)*(1.0/6.0);
307 fvdw = (br*vvdwexp-vvdw6)*rinvsq00;
309 /* Update potential sums from outer loop */
314 /* Calculate temporary vectorial force */
319 /* Update vectorial force */
323 f[j_coord_offset+DIM*0+XX] -= tx;
324 f[j_coord_offset+DIM*0+YY] -= ty;
325 f[j_coord_offset+DIM*0+ZZ] -= tz;
329 /**************************
330 * CALCULATE INTERACTIONS *
331 **************************/
336 /* REACTION-FIELD ELECTROSTATICS */
337 velec = qq11*(rinv11+krf*rsq11-crf);
338 felec = qq11*(rinv11*rinvsq11-krf2);
340 /* Update potential sums from outer loop */
345 /* Calculate temporary vectorial force */
350 /* Update vectorial force */
354 f[j_coord_offset+DIM*1+XX] -= tx;
355 f[j_coord_offset+DIM*1+YY] -= ty;
356 f[j_coord_offset+DIM*1+ZZ] -= tz;
360 /**************************
361 * CALCULATE INTERACTIONS *
362 **************************/
367 /* REACTION-FIELD ELECTROSTATICS */
368 velec = qq12*(rinv12+krf*rsq12-crf);
369 felec = qq12*(rinv12*rinvsq12-krf2);
371 /* Update potential sums from outer loop */
376 /* Calculate temporary vectorial force */
381 /* Update vectorial force */
385 f[j_coord_offset+DIM*2+XX] -= tx;
386 f[j_coord_offset+DIM*2+YY] -= ty;
387 f[j_coord_offset+DIM*2+ZZ] -= tz;
391 /**************************
392 * CALCULATE INTERACTIONS *
393 **************************/
398 /* REACTION-FIELD ELECTROSTATICS */
399 velec = qq13*(rinv13+krf*rsq13-crf);
400 felec = qq13*(rinv13*rinvsq13-krf2);
402 /* Update potential sums from outer loop */
407 /* Calculate temporary vectorial force */
412 /* Update vectorial force */
416 f[j_coord_offset+DIM*3+XX] -= tx;
417 f[j_coord_offset+DIM*3+YY] -= ty;
418 f[j_coord_offset+DIM*3+ZZ] -= tz;
422 /**************************
423 * CALCULATE INTERACTIONS *
424 **************************/
429 /* REACTION-FIELD ELECTROSTATICS */
430 velec = qq21*(rinv21+krf*rsq21-crf);
431 felec = qq21*(rinv21*rinvsq21-krf2);
433 /* Update potential sums from outer loop */
438 /* Calculate temporary vectorial force */
443 /* Update vectorial force */
447 f[j_coord_offset+DIM*1+XX] -= tx;
448 f[j_coord_offset+DIM*1+YY] -= ty;
449 f[j_coord_offset+DIM*1+ZZ] -= tz;
453 /**************************
454 * CALCULATE INTERACTIONS *
455 **************************/
460 /* REACTION-FIELD ELECTROSTATICS */
461 velec = qq22*(rinv22+krf*rsq22-crf);
462 felec = qq22*(rinv22*rinvsq22-krf2);
464 /* Update potential sums from outer loop */
469 /* Calculate temporary vectorial force */
474 /* Update vectorial force */
478 f[j_coord_offset+DIM*2+XX] -= tx;
479 f[j_coord_offset+DIM*2+YY] -= ty;
480 f[j_coord_offset+DIM*2+ZZ] -= tz;
484 /**************************
485 * CALCULATE INTERACTIONS *
486 **************************/
491 /* REACTION-FIELD ELECTROSTATICS */
492 velec = qq23*(rinv23+krf*rsq23-crf);
493 felec = qq23*(rinv23*rinvsq23-krf2);
495 /* Update potential sums from outer loop */
500 /* Calculate temporary vectorial force */
505 /* Update vectorial force */
509 f[j_coord_offset+DIM*3+XX] -= tx;
510 f[j_coord_offset+DIM*3+YY] -= ty;
511 f[j_coord_offset+DIM*3+ZZ] -= tz;
515 /**************************
516 * CALCULATE INTERACTIONS *
517 **************************/
522 /* REACTION-FIELD ELECTROSTATICS */
523 velec = qq31*(rinv31+krf*rsq31-crf);
524 felec = qq31*(rinv31*rinvsq31-krf2);
526 /* Update potential sums from outer loop */
531 /* Calculate temporary vectorial force */
536 /* Update vectorial force */
540 f[j_coord_offset+DIM*1+XX] -= tx;
541 f[j_coord_offset+DIM*1+YY] -= ty;
542 f[j_coord_offset+DIM*1+ZZ] -= tz;
546 /**************************
547 * CALCULATE INTERACTIONS *
548 **************************/
553 /* REACTION-FIELD ELECTROSTATICS */
554 velec = qq32*(rinv32+krf*rsq32-crf);
555 felec = qq32*(rinv32*rinvsq32-krf2);
557 /* Update potential sums from outer loop */
562 /* Calculate temporary vectorial force */
567 /* Update vectorial force */
571 f[j_coord_offset+DIM*2+XX] -= tx;
572 f[j_coord_offset+DIM*2+YY] -= ty;
573 f[j_coord_offset+DIM*2+ZZ] -= tz;
577 /**************************
578 * CALCULATE INTERACTIONS *
579 **************************/
584 /* REACTION-FIELD ELECTROSTATICS */
585 velec = qq33*(rinv33+krf*rsq33-crf);
586 felec = qq33*(rinv33*rinvsq33-krf2);
588 /* Update potential sums from outer loop */
593 /* Calculate temporary vectorial force */
598 /* Update vectorial force */
602 f[j_coord_offset+DIM*3+XX] -= tx;
603 f[j_coord_offset+DIM*3+YY] -= ty;
604 f[j_coord_offset+DIM*3+ZZ] -= tz;
608 /* Inner loop uses 371 flops */
610 /* End of innermost loop */
613 f[i_coord_offset+DIM*0+XX] += fix0;
614 f[i_coord_offset+DIM*0+YY] += fiy0;
615 f[i_coord_offset+DIM*0+ZZ] += fiz0;
619 f[i_coord_offset+DIM*1+XX] += fix1;
620 f[i_coord_offset+DIM*1+YY] += fiy1;
621 f[i_coord_offset+DIM*1+ZZ] += fiz1;
625 f[i_coord_offset+DIM*2+XX] += fix2;
626 f[i_coord_offset+DIM*2+YY] += fiy2;
627 f[i_coord_offset+DIM*2+ZZ] += fiz2;
631 f[i_coord_offset+DIM*3+XX] += fix3;
632 f[i_coord_offset+DIM*3+YY] += fiy3;
633 f[i_coord_offset+DIM*3+ZZ] += fiz3;
637 fshift[i_shift_offset+XX] += tx;
638 fshift[i_shift_offset+YY] += ty;
639 fshift[i_shift_offset+ZZ] += tz;
642 /* Update potential energies */
643 kernel_data->energygrp_elec[ggid] += velecsum;
644 kernel_data->energygrp_vdw[ggid] += vvdwsum;
646 /* Increment number of inner iterations */
647 inneriter += j_index_end - j_index_start;
649 /* Outer loop uses 41 flops */
652 /* Increment number of outer iterations */
655 /* Update outer/inner flops */
657 inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4W4_VF,outeriter*41 + inneriter*371);
660 * Gromacs nonbonded kernel: nb_kernel_ElecRFCut_VdwBhamSh_GeomW4W4_F_c
661 * Electrostatics interaction: ReactionField
662 * VdW interaction: Buckingham
663 * Geometry: Water4-Water4
664 * Calculate force/pot: Force
667 nb_kernel_ElecRFCut_VdwBhamSh_GeomW4W4_F_c
668 (t_nblist * gmx_restrict nlist,
669 rvec * gmx_restrict xx,
670 rvec * gmx_restrict ff,
671 t_forcerec * gmx_restrict fr,
672 t_mdatoms * gmx_restrict mdatoms,
673 nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
674 t_nrnb * gmx_restrict nrnb)
676 int i_shift_offset,i_coord_offset,j_coord_offset;
677 int j_index_start,j_index_end;
678 int nri,inr,ggid,iidx,jidx,jnr,outeriter,inneriter;
679 real shX,shY,shZ,tx,ty,tz,fscal,rcutoff,rcutoff2;
680 int *iinr,*jindex,*jjnr,*shiftidx,*gid;
681 real *shiftvec,*fshift,*x,*f;
683 real ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
685 real ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
687 real ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
689 real ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
691 real jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
693 real jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
695 real jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
697 real jx3,jy3,jz3,fjx3,fjy3,fjz3,jq3,isaj3;
698 real dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00,cexp1_00,cexp2_00;
699 real dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11,cexp1_11,cexp2_11;
700 real dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12,cexp1_12,cexp2_12;
701 real dx13,dy13,dz13,rsq13,rinv13,rinvsq13,r13,qq13,c6_13,c12_13,cexp1_13,cexp2_13;
702 real dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21,cexp1_21,cexp2_21;
703 real dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22,cexp1_22,cexp2_22;
704 real dx23,dy23,dz23,rsq23,rinv23,rinvsq23,r23,qq23,c6_23,c12_23,cexp1_23,cexp2_23;
705 real dx31,dy31,dz31,rsq31,rinv31,rinvsq31,r31,qq31,c6_31,c12_31,cexp1_31,cexp2_31;
706 real dx32,dy32,dz32,rsq32,rinv32,rinvsq32,r32,qq32,c6_32,c12_32,cexp1_32,cexp2_32;
707 real dx33,dy33,dz33,rsq33,rinv33,rinvsq33,r33,qq33,c6_33,c12_33,cexp1_33,cexp2_33;
708 real velec,felec,velecsum,facel,crf,krf,krf2;
711 real rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,br,vvdwexp,sh_vdw_invrcut6;
720 jindex = nlist->jindex;
722 shiftidx = nlist->shift;
724 shiftvec = fr->shift_vec[0];
725 fshift = fr->fshift[0];
727 charge = mdatoms->chargeA;
731 nvdwtype = fr->ntype;
733 vdwtype = mdatoms->typeA;
735 /* Setup water-specific parameters */
736 inr = nlist->iinr[0];
737 iq1 = facel*charge[inr+1];
738 iq2 = facel*charge[inr+2];
739 iq3 = facel*charge[inr+3];
740 vdwioffset0 = 3*nvdwtype*vdwtype[inr+0];
745 vdwjidx0 = 3*vdwtype[inr+0];
746 c6_00 = vdwparam[vdwioffset0+vdwjidx0];
747 cexp1_00 = vdwparam[vdwioffset0+vdwjidx0+1];
748 cexp2_00 = vdwparam[vdwioffset0+vdwjidx0+2];
759 /* When we use explicit cutoffs the value must be identical for elec and VdW, so use elec as an arbitrary choice */
760 rcutoff = fr->rcoulomb;
761 rcutoff2 = rcutoff*rcutoff;
763 sh_vdw_invrcut6 = fr->ic->sh_invrc6;
769 /* Start outer loop over neighborlists */
770 for(iidx=0; iidx<nri; iidx++)
772 /* Load shift vector for this list */
773 i_shift_offset = DIM*shiftidx[iidx];
774 shX = shiftvec[i_shift_offset+XX];
775 shY = shiftvec[i_shift_offset+YY];
776 shZ = shiftvec[i_shift_offset+ZZ];
778 /* Load limits for loop over neighbors */
779 j_index_start = jindex[iidx];
780 j_index_end = jindex[iidx+1];
782 /* Get outer coordinate index */
784 i_coord_offset = DIM*inr;
786 /* Load i particle coords and add shift vector */
787 ix0 = shX + x[i_coord_offset+DIM*0+XX];
788 iy0 = shY + x[i_coord_offset+DIM*0+YY];
789 iz0 = shZ + x[i_coord_offset+DIM*0+ZZ];
790 ix1 = shX + x[i_coord_offset+DIM*1+XX];
791 iy1 = shY + x[i_coord_offset+DIM*1+YY];
792 iz1 = shZ + x[i_coord_offset+DIM*1+ZZ];
793 ix2 = shX + x[i_coord_offset+DIM*2+XX];
794 iy2 = shY + x[i_coord_offset+DIM*2+YY];
795 iz2 = shZ + x[i_coord_offset+DIM*2+ZZ];
796 ix3 = shX + x[i_coord_offset+DIM*3+XX];
797 iy3 = shY + x[i_coord_offset+DIM*3+YY];
798 iz3 = shZ + x[i_coord_offset+DIM*3+ZZ];
813 /* Start inner kernel loop */
814 for(jidx=j_index_start; jidx<j_index_end; jidx++)
816 /* Get j neighbor index, and coordinate index */
818 j_coord_offset = DIM*jnr;
820 /* load j atom coordinates */
821 jx0 = x[j_coord_offset+DIM*0+XX];
822 jy0 = x[j_coord_offset+DIM*0+YY];
823 jz0 = x[j_coord_offset+DIM*0+ZZ];
824 jx1 = x[j_coord_offset+DIM*1+XX];
825 jy1 = x[j_coord_offset+DIM*1+YY];
826 jz1 = x[j_coord_offset+DIM*1+ZZ];
827 jx2 = x[j_coord_offset+DIM*2+XX];
828 jy2 = x[j_coord_offset+DIM*2+YY];
829 jz2 = x[j_coord_offset+DIM*2+ZZ];
830 jx3 = x[j_coord_offset+DIM*3+XX];
831 jy3 = x[j_coord_offset+DIM*3+YY];
832 jz3 = x[j_coord_offset+DIM*3+ZZ];
834 /* Calculate displacement vector */
866 /* Calculate squared distance and things based on it */
867 rsq00 = dx00*dx00+dy00*dy00+dz00*dz00;
868 rsq11 = dx11*dx11+dy11*dy11+dz11*dz11;
869 rsq12 = dx12*dx12+dy12*dy12+dz12*dz12;
870 rsq13 = dx13*dx13+dy13*dy13+dz13*dz13;
871 rsq21 = dx21*dx21+dy21*dy21+dz21*dz21;
872 rsq22 = dx22*dx22+dy22*dy22+dz22*dz22;
873 rsq23 = dx23*dx23+dy23*dy23+dz23*dz23;
874 rsq31 = dx31*dx31+dy31*dy31+dz31*dz31;
875 rsq32 = dx32*dx32+dy32*dy32+dz32*dz32;
876 rsq33 = dx33*dx33+dy33*dy33+dz33*dz33;
878 rinv00 = gmx_invsqrt(rsq00);
879 rinv11 = gmx_invsqrt(rsq11);
880 rinv12 = gmx_invsqrt(rsq12);
881 rinv13 = gmx_invsqrt(rsq13);
882 rinv21 = gmx_invsqrt(rsq21);
883 rinv22 = gmx_invsqrt(rsq22);
884 rinv23 = gmx_invsqrt(rsq23);
885 rinv31 = gmx_invsqrt(rsq31);
886 rinv32 = gmx_invsqrt(rsq32);
887 rinv33 = gmx_invsqrt(rsq33);
889 rinvsq00 = rinv00*rinv00;
890 rinvsq11 = rinv11*rinv11;
891 rinvsq12 = rinv12*rinv12;
892 rinvsq13 = rinv13*rinv13;
893 rinvsq21 = rinv21*rinv21;
894 rinvsq22 = rinv22*rinv22;
895 rinvsq23 = rinv23*rinv23;
896 rinvsq31 = rinv31*rinv31;
897 rinvsq32 = rinv32*rinv32;
898 rinvsq33 = rinv33*rinv33;
900 /**************************
901 * CALCULATE INTERACTIONS *
902 **************************/
909 /* BUCKINGHAM DISPERSION/REPULSION */
910 rinvsix = rinvsq00*rinvsq00*rinvsq00;
911 vvdw6 = c6_00*rinvsix;
913 vvdwexp = cexp1_00*exp(-br);
914 fvdw = (br*vvdwexp-vvdw6)*rinvsq00;
918 /* Calculate temporary vectorial force */
923 /* Update vectorial force */
927 f[j_coord_offset+DIM*0+XX] -= tx;
928 f[j_coord_offset+DIM*0+YY] -= ty;
929 f[j_coord_offset+DIM*0+ZZ] -= tz;
933 /**************************
934 * CALCULATE INTERACTIONS *
935 **************************/
940 /* REACTION-FIELD ELECTROSTATICS */
941 felec = qq11*(rinv11*rinvsq11-krf2);
945 /* Calculate temporary vectorial force */
950 /* Update vectorial force */
954 f[j_coord_offset+DIM*1+XX] -= tx;
955 f[j_coord_offset+DIM*1+YY] -= ty;
956 f[j_coord_offset+DIM*1+ZZ] -= tz;
960 /**************************
961 * CALCULATE INTERACTIONS *
962 **************************/
967 /* REACTION-FIELD ELECTROSTATICS */
968 felec = qq12*(rinv12*rinvsq12-krf2);
972 /* Calculate temporary vectorial force */
977 /* Update vectorial force */
981 f[j_coord_offset+DIM*2+XX] -= tx;
982 f[j_coord_offset+DIM*2+YY] -= ty;
983 f[j_coord_offset+DIM*2+ZZ] -= tz;
987 /**************************
988 * CALCULATE INTERACTIONS *
989 **************************/
994 /* REACTION-FIELD ELECTROSTATICS */
995 felec = qq13*(rinv13*rinvsq13-krf2);
999 /* Calculate temporary vectorial force */
1004 /* Update vectorial force */
1008 f[j_coord_offset+DIM*3+XX] -= tx;
1009 f[j_coord_offset+DIM*3+YY] -= ty;
1010 f[j_coord_offset+DIM*3+ZZ] -= tz;
1014 /**************************
1015 * CALCULATE INTERACTIONS *
1016 **************************/
1021 /* REACTION-FIELD ELECTROSTATICS */
1022 felec = qq21*(rinv21*rinvsq21-krf2);
1026 /* Calculate temporary vectorial force */
1031 /* Update vectorial force */
1035 f[j_coord_offset+DIM*1+XX] -= tx;
1036 f[j_coord_offset+DIM*1+YY] -= ty;
1037 f[j_coord_offset+DIM*1+ZZ] -= tz;
1041 /**************************
1042 * CALCULATE INTERACTIONS *
1043 **************************/
1048 /* REACTION-FIELD ELECTROSTATICS */
1049 felec = qq22*(rinv22*rinvsq22-krf2);
1053 /* Calculate temporary vectorial force */
1058 /* Update vectorial force */
1062 f[j_coord_offset+DIM*2+XX] -= tx;
1063 f[j_coord_offset+DIM*2+YY] -= ty;
1064 f[j_coord_offset+DIM*2+ZZ] -= tz;
1068 /**************************
1069 * CALCULATE INTERACTIONS *
1070 **************************/
1075 /* REACTION-FIELD ELECTROSTATICS */
1076 felec = qq23*(rinv23*rinvsq23-krf2);
1080 /* Calculate temporary vectorial force */
1085 /* Update vectorial force */
1089 f[j_coord_offset+DIM*3+XX] -= tx;
1090 f[j_coord_offset+DIM*3+YY] -= ty;
1091 f[j_coord_offset+DIM*3+ZZ] -= tz;
1095 /**************************
1096 * CALCULATE INTERACTIONS *
1097 **************************/
1102 /* REACTION-FIELD ELECTROSTATICS */
1103 felec = qq31*(rinv31*rinvsq31-krf2);
1107 /* Calculate temporary vectorial force */
1112 /* Update vectorial force */
1116 f[j_coord_offset+DIM*1+XX] -= tx;
1117 f[j_coord_offset+DIM*1+YY] -= ty;
1118 f[j_coord_offset+DIM*1+ZZ] -= tz;
1122 /**************************
1123 * CALCULATE INTERACTIONS *
1124 **************************/
1129 /* REACTION-FIELD ELECTROSTATICS */
1130 felec = qq32*(rinv32*rinvsq32-krf2);
1134 /* Calculate temporary vectorial force */
1139 /* Update vectorial force */
1143 f[j_coord_offset+DIM*2+XX] -= tx;
1144 f[j_coord_offset+DIM*2+YY] -= ty;
1145 f[j_coord_offset+DIM*2+ZZ] -= tz;
1149 /**************************
1150 * CALCULATE INTERACTIONS *
1151 **************************/
1156 /* REACTION-FIELD ELECTROSTATICS */
1157 felec = qq33*(rinv33*rinvsq33-krf2);
1161 /* Calculate temporary vectorial force */
1166 /* Update vectorial force */
1170 f[j_coord_offset+DIM*3+XX] -= tx;
1171 f[j_coord_offset+DIM*3+YY] -= ty;
1172 f[j_coord_offset+DIM*3+ZZ] -= tz;
1176 /* Inner loop uses 292 flops */
1178 /* End of innermost loop */
1181 f[i_coord_offset+DIM*0+XX] += fix0;
1182 f[i_coord_offset+DIM*0+YY] += fiy0;
1183 f[i_coord_offset+DIM*0+ZZ] += fiz0;
1187 f[i_coord_offset+DIM*1+XX] += fix1;
1188 f[i_coord_offset+DIM*1+YY] += fiy1;
1189 f[i_coord_offset+DIM*1+ZZ] += fiz1;
1193 f[i_coord_offset+DIM*2+XX] += fix2;
1194 f[i_coord_offset+DIM*2+YY] += fiy2;
1195 f[i_coord_offset+DIM*2+ZZ] += fiz2;
1199 f[i_coord_offset+DIM*3+XX] += fix3;
1200 f[i_coord_offset+DIM*3+YY] += fiy3;
1201 f[i_coord_offset+DIM*3+ZZ] += fiz3;
1205 fshift[i_shift_offset+XX] += tx;
1206 fshift[i_shift_offset+YY] += ty;
1207 fshift[i_shift_offset+ZZ] += tz;
1209 /* Increment number of inner iterations */
1210 inneriter += j_index_end - j_index_start;
1212 /* Outer loop uses 39 flops */
1215 /* Increment number of outer iterations */
1218 /* Update outer/inner flops */
1220 inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4W4_F,outeriter*39 + inneriter*292);