2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013, 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 "types/simple.h"
50 * Gromacs nonbonded kernel: nb_kernel_ElecCoul_VdwLJ_GeomW4W4_VF_c
51 * Electrostatics interaction: Coulomb
52 * VdW interaction: LennardJones
53 * Geometry: Water4-Water4
54 * Calculate force/pot: PotentialAndForce
57 nb_kernel_ElecCoul_VdwLJ_GeomW4W4_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 ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
81 real jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
83 real jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
85 real jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
87 real jx3,jy3,jz3,fjx3,fjy3,fjz3,jq3,isaj3;
88 real dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00,cexp1_00,cexp2_00;
89 real dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11,cexp1_11,cexp2_11;
90 real dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12,cexp1_12,cexp2_12;
91 real dx13,dy13,dz13,rsq13,rinv13,rinvsq13,r13,qq13,c6_13,c12_13,cexp1_13,cexp2_13;
92 real dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21,cexp1_21,cexp2_21;
93 real dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22,cexp1_22,cexp2_22;
94 real dx23,dy23,dz23,rsq23,rinv23,rinvsq23,r23,qq23,c6_23,c12_23,cexp1_23,cexp2_23;
95 real dx31,dy31,dz31,rsq31,rinv31,rinvsq31,r31,qq31,c6_31,c12_31,cexp1_31,cexp2_31;
96 real dx32,dy32,dz32,rsq32,rinv32,rinvsq32,r32,qq32,c6_32,c12_32,cexp1_32,cexp2_32;
97 real dx33,dy33,dz33,rsq33,rinv33,rinvsq33,r33,qq33,c6_33,c12_33,cexp1_33,cexp2_33;
98 real velec,felec,velecsum,facel,crf,krf,krf2;
101 real rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,br,vvdwexp,sh_vdw_invrcut6;
110 jindex = nlist->jindex;
112 shiftidx = nlist->shift;
114 shiftvec = fr->shift_vec[0];
115 fshift = fr->fshift[0];
117 charge = mdatoms->chargeA;
118 nvdwtype = fr->ntype;
120 vdwtype = mdatoms->typeA;
122 /* Setup water-specific parameters */
123 inr = nlist->iinr[0];
124 iq1 = facel*charge[inr+1];
125 iq2 = facel*charge[inr+2];
126 iq3 = facel*charge[inr+3];
127 vdwioffset0 = 2*nvdwtype*vdwtype[inr+0];
132 vdwjidx0 = 2*vdwtype[inr+0];
133 c6_00 = vdwparam[vdwioffset0+vdwjidx0];
134 c12_00 = vdwparam[vdwioffset0+vdwjidx0+1];
148 /* Start outer loop over neighborlists */
149 for(iidx=0; iidx<nri; iidx++)
151 /* Load shift vector for this list */
152 i_shift_offset = DIM*shiftidx[iidx];
153 shX = shiftvec[i_shift_offset+XX];
154 shY = shiftvec[i_shift_offset+YY];
155 shZ = shiftvec[i_shift_offset+ZZ];
157 /* Load limits for loop over neighbors */
158 j_index_start = jindex[iidx];
159 j_index_end = jindex[iidx+1];
161 /* Get outer coordinate index */
163 i_coord_offset = DIM*inr;
165 /* Load i particle coords and add shift vector */
166 ix0 = shX + x[i_coord_offset+DIM*0+XX];
167 iy0 = shY + x[i_coord_offset+DIM*0+YY];
168 iz0 = shZ + x[i_coord_offset+DIM*0+ZZ];
169 ix1 = shX + x[i_coord_offset+DIM*1+XX];
170 iy1 = shY + x[i_coord_offset+DIM*1+YY];
171 iz1 = shZ + x[i_coord_offset+DIM*1+ZZ];
172 ix2 = shX + x[i_coord_offset+DIM*2+XX];
173 iy2 = shY + x[i_coord_offset+DIM*2+YY];
174 iz2 = shZ + x[i_coord_offset+DIM*2+ZZ];
175 ix3 = shX + x[i_coord_offset+DIM*3+XX];
176 iy3 = shY + x[i_coord_offset+DIM*3+YY];
177 iz3 = shZ + x[i_coord_offset+DIM*3+ZZ];
192 /* Reset potential sums */
196 /* Start inner kernel loop */
197 for(jidx=j_index_start; jidx<j_index_end; jidx++)
199 /* Get j neighbor index, and coordinate index */
201 j_coord_offset = DIM*jnr;
203 /* load j atom coordinates */
204 jx0 = x[j_coord_offset+DIM*0+XX];
205 jy0 = x[j_coord_offset+DIM*0+YY];
206 jz0 = x[j_coord_offset+DIM*0+ZZ];
207 jx1 = x[j_coord_offset+DIM*1+XX];
208 jy1 = x[j_coord_offset+DIM*1+YY];
209 jz1 = x[j_coord_offset+DIM*1+ZZ];
210 jx2 = x[j_coord_offset+DIM*2+XX];
211 jy2 = x[j_coord_offset+DIM*2+YY];
212 jz2 = x[j_coord_offset+DIM*2+ZZ];
213 jx3 = x[j_coord_offset+DIM*3+XX];
214 jy3 = x[j_coord_offset+DIM*3+YY];
215 jz3 = x[j_coord_offset+DIM*3+ZZ];
217 /* Calculate displacement vector */
249 /* Calculate squared distance and things based on it */
250 rsq00 = dx00*dx00+dy00*dy00+dz00*dz00;
251 rsq11 = dx11*dx11+dy11*dy11+dz11*dz11;
252 rsq12 = dx12*dx12+dy12*dy12+dz12*dz12;
253 rsq13 = dx13*dx13+dy13*dy13+dz13*dz13;
254 rsq21 = dx21*dx21+dy21*dy21+dz21*dz21;
255 rsq22 = dx22*dx22+dy22*dy22+dz22*dz22;
256 rsq23 = dx23*dx23+dy23*dy23+dz23*dz23;
257 rsq31 = dx31*dx31+dy31*dy31+dz31*dz31;
258 rsq32 = dx32*dx32+dy32*dy32+dz32*dz32;
259 rsq33 = dx33*dx33+dy33*dy33+dz33*dz33;
261 rinv11 = gmx_invsqrt(rsq11);
262 rinv12 = gmx_invsqrt(rsq12);
263 rinv13 = gmx_invsqrt(rsq13);
264 rinv21 = gmx_invsqrt(rsq21);
265 rinv22 = gmx_invsqrt(rsq22);
266 rinv23 = gmx_invsqrt(rsq23);
267 rinv31 = gmx_invsqrt(rsq31);
268 rinv32 = gmx_invsqrt(rsq32);
269 rinv33 = gmx_invsqrt(rsq33);
271 rinvsq00 = 1.0/rsq00;
272 rinvsq11 = rinv11*rinv11;
273 rinvsq12 = rinv12*rinv12;
274 rinvsq13 = rinv13*rinv13;
275 rinvsq21 = rinv21*rinv21;
276 rinvsq22 = rinv22*rinv22;
277 rinvsq23 = rinv23*rinv23;
278 rinvsq31 = rinv31*rinv31;
279 rinvsq32 = rinv32*rinv32;
280 rinvsq33 = rinv33*rinv33;
282 /**************************
283 * CALCULATE INTERACTIONS *
284 **************************/
286 /* LENNARD-JONES DISPERSION/REPULSION */
288 rinvsix = rinvsq00*rinvsq00*rinvsq00;
289 vvdw6 = c6_00*rinvsix;
290 vvdw12 = c12_00*rinvsix*rinvsix;
291 vvdw = vvdw12*(1.0/12.0) - vvdw6*(1.0/6.0);
292 fvdw = (vvdw12-vvdw6)*rinvsq00;
294 /* Update potential sums from outer loop */
299 /* Calculate temporary vectorial force */
304 /* Update vectorial force */
308 f[j_coord_offset+DIM*0+XX] -= tx;
309 f[j_coord_offset+DIM*0+YY] -= ty;
310 f[j_coord_offset+DIM*0+ZZ] -= tz;
312 /**************************
313 * CALCULATE INTERACTIONS *
314 **************************/
316 /* COULOMB ELECTROSTATICS */
318 felec = velec*rinvsq11;
320 /* Update potential sums from outer loop */
325 /* Calculate temporary vectorial force */
330 /* Update vectorial force */
334 f[j_coord_offset+DIM*1+XX] -= tx;
335 f[j_coord_offset+DIM*1+YY] -= ty;
336 f[j_coord_offset+DIM*1+ZZ] -= tz;
338 /**************************
339 * CALCULATE INTERACTIONS *
340 **************************/
342 /* COULOMB ELECTROSTATICS */
344 felec = velec*rinvsq12;
346 /* Update potential sums from outer loop */
351 /* Calculate temporary vectorial force */
356 /* Update vectorial force */
360 f[j_coord_offset+DIM*2+XX] -= tx;
361 f[j_coord_offset+DIM*2+YY] -= ty;
362 f[j_coord_offset+DIM*2+ZZ] -= tz;
364 /**************************
365 * CALCULATE INTERACTIONS *
366 **************************/
368 /* COULOMB ELECTROSTATICS */
370 felec = velec*rinvsq13;
372 /* Update potential sums from outer loop */
377 /* Calculate temporary vectorial force */
382 /* Update vectorial force */
386 f[j_coord_offset+DIM*3+XX] -= tx;
387 f[j_coord_offset+DIM*3+YY] -= ty;
388 f[j_coord_offset+DIM*3+ZZ] -= tz;
390 /**************************
391 * CALCULATE INTERACTIONS *
392 **************************/
394 /* COULOMB ELECTROSTATICS */
396 felec = velec*rinvsq21;
398 /* Update potential sums from outer loop */
403 /* Calculate temporary vectorial force */
408 /* Update vectorial force */
412 f[j_coord_offset+DIM*1+XX] -= tx;
413 f[j_coord_offset+DIM*1+YY] -= ty;
414 f[j_coord_offset+DIM*1+ZZ] -= tz;
416 /**************************
417 * CALCULATE INTERACTIONS *
418 **************************/
420 /* COULOMB ELECTROSTATICS */
422 felec = velec*rinvsq22;
424 /* Update potential sums from outer loop */
429 /* Calculate temporary vectorial force */
434 /* Update vectorial force */
438 f[j_coord_offset+DIM*2+XX] -= tx;
439 f[j_coord_offset+DIM*2+YY] -= ty;
440 f[j_coord_offset+DIM*2+ZZ] -= tz;
442 /**************************
443 * CALCULATE INTERACTIONS *
444 **************************/
446 /* COULOMB ELECTROSTATICS */
448 felec = velec*rinvsq23;
450 /* Update potential sums from outer loop */
455 /* Calculate temporary vectorial force */
460 /* Update vectorial force */
464 f[j_coord_offset+DIM*3+XX] -= tx;
465 f[j_coord_offset+DIM*3+YY] -= ty;
466 f[j_coord_offset+DIM*3+ZZ] -= tz;
468 /**************************
469 * CALCULATE INTERACTIONS *
470 **************************/
472 /* COULOMB ELECTROSTATICS */
474 felec = velec*rinvsq31;
476 /* Update potential sums from outer loop */
481 /* Calculate temporary vectorial force */
486 /* Update vectorial force */
490 f[j_coord_offset+DIM*1+XX] -= tx;
491 f[j_coord_offset+DIM*1+YY] -= ty;
492 f[j_coord_offset+DIM*1+ZZ] -= tz;
494 /**************************
495 * CALCULATE INTERACTIONS *
496 **************************/
498 /* COULOMB ELECTROSTATICS */
500 felec = velec*rinvsq32;
502 /* Update potential sums from outer loop */
507 /* Calculate temporary vectorial force */
512 /* Update vectorial force */
516 f[j_coord_offset+DIM*2+XX] -= tx;
517 f[j_coord_offset+DIM*2+YY] -= ty;
518 f[j_coord_offset+DIM*2+ZZ] -= tz;
520 /**************************
521 * CALCULATE INTERACTIONS *
522 **************************/
524 /* COULOMB ELECTROSTATICS */
526 felec = velec*rinvsq33;
528 /* Update potential sums from outer loop */
533 /* Calculate temporary vectorial force */
538 /* Update vectorial force */
542 f[j_coord_offset+DIM*3+XX] -= tx;
543 f[j_coord_offset+DIM*3+YY] -= ty;
544 f[j_coord_offset+DIM*3+ZZ] -= tz;
546 /* Inner loop uses 275 flops */
548 /* End of innermost loop */
551 f[i_coord_offset+DIM*0+XX] += fix0;
552 f[i_coord_offset+DIM*0+YY] += fiy0;
553 f[i_coord_offset+DIM*0+ZZ] += fiz0;
557 f[i_coord_offset+DIM*1+XX] += fix1;
558 f[i_coord_offset+DIM*1+YY] += fiy1;
559 f[i_coord_offset+DIM*1+ZZ] += fiz1;
563 f[i_coord_offset+DIM*2+XX] += fix2;
564 f[i_coord_offset+DIM*2+YY] += fiy2;
565 f[i_coord_offset+DIM*2+ZZ] += fiz2;
569 f[i_coord_offset+DIM*3+XX] += fix3;
570 f[i_coord_offset+DIM*3+YY] += fiy3;
571 f[i_coord_offset+DIM*3+ZZ] += fiz3;
575 fshift[i_shift_offset+XX] += tx;
576 fshift[i_shift_offset+YY] += ty;
577 fshift[i_shift_offset+ZZ] += tz;
580 /* Update potential energies */
581 kernel_data->energygrp_elec[ggid] += velecsum;
582 kernel_data->energygrp_vdw[ggid] += vvdwsum;
584 /* Increment number of inner iterations */
585 inneriter += j_index_end - j_index_start;
587 /* Outer loop uses 41 flops */
590 /* Increment number of outer iterations */
593 /* Update outer/inner flops */
595 inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4W4_VF,outeriter*41 + inneriter*275);
598 * Gromacs nonbonded kernel: nb_kernel_ElecCoul_VdwLJ_GeomW4W4_F_c
599 * Electrostatics interaction: Coulomb
600 * VdW interaction: LennardJones
601 * Geometry: Water4-Water4
602 * Calculate force/pot: Force
605 nb_kernel_ElecCoul_VdwLJ_GeomW4W4_F_c
606 (t_nblist * gmx_restrict nlist,
607 rvec * gmx_restrict xx,
608 rvec * gmx_restrict ff,
609 t_forcerec * gmx_restrict fr,
610 t_mdatoms * gmx_restrict mdatoms,
611 nb_kernel_data_t gmx_unused * gmx_restrict kernel_data,
612 t_nrnb * gmx_restrict nrnb)
614 int i_shift_offset,i_coord_offset,j_coord_offset;
615 int j_index_start,j_index_end;
616 int nri,inr,ggid,iidx,jidx,jnr,outeriter,inneriter;
617 real shX,shY,shZ,tx,ty,tz,fscal,rcutoff,rcutoff2;
618 int *iinr,*jindex,*jjnr,*shiftidx,*gid;
619 real *shiftvec,*fshift,*x,*f;
621 real ix0,iy0,iz0,fix0,fiy0,fiz0,iq0,isai0;
623 real ix1,iy1,iz1,fix1,fiy1,fiz1,iq1,isai1;
625 real ix2,iy2,iz2,fix2,fiy2,fiz2,iq2,isai2;
627 real ix3,iy3,iz3,fix3,fiy3,fiz3,iq3,isai3;
629 real jx0,jy0,jz0,fjx0,fjy0,fjz0,jq0,isaj0;
631 real jx1,jy1,jz1,fjx1,fjy1,fjz1,jq1,isaj1;
633 real jx2,jy2,jz2,fjx2,fjy2,fjz2,jq2,isaj2;
635 real jx3,jy3,jz3,fjx3,fjy3,fjz3,jq3,isaj3;
636 real dx00,dy00,dz00,rsq00,rinv00,rinvsq00,r00,qq00,c6_00,c12_00,cexp1_00,cexp2_00;
637 real dx11,dy11,dz11,rsq11,rinv11,rinvsq11,r11,qq11,c6_11,c12_11,cexp1_11,cexp2_11;
638 real dx12,dy12,dz12,rsq12,rinv12,rinvsq12,r12,qq12,c6_12,c12_12,cexp1_12,cexp2_12;
639 real dx13,dy13,dz13,rsq13,rinv13,rinvsq13,r13,qq13,c6_13,c12_13,cexp1_13,cexp2_13;
640 real dx21,dy21,dz21,rsq21,rinv21,rinvsq21,r21,qq21,c6_21,c12_21,cexp1_21,cexp2_21;
641 real dx22,dy22,dz22,rsq22,rinv22,rinvsq22,r22,qq22,c6_22,c12_22,cexp1_22,cexp2_22;
642 real dx23,dy23,dz23,rsq23,rinv23,rinvsq23,r23,qq23,c6_23,c12_23,cexp1_23,cexp2_23;
643 real dx31,dy31,dz31,rsq31,rinv31,rinvsq31,r31,qq31,c6_31,c12_31,cexp1_31,cexp2_31;
644 real dx32,dy32,dz32,rsq32,rinv32,rinvsq32,r32,qq32,c6_32,c12_32,cexp1_32,cexp2_32;
645 real dx33,dy33,dz33,rsq33,rinv33,rinvsq33,r33,qq33,c6_33,c12_33,cexp1_33,cexp2_33;
646 real velec,felec,velecsum,facel,crf,krf,krf2;
649 real rinvsix,rvdw,vvdw,vvdw6,vvdw12,fvdw,fvdw6,fvdw12,vvdwsum,br,vvdwexp,sh_vdw_invrcut6;
658 jindex = nlist->jindex;
660 shiftidx = nlist->shift;
662 shiftvec = fr->shift_vec[0];
663 fshift = fr->fshift[0];
665 charge = mdatoms->chargeA;
666 nvdwtype = fr->ntype;
668 vdwtype = mdatoms->typeA;
670 /* Setup water-specific parameters */
671 inr = nlist->iinr[0];
672 iq1 = facel*charge[inr+1];
673 iq2 = facel*charge[inr+2];
674 iq3 = facel*charge[inr+3];
675 vdwioffset0 = 2*nvdwtype*vdwtype[inr+0];
680 vdwjidx0 = 2*vdwtype[inr+0];
681 c6_00 = vdwparam[vdwioffset0+vdwjidx0];
682 c12_00 = vdwparam[vdwioffset0+vdwjidx0+1];
696 /* Start outer loop over neighborlists */
697 for(iidx=0; iidx<nri; iidx++)
699 /* Load shift vector for this list */
700 i_shift_offset = DIM*shiftidx[iidx];
701 shX = shiftvec[i_shift_offset+XX];
702 shY = shiftvec[i_shift_offset+YY];
703 shZ = shiftvec[i_shift_offset+ZZ];
705 /* Load limits for loop over neighbors */
706 j_index_start = jindex[iidx];
707 j_index_end = jindex[iidx+1];
709 /* Get outer coordinate index */
711 i_coord_offset = DIM*inr;
713 /* Load i particle coords and add shift vector */
714 ix0 = shX + x[i_coord_offset+DIM*0+XX];
715 iy0 = shY + x[i_coord_offset+DIM*0+YY];
716 iz0 = shZ + x[i_coord_offset+DIM*0+ZZ];
717 ix1 = shX + x[i_coord_offset+DIM*1+XX];
718 iy1 = shY + x[i_coord_offset+DIM*1+YY];
719 iz1 = shZ + x[i_coord_offset+DIM*1+ZZ];
720 ix2 = shX + x[i_coord_offset+DIM*2+XX];
721 iy2 = shY + x[i_coord_offset+DIM*2+YY];
722 iz2 = shZ + x[i_coord_offset+DIM*2+ZZ];
723 ix3 = shX + x[i_coord_offset+DIM*3+XX];
724 iy3 = shY + x[i_coord_offset+DIM*3+YY];
725 iz3 = shZ + x[i_coord_offset+DIM*3+ZZ];
740 /* Start inner kernel loop */
741 for(jidx=j_index_start; jidx<j_index_end; jidx++)
743 /* Get j neighbor index, and coordinate index */
745 j_coord_offset = DIM*jnr;
747 /* load j atom coordinates */
748 jx0 = x[j_coord_offset+DIM*0+XX];
749 jy0 = x[j_coord_offset+DIM*0+YY];
750 jz0 = x[j_coord_offset+DIM*0+ZZ];
751 jx1 = x[j_coord_offset+DIM*1+XX];
752 jy1 = x[j_coord_offset+DIM*1+YY];
753 jz1 = x[j_coord_offset+DIM*1+ZZ];
754 jx2 = x[j_coord_offset+DIM*2+XX];
755 jy2 = x[j_coord_offset+DIM*2+YY];
756 jz2 = x[j_coord_offset+DIM*2+ZZ];
757 jx3 = x[j_coord_offset+DIM*3+XX];
758 jy3 = x[j_coord_offset+DIM*3+YY];
759 jz3 = x[j_coord_offset+DIM*3+ZZ];
761 /* Calculate displacement vector */
793 /* Calculate squared distance and things based on it */
794 rsq00 = dx00*dx00+dy00*dy00+dz00*dz00;
795 rsq11 = dx11*dx11+dy11*dy11+dz11*dz11;
796 rsq12 = dx12*dx12+dy12*dy12+dz12*dz12;
797 rsq13 = dx13*dx13+dy13*dy13+dz13*dz13;
798 rsq21 = dx21*dx21+dy21*dy21+dz21*dz21;
799 rsq22 = dx22*dx22+dy22*dy22+dz22*dz22;
800 rsq23 = dx23*dx23+dy23*dy23+dz23*dz23;
801 rsq31 = dx31*dx31+dy31*dy31+dz31*dz31;
802 rsq32 = dx32*dx32+dy32*dy32+dz32*dz32;
803 rsq33 = dx33*dx33+dy33*dy33+dz33*dz33;
805 rinv11 = gmx_invsqrt(rsq11);
806 rinv12 = gmx_invsqrt(rsq12);
807 rinv13 = gmx_invsqrt(rsq13);
808 rinv21 = gmx_invsqrt(rsq21);
809 rinv22 = gmx_invsqrt(rsq22);
810 rinv23 = gmx_invsqrt(rsq23);
811 rinv31 = gmx_invsqrt(rsq31);
812 rinv32 = gmx_invsqrt(rsq32);
813 rinv33 = gmx_invsqrt(rsq33);
815 rinvsq00 = 1.0/rsq00;
816 rinvsq11 = rinv11*rinv11;
817 rinvsq12 = rinv12*rinv12;
818 rinvsq13 = rinv13*rinv13;
819 rinvsq21 = rinv21*rinv21;
820 rinvsq22 = rinv22*rinv22;
821 rinvsq23 = rinv23*rinv23;
822 rinvsq31 = rinv31*rinv31;
823 rinvsq32 = rinv32*rinv32;
824 rinvsq33 = rinv33*rinv33;
826 /**************************
827 * CALCULATE INTERACTIONS *
828 **************************/
830 /* LENNARD-JONES DISPERSION/REPULSION */
832 rinvsix = rinvsq00*rinvsq00*rinvsq00;
833 fvdw = (c12_00*rinvsix-c6_00)*rinvsix*rinvsq00;
837 /* Calculate temporary vectorial force */
842 /* Update vectorial force */
846 f[j_coord_offset+DIM*0+XX] -= tx;
847 f[j_coord_offset+DIM*0+YY] -= ty;
848 f[j_coord_offset+DIM*0+ZZ] -= tz;
850 /**************************
851 * CALCULATE INTERACTIONS *
852 **************************/
854 /* COULOMB ELECTROSTATICS */
856 felec = velec*rinvsq11;
860 /* Calculate temporary vectorial force */
865 /* Update vectorial force */
869 f[j_coord_offset+DIM*1+XX] -= tx;
870 f[j_coord_offset+DIM*1+YY] -= ty;
871 f[j_coord_offset+DIM*1+ZZ] -= tz;
873 /**************************
874 * CALCULATE INTERACTIONS *
875 **************************/
877 /* COULOMB ELECTROSTATICS */
879 felec = velec*rinvsq12;
883 /* Calculate temporary vectorial force */
888 /* Update vectorial force */
892 f[j_coord_offset+DIM*2+XX] -= tx;
893 f[j_coord_offset+DIM*2+YY] -= ty;
894 f[j_coord_offset+DIM*2+ZZ] -= tz;
896 /**************************
897 * CALCULATE INTERACTIONS *
898 **************************/
900 /* COULOMB ELECTROSTATICS */
902 felec = velec*rinvsq13;
906 /* Calculate temporary vectorial force */
911 /* Update vectorial force */
915 f[j_coord_offset+DIM*3+XX] -= tx;
916 f[j_coord_offset+DIM*3+YY] -= ty;
917 f[j_coord_offset+DIM*3+ZZ] -= tz;
919 /**************************
920 * CALCULATE INTERACTIONS *
921 **************************/
923 /* COULOMB ELECTROSTATICS */
925 felec = velec*rinvsq21;
929 /* Calculate temporary vectorial force */
934 /* Update vectorial force */
938 f[j_coord_offset+DIM*1+XX] -= tx;
939 f[j_coord_offset+DIM*1+YY] -= ty;
940 f[j_coord_offset+DIM*1+ZZ] -= tz;
942 /**************************
943 * CALCULATE INTERACTIONS *
944 **************************/
946 /* COULOMB ELECTROSTATICS */
948 felec = velec*rinvsq22;
952 /* Calculate temporary vectorial force */
957 /* Update vectorial force */
961 f[j_coord_offset+DIM*2+XX] -= tx;
962 f[j_coord_offset+DIM*2+YY] -= ty;
963 f[j_coord_offset+DIM*2+ZZ] -= tz;
965 /**************************
966 * CALCULATE INTERACTIONS *
967 **************************/
969 /* COULOMB ELECTROSTATICS */
971 felec = velec*rinvsq23;
975 /* Calculate temporary vectorial force */
980 /* Update vectorial force */
984 f[j_coord_offset+DIM*3+XX] -= tx;
985 f[j_coord_offset+DIM*3+YY] -= ty;
986 f[j_coord_offset+DIM*3+ZZ] -= tz;
988 /**************************
989 * CALCULATE INTERACTIONS *
990 **************************/
992 /* COULOMB ELECTROSTATICS */
994 felec = velec*rinvsq31;
998 /* Calculate temporary vectorial force */
1003 /* Update vectorial force */
1007 f[j_coord_offset+DIM*1+XX] -= tx;
1008 f[j_coord_offset+DIM*1+YY] -= ty;
1009 f[j_coord_offset+DIM*1+ZZ] -= tz;
1011 /**************************
1012 * CALCULATE INTERACTIONS *
1013 **************************/
1015 /* COULOMB ELECTROSTATICS */
1016 velec = qq32*rinv32;
1017 felec = velec*rinvsq32;
1021 /* Calculate temporary vectorial force */
1026 /* Update vectorial force */
1030 f[j_coord_offset+DIM*2+XX] -= tx;
1031 f[j_coord_offset+DIM*2+YY] -= ty;
1032 f[j_coord_offset+DIM*2+ZZ] -= tz;
1034 /**************************
1035 * CALCULATE INTERACTIONS *
1036 **************************/
1038 /* COULOMB ELECTROSTATICS */
1039 velec = qq33*rinv33;
1040 felec = velec*rinvsq33;
1044 /* Calculate temporary vectorial force */
1049 /* Update vectorial force */
1053 f[j_coord_offset+DIM*3+XX] -= tx;
1054 f[j_coord_offset+DIM*3+YY] -= ty;
1055 f[j_coord_offset+DIM*3+ZZ] -= tz;
1057 /* Inner loop uses 261 flops */
1059 /* End of innermost loop */
1062 f[i_coord_offset+DIM*0+XX] += fix0;
1063 f[i_coord_offset+DIM*0+YY] += fiy0;
1064 f[i_coord_offset+DIM*0+ZZ] += fiz0;
1068 f[i_coord_offset+DIM*1+XX] += fix1;
1069 f[i_coord_offset+DIM*1+YY] += fiy1;
1070 f[i_coord_offset+DIM*1+ZZ] += fiz1;
1074 f[i_coord_offset+DIM*2+XX] += fix2;
1075 f[i_coord_offset+DIM*2+YY] += fiy2;
1076 f[i_coord_offset+DIM*2+ZZ] += fiz2;
1080 f[i_coord_offset+DIM*3+XX] += fix3;
1081 f[i_coord_offset+DIM*3+YY] += fiy3;
1082 f[i_coord_offset+DIM*3+ZZ] += fiz3;
1086 fshift[i_shift_offset+XX] += tx;
1087 fshift[i_shift_offset+YY] += ty;
1088 fshift[i_shift_offset+ZZ] += tz;
1090 /* Increment number of inner iterations */
1091 inneriter += j_index_end - j_index_start;
1093 /* Outer loop uses 39 flops */
1096 /* Increment number of outer iterations */
1099 /* Update outer/inner flops */
1101 inc_nrnb(nrnb,eNR_NBKERNEL_ELEC_VDW_W4W4_F,outeriter*39 + inneriter*261);