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 "types/simple.h"
44 #include "gromacs/math/vec.h"
48 * Gromacs nonbonded kernel: nb_kernel_ElecCSTab_VdwCSTab_GeomW4W4_VF_c
49 * Electrostatics interaction: CubicSplineTable
50 * VdW interaction: CubicSplineTable
51 * Geometry: Water4-Water4
52 * Calculate force/pot: PotentialAndForce
55 nb_kernel_ElecCSTab_VdwCSTab_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
;
103 real rt
,vfeps
,vftabscale
,Y
,F
,Geps
,Heps2
,Fp
,VV
,FF
;
111 jindex
= nlist
->jindex
;
113 shiftidx
= nlist
->shift
;
115 shiftvec
= fr
->shift_vec
[0];
116 fshift
= fr
->fshift
[0];
118 charge
= mdatoms
->chargeA
;
119 nvdwtype
= fr
->ntype
;
121 vdwtype
= mdatoms
->typeA
;
123 vftab
= kernel_data
->table_elec_vdw
->data
;
124 vftabscale
= kernel_data
->table_elec_vdw
->scale
;
126 /* Setup water-specific parameters */
127 inr
= nlist
->iinr
[0];
128 iq1
= facel
*charge
[inr
+1];
129 iq2
= facel
*charge
[inr
+2];
130 iq3
= facel
*charge
[inr
+3];
131 vdwioffset0
= 2*nvdwtype
*vdwtype
[inr
+0];
136 vdwjidx0
= 2*vdwtype
[inr
+0];
137 c6_00
= vdwparam
[vdwioffset0
+vdwjidx0
];
138 c12_00
= vdwparam
[vdwioffset0
+vdwjidx0
+1];
152 /* Start outer loop over neighborlists */
153 for(iidx
=0; iidx
<nri
; iidx
++)
155 /* Load shift vector for this list */
156 i_shift_offset
= DIM
*shiftidx
[iidx
];
157 shX
= shiftvec
[i_shift_offset
+XX
];
158 shY
= shiftvec
[i_shift_offset
+YY
];
159 shZ
= shiftvec
[i_shift_offset
+ZZ
];
161 /* Load limits for loop over neighbors */
162 j_index_start
= jindex
[iidx
];
163 j_index_end
= jindex
[iidx
+1];
165 /* Get outer coordinate index */
167 i_coord_offset
= DIM
*inr
;
169 /* Load i particle coords and add shift vector */
170 ix0
= shX
+ x
[i_coord_offset
+DIM
*0+XX
];
171 iy0
= shY
+ x
[i_coord_offset
+DIM
*0+YY
];
172 iz0
= shZ
+ x
[i_coord_offset
+DIM
*0+ZZ
];
173 ix1
= shX
+ x
[i_coord_offset
+DIM
*1+XX
];
174 iy1
= shY
+ x
[i_coord_offset
+DIM
*1+YY
];
175 iz1
= shZ
+ x
[i_coord_offset
+DIM
*1+ZZ
];
176 ix2
= shX
+ x
[i_coord_offset
+DIM
*2+XX
];
177 iy2
= shY
+ x
[i_coord_offset
+DIM
*2+YY
];
178 iz2
= shZ
+ x
[i_coord_offset
+DIM
*2+ZZ
];
179 ix3
= shX
+ x
[i_coord_offset
+DIM
*3+XX
];
180 iy3
= shY
+ x
[i_coord_offset
+DIM
*3+YY
];
181 iz3
= shZ
+ x
[i_coord_offset
+DIM
*3+ZZ
];
196 /* Reset potential sums */
200 /* Start inner kernel loop */
201 for(jidx
=j_index_start
; jidx
<j_index_end
; jidx
++)
203 /* Get j neighbor index, and coordinate index */
205 j_coord_offset
= DIM
*jnr
;
207 /* load j atom coordinates */
208 jx0
= x
[j_coord_offset
+DIM
*0+XX
];
209 jy0
= x
[j_coord_offset
+DIM
*0+YY
];
210 jz0
= x
[j_coord_offset
+DIM
*0+ZZ
];
211 jx1
= x
[j_coord_offset
+DIM
*1+XX
];
212 jy1
= x
[j_coord_offset
+DIM
*1+YY
];
213 jz1
= x
[j_coord_offset
+DIM
*1+ZZ
];
214 jx2
= x
[j_coord_offset
+DIM
*2+XX
];
215 jy2
= x
[j_coord_offset
+DIM
*2+YY
];
216 jz2
= x
[j_coord_offset
+DIM
*2+ZZ
];
217 jx3
= x
[j_coord_offset
+DIM
*3+XX
];
218 jy3
= x
[j_coord_offset
+DIM
*3+YY
];
219 jz3
= x
[j_coord_offset
+DIM
*3+ZZ
];
221 /* Calculate displacement vector */
253 /* Calculate squared distance and things based on it */
254 rsq00
= dx00
*dx00
+dy00
*dy00
+dz00
*dz00
;
255 rsq11
= dx11
*dx11
+dy11
*dy11
+dz11
*dz11
;
256 rsq12
= dx12
*dx12
+dy12
*dy12
+dz12
*dz12
;
257 rsq13
= dx13
*dx13
+dy13
*dy13
+dz13
*dz13
;
258 rsq21
= dx21
*dx21
+dy21
*dy21
+dz21
*dz21
;
259 rsq22
= dx22
*dx22
+dy22
*dy22
+dz22
*dz22
;
260 rsq23
= dx23
*dx23
+dy23
*dy23
+dz23
*dz23
;
261 rsq31
= dx31
*dx31
+dy31
*dy31
+dz31
*dz31
;
262 rsq32
= dx32
*dx32
+dy32
*dy32
+dz32
*dz32
;
263 rsq33
= dx33
*dx33
+dy33
*dy33
+dz33
*dz33
;
265 rinv00
= gmx_invsqrt(rsq00
);
266 rinv11
= gmx_invsqrt(rsq11
);
267 rinv12
= gmx_invsqrt(rsq12
);
268 rinv13
= gmx_invsqrt(rsq13
);
269 rinv21
= gmx_invsqrt(rsq21
);
270 rinv22
= gmx_invsqrt(rsq22
);
271 rinv23
= gmx_invsqrt(rsq23
);
272 rinv31
= gmx_invsqrt(rsq31
);
273 rinv32
= gmx_invsqrt(rsq32
);
274 rinv33
= gmx_invsqrt(rsq33
);
276 /**************************
277 * CALCULATE INTERACTIONS *
278 **************************/
282 /* Calculate table index by multiplying r with table scale and truncate to integer */
288 /* CUBIC SPLINE TABLE DISPERSION */
292 Geps
= vfeps
*vftab
[vfitab
+2];
293 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
297 FF
= Fp
+Geps
+2.0*Heps2
;
300 /* CUBIC SPLINE TABLE REPULSION */
303 Geps
= vfeps
*vftab
[vfitab
+6];
304 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+7];
308 FF
= Fp
+Geps
+2.0*Heps2
;
311 fvdw
= -(fvdw6
+fvdw12
)*vftabscale
*rinv00
;
313 /* Update potential sums from outer loop */
318 /* Calculate temporary vectorial force */
323 /* Update vectorial force */
327 f
[j_coord_offset
+DIM
*0+XX
] -= tx
;
328 f
[j_coord_offset
+DIM
*0+YY
] -= ty
;
329 f
[j_coord_offset
+DIM
*0+ZZ
] -= tz
;
331 /**************************
332 * CALCULATE INTERACTIONS *
333 **************************/
337 /* Calculate table index by multiplying r with table scale and truncate to integer */
343 /* CUBIC SPLINE TABLE ELECTROSTATICS */
346 Geps
= vfeps
*vftab
[vfitab
+2];
347 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
351 FF
= Fp
+Geps
+2.0*Heps2
;
352 felec
= -qq11
*FF
*vftabscale
*rinv11
;
354 /* Update potential sums from outer loop */
359 /* Calculate temporary vectorial force */
364 /* Update vectorial force */
368 f
[j_coord_offset
+DIM
*1+XX
] -= tx
;
369 f
[j_coord_offset
+DIM
*1+YY
] -= ty
;
370 f
[j_coord_offset
+DIM
*1+ZZ
] -= tz
;
372 /**************************
373 * CALCULATE INTERACTIONS *
374 **************************/
378 /* Calculate table index by multiplying r with table scale and truncate to integer */
384 /* CUBIC SPLINE TABLE ELECTROSTATICS */
387 Geps
= vfeps
*vftab
[vfitab
+2];
388 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
392 FF
= Fp
+Geps
+2.0*Heps2
;
393 felec
= -qq12
*FF
*vftabscale
*rinv12
;
395 /* Update potential sums from outer loop */
400 /* Calculate temporary vectorial force */
405 /* Update vectorial force */
409 f
[j_coord_offset
+DIM
*2+XX
] -= tx
;
410 f
[j_coord_offset
+DIM
*2+YY
] -= ty
;
411 f
[j_coord_offset
+DIM
*2+ZZ
] -= tz
;
413 /**************************
414 * CALCULATE INTERACTIONS *
415 **************************/
419 /* Calculate table index by multiplying r with table scale and truncate to integer */
425 /* CUBIC SPLINE TABLE ELECTROSTATICS */
428 Geps
= vfeps
*vftab
[vfitab
+2];
429 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
433 FF
= Fp
+Geps
+2.0*Heps2
;
434 felec
= -qq13
*FF
*vftabscale
*rinv13
;
436 /* Update potential sums from outer loop */
441 /* Calculate temporary vectorial force */
446 /* Update vectorial force */
450 f
[j_coord_offset
+DIM
*3+XX
] -= tx
;
451 f
[j_coord_offset
+DIM
*3+YY
] -= ty
;
452 f
[j_coord_offset
+DIM
*3+ZZ
] -= tz
;
454 /**************************
455 * CALCULATE INTERACTIONS *
456 **************************/
460 /* Calculate table index by multiplying r with table scale and truncate to integer */
466 /* CUBIC SPLINE TABLE ELECTROSTATICS */
469 Geps
= vfeps
*vftab
[vfitab
+2];
470 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
474 FF
= Fp
+Geps
+2.0*Heps2
;
475 felec
= -qq21
*FF
*vftabscale
*rinv21
;
477 /* Update potential sums from outer loop */
482 /* Calculate temporary vectorial force */
487 /* Update vectorial force */
491 f
[j_coord_offset
+DIM
*1+XX
] -= tx
;
492 f
[j_coord_offset
+DIM
*1+YY
] -= ty
;
493 f
[j_coord_offset
+DIM
*1+ZZ
] -= tz
;
495 /**************************
496 * CALCULATE INTERACTIONS *
497 **************************/
501 /* Calculate table index by multiplying r with table scale and truncate to integer */
507 /* CUBIC SPLINE TABLE ELECTROSTATICS */
510 Geps
= vfeps
*vftab
[vfitab
+2];
511 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
515 FF
= Fp
+Geps
+2.0*Heps2
;
516 felec
= -qq22
*FF
*vftabscale
*rinv22
;
518 /* Update potential sums from outer loop */
523 /* Calculate temporary vectorial force */
528 /* Update vectorial force */
532 f
[j_coord_offset
+DIM
*2+XX
] -= tx
;
533 f
[j_coord_offset
+DIM
*2+YY
] -= ty
;
534 f
[j_coord_offset
+DIM
*2+ZZ
] -= tz
;
536 /**************************
537 * CALCULATE INTERACTIONS *
538 **************************/
542 /* Calculate table index by multiplying r with table scale and truncate to integer */
548 /* CUBIC SPLINE TABLE ELECTROSTATICS */
551 Geps
= vfeps
*vftab
[vfitab
+2];
552 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
556 FF
= Fp
+Geps
+2.0*Heps2
;
557 felec
= -qq23
*FF
*vftabscale
*rinv23
;
559 /* Update potential sums from outer loop */
564 /* Calculate temporary vectorial force */
569 /* Update vectorial force */
573 f
[j_coord_offset
+DIM
*3+XX
] -= tx
;
574 f
[j_coord_offset
+DIM
*3+YY
] -= ty
;
575 f
[j_coord_offset
+DIM
*3+ZZ
] -= tz
;
577 /**************************
578 * CALCULATE INTERACTIONS *
579 **************************/
583 /* Calculate table index by multiplying r with table scale and truncate to integer */
589 /* CUBIC SPLINE TABLE ELECTROSTATICS */
592 Geps
= vfeps
*vftab
[vfitab
+2];
593 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
597 FF
= Fp
+Geps
+2.0*Heps2
;
598 felec
= -qq31
*FF
*vftabscale
*rinv31
;
600 /* Update potential sums from outer loop */
605 /* Calculate temporary vectorial force */
610 /* Update vectorial force */
614 f
[j_coord_offset
+DIM
*1+XX
] -= tx
;
615 f
[j_coord_offset
+DIM
*1+YY
] -= ty
;
616 f
[j_coord_offset
+DIM
*1+ZZ
] -= tz
;
618 /**************************
619 * CALCULATE INTERACTIONS *
620 **************************/
624 /* Calculate table index by multiplying r with table scale and truncate to integer */
630 /* CUBIC SPLINE TABLE ELECTROSTATICS */
633 Geps
= vfeps
*vftab
[vfitab
+2];
634 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
638 FF
= Fp
+Geps
+2.0*Heps2
;
639 felec
= -qq32
*FF
*vftabscale
*rinv32
;
641 /* Update potential sums from outer loop */
646 /* Calculate temporary vectorial force */
651 /* Update vectorial force */
655 f
[j_coord_offset
+DIM
*2+XX
] -= tx
;
656 f
[j_coord_offset
+DIM
*2+YY
] -= ty
;
657 f
[j_coord_offset
+DIM
*2+ZZ
] -= tz
;
659 /**************************
660 * CALCULATE INTERACTIONS *
661 **************************/
665 /* Calculate table index by multiplying r with table scale and truncate to integer */
671 /* CUBIC SPLINE TABLE ELECTROSTATICS */
674 Geps
= vfeps
*vftab
[vfitab
+2];
675 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
679 FF
= Fp
+Geps
+2.0*Heps2
;
680 felec
= -qq33
*FF
*vftabscale
*rinv33
;
682 /* Update potential sums from outer loop */
687 /* Calculate temporary vectorial force */
692 /* Update vectorial force */
696 f
[j_coord_offset
+DIM
*3+XX
] -= tx
;
697 f
[j_coord_offset
+DIM
*3+YY
] -= ty
;
698 f
[j_coord_offset
+DIM
*3+ZZ
] -= tz
;
700 /* Inner loop uses 424 flops */
702 /* End of innermost loop */
705 f
[i_coord_offset
+DIM
*0+XX
] += fix0
;
706 f
[i_coord_offset
+DIM
*0+YY
] += fiy0
;
707 f
[i_coord_offset
+DIM
*0+ZZ
] += fiz0
;
711 f
[i_coord_offset
+DIM
*1+XX
] += fix1
;
712 f
[i_coord_offset
+DIM
*1+YY
] += fiy1
;
713 f
[i_coord_offset
+DIM
*1+ZZ
] += fiz1
;
717 f
[i_coord_offset
+DIM
*2+XX
] += fix2
;
718 f
[i_coord_offset
+DIM
*2+YY
] += fiy2
;
719 f
[i_coord_offset
+DIM
*2+ZZ
] += fiz2
;
723 f
[i_coord_offset
+DIM
*3+XX
] += fix3
;
724 f
[i_coord_offset
+DIM
*3+YY
] += fiy3
;
725 f
[i_coord_offset
+DIM
*3+ZZ
] += fiz3
;
729 fshift
[i_shift_offset
+XX
] += tx
;
730 fshift
[i_shift_offset
+YY
] += ty
;
731 fshift
[i_shift_offset
+ZZ
] += tz
;
734 /* Update potential energies */
735 kernel_data
->energygrp_elec
[ggid
] += velecsum
;
736 kernel_data
->energygrp_vdw
[ggid
] += vvdwsum
;
738 /* Increment number of inner iterations */
739 inneriter
+= j_index_end
- j_index_start
;
741 /* Outer loop uses 41 flops */
744 /* Increment number of outer iterations */
747 /* Update outer/inner flops */
749 inc_nrnb(nrnb
,eNR_NBKERNEL_ELEC_VDW_W4W4_VF
,outeriter
*41 + inneriter
*424);
752 * Gromacs nonbonded kernel: nb_kernel_ElecCSTab_VdwCSTab_GeomW4W4_F_c
753 * Electrostatics interaction: CubicSplineTable
754 * VdW interaction: CubicSplineTable
755 * Geometry: Water4-Water4
756 * Calculate force/pot: Force
759 nb_kernel_ElecCSTab_VdwCSTab_GeomW4W4_F_c
760 (t_nblist
* gmx_restrict nlist
,
761 rvec
* gmx_restrict xx
,
762 rvec
* gmx_restrict ff
,
763 t_forcerec
* gmx_restrict fr
,
764 t_mdatoms
* gmx_restrict mdatoms
,
765 nb_kernel_data_t gmx_unused
* gmx_restrict kernel_data
,
766 t_nrnb
* gmx_restrict nrnb
)
768 int i_shift_offset
,i_coord_offset
,j_coord_offset
;
769 int j_index_start
,j_index_end
;
770 int nri
,inr
,ggid
,iidx
,jidx
,jnr
,outeriter
,inneriter
;
771 real shX
,shY
,shZ
,tx
,ty
,tz
,fscal
,rcutoff
,rcutoff2
;
772 int *iinr
,*jindex
,*jjnr
,*shiftidx
,*gid
;
773 real
*shiftvec
,*fshift
,*x
,*f
;
775 real ix0
,iy0
,iz0
,fix0
,fiy0
,fiz0
,iq0
,isai0
;
777 real ix1
,iy1
,iz1
,fix1
,fiy1
,fiz1
,iq1
,isai1
;
779 real ix2
,iy2
,iz2
,fix2
,fiy2
,fiz2
,iq2
,isai2
;
781 real ix3
,iy3
,iz3
,fix3
,fiy3
,fiz3
,iq3
,isai3
;
783 real jx0
,jy0
,jz0
,fjx0
,fjy0
,fjz0
,jq0
,isaj0
;
785 real jx1
,jy1
,jz1
,fjx1
,fjy1
,fjz1
,jq1
,isaj1
;
787 real jx2
,jy2
,jz2
,fjx2
,fjy2
,fjz2
,jq2
,isaj2
;
789 real jx3
,jy3
,jz3
,fjx3
,fjy3
,fjz3
,jq3
,isaj3
;
790 real dx00
,dy00
,dz00
,rsq00
,rinv00
,rinvsq00
,r00
,qq00
,c6_00
,c12_00
,cexp1_00
,cexp2_00
;
791 real dx11
,dy11
,dz11
,rsq11
,rinv11
,rinvsq11
,r11
,qq11
,c6_11
,c12_11
,cexp1_11
,cexp2_11
;
792 real dx12
,dy12
,dz12
,rsq12
,rinv12
,rinvsq12
,r12
,qq12
,c6_12
,c12_12
,cexp1_12
,cexp2_12
;
793 real dx13
,dy13
,dz13
,rsq13
,rinv13
,rinvsq13
,r13
,qq13
,c6_13
,c12_13
,cexp1_13
,cexp2_13
;
794 real dx21
,dy21
,dz21
,rsq21
,rinv21
,rinvsq21
,r21
,qq21
,c6_21
,c12_21
,cexp1_21
,cexp2_21
;
795 real dx22
,dy22
,dz22
,rsq22
,rinv22
,rinvsq22
,r22
,qq22
,c6_22
,c12_22
,cexp1_22
,cexp2_22
;
796 real dx23
,dy23
,dz23
,rsq23
,rinv23
,rinvsq23
,r23
,qq23
,c6_23
,c12_23
,cexp1_23
,cexp2_23
;
797 real dx31
,dy31
,dz31
,rsq31
,rinv31
,rinvsq31
,r31
,qq31
,c6_31
,c12_31
,cexp1_31
,cexp2_31
;
798 real dx32
,dy32
,dz32
,rsq32
,rinv32
,rinvsq32
,r32
,qq32
,c6_32
,c12_32
,cexp1_32
,cexp2_32
;
799 real dx33
,dy33
,dz33
,rsq33
,rinv33
,rinvsq33
,r33
,qq33
,c6_33
,c12_33
,cexp1_33
,cexp2_33
;
800 real velec
,felec
,velecsum
,facel
,crf
,krf
,krf2
;
803 real rinvsix
,rvdw
,vvdw
,vvdw6
,vvdw12
,fvdw
,fvdw6
,fvdw12
,vvdwsum
,br
,vvdwexp
,sh_vdw_invrcut6
;
807 real rt
,vfeps
,vftabscale
,Y
,F
,Geps
,Heps2
,Fp
,VV
,FF
;
815 jindex
= nlist
->jindex
;
817 shiftidx
= nlist
->shift
;
819 shiftvec
= fr
->shift_vec
[0];
820 fshift
= fr
->fshift
[0];
822 charge
= mdatoms
->chargeA
;
823 nvdwtype
= fr
->ntype
;
825 vdwtype
= mdatoms
->typeA
;
827 vftab
= kernel_data
->table_elec_vdw
->data
;
828 vftabscale
= kernel_data
->table_elec_vdw
->scale
;
830 /* Setup water-specific parameters */
831 inr
= nlist
->iinr
[0];
832 iq1
= facel
*charge
[inr
+1];
833 iq2
= facel
*charge
[inr
+2];
834 iq3
= facel
*charge
[inr
+3];
835 vdwioffset0
= 2*nvdwtype
*vdwtype
[inr
+0];
840 vdwjidx0
= 2*vdwtype
[inr
+0];
841 c6_00
= vdwparam
[vdwioffset0
+vdwjidx0
];
842 c12_00
= vdwparam
[vdwioffset0
+vdwjidx0
+1];
856 /* Start outer loop over neighborlists */
857 for(iidx
=0; iidx
<nri
; iidx
++)
859 /* Load shift vector for this list */
860 i_shift_offset
= DIM
*shiftidx
[iidx
];
861 shX
= shiftvec
[i_shift_offset
+XX
];
862 shY
= shiftvec
[i_shift_offset
+YY
];
863 shZ
= shiftvec
[i_shift_offset
+ZZ
];
865 /* Load limits for loop over neighbors */
866 j_index_start
= jindex
[iidx
];
867 j_index_end
= jindex
[iidx
+1];
869 /* Get outer coordinate index */
871 i_coord_offset
= DIM
*inr
;
873 /* Load i particle coords and add shift vector */
874 ix0
= shX
+ x
[i_coord_offset
+DIM
*0+XX
];
875 iy0
= shY
+ x
[i_coord_offset
+DIM
*0+YY
];
876 iz0
= shZ
+ x
[i_coord_offset
+DIM
*0+ZZ
];
877 ix1
= shX
+ x
[i_coord_offset
+DIM
*1+XX
];
878 iy1
= shY
+ x
[i_coord_offset
+DIM
*1+YY
];
879 iz1
= shZ
+ x
[i_coord_offset
+DIM
*1+ZZ
];
880 ix2
= shX
+ x
[i_coord_offset
+DIM
*2+XX
];
881 iy2
= shY
+ x
[i_coord_offset
+DIM
*2+YY
];
882 iz2
= shZ
+ x
[i_coord_offset
+DIM
*2+ZZ
];
883 ix3
= shX
+ x
[i_coord_offset
+DIM
*3+XX
];
884 iy3
= shY
+ x
[i_coord_offset
+DIM
*3+YY
];
885 iz3
= shZ
+ x
[i_coord_offset
+DIM
*3+ZZ
];
900 /* Start inner kernel loop */
901 for(jidx
=j_index_start
; jidx
<j_index_end
; jidx
++)
903 /* Get j neighbor index, and coordinate index */
905 j_coord_offset
= DIM
*jnr
;
907 /* load j atom coordinates */
908 jx0
= x
[j_coord_offset
+DIM
*0+XX
];
909 jy0
= x
[j_coord_offset
+DIM
*0+YY
];
910 jz0
= x
[j_coord_offset
+DIM
*0+ZZ
];
911 jx1
= x
[j_coord_offset
+DIM
*1+XX
];
912 jy1
= x
[j_coord_offset
+DIM
*1+YY
];
913 jz1
= x
[j_coord_offset
+DIM
*1+ZZ
];
914 jx2
= x
[j_coord_offset
+DIM
*2+XX
];
915 jy2
= x
[j_coord_offset
+DIM
*2+YY
];
916 jz2
= x
[j_coord_offset
+DIM
*2+ZZ
];
917 jx3
= x
[j_coord_offset
+DIM
*3+XX
];
918 jy3
= x
[j_coord_offset
+DIM
*3+YY
];
919 jz3
= x
[j_coord_offset
+DIM
*3+ZZ
];
921 /* Calculate displacement vector */
953 /* Calculate squared distance and things based on it */
954 rsq00
= dx00
*dx00
+dy00
*dy00
+dz00
*dz00
;
955 rsq11
= dx11
*dx11
+dy11
*dy11
+dz11
*dz11
;
956 rsq12
= dx12
*dx12
+dy12
*dy12
+dz12
*dz12
;
957 rsq13
= dx13
*dx13
+dy13
*dy13
+dz13
*dz13
;
958 rsq21
= dx21
*dx21
+dy21
*dy21
+dz21
*dz21
;
959 rsq22
= dx22
*dx22
+dy22
*dy22
+dz22
*dz22
;
960 rsq23
= dx23
*dx23
+dy23
*dy23
+dz23
*dz23
;
961 rsq31
= dx31
*dx31
+dy31
*dy31
+dz31
*dz31
;
962 rsq32
= dx32
*dx32
+dy32
*dy32
+dz32
*dz32
;
963 rsq33
= dx33
*dx33
+dy33
*dy33
+dz33
*dz33
;
965 rinv00
= gmx_invsqrt(rsq00
);
966 rinv11
= gmx_invsqrt(rsq11
);
967 rinv12
= gmx_invsqrt(rsq12
);
968 rinv13
= gmx_invsqrt(rsq13
);
969 rinv21
= gmx_invsqrt(rsq21
);
970 rinv22
= gmx_invsqrt(rsq22
);
971 rinv23
= gmx_invsqrt(rsq23
);
972 rinv31
= gmx_invsqrt(rsq31
);
973 rinv32
= gmx_invsqrt(rsq32
);
974 rinv33
= gmx_invsqrt(rsq33
);
976 /**************************
977 * CALCULATE INTERACTIONS *
978 **************************/
982 /* Calculate table index by multiplying r with table scale and truncate to integer */
988 /* CUBIC SPLINE TABLE DISPERSION */
991 Geps
= vfeps
*vftab
[vfitab
+2];
992 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
994 FF
= Fp
+Geps
+2.0*Heps2
;
997 /* CUBIC SPLINE TABLE REPULSION */
999 Geps
= vfeps
*vftab
[vfitab
+6];
1000 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+7];
1002 FF
= Fp
+Geps
+2.0*Heps2
;
1004 fvdw
= -(fvdw6
+fvdw12
)*vftabscale
*rinv00
;
1008 /* Calculate temporary vectorial force */
1013 /* Update vectorial force */
1017 f
[j_coord_offset
+DIM
*0+XX
] -= tx
;
1018 f
[j_coord_offset
+DIM
*0+YY
] -= ty
;
1019 f
[j_coord_offset
+DIM
*0+ZZ
] -= tz
;
1021 /**************************
1022 * CALCULATE INTERACTIONS *
1023 **************************/
1027 /* Calculate table index by multiplying r with table scale and truncate to integer */
1028 rt
= r11
*vftabscale
;
1031 vfitab
= 3*4*vfitab
;
1033 /* CUBIC SPLINE TABLE ELECTROSTATICS */
1034 F
= vftab
[vfitab
+1];
1035 Geps
= vfeps
*vftab
[vfitab
+2];
1036 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
1038 FF
= Fp
+Geps
+2.0*Heps2
;
1039 felec
= -qq11
*FF
*vftabscale
*rinv11
;
1043 /* Calculate temporary vectorial force */
1048 /* Update vectorial force */
1052 f
[j_coord_offset
+DIM
*1+XX
] -= tx
;
1053 f
[j_coord_offset
+DIM
*1+YY
] -= ty
;
1054 f
[j_coord_offset
+DIM
*1+ZZ
] -= tz
;
1056 /**************************
1057 * CALCULATE INTERACTIONS *
1058 **************************/
1062 /* Calculate table index by multiplying r with table scale and truncate to integer */
1063 rt
= r12
*vftabscale
;
1066 vfitab
= 3*4*vfitab
;
1068 /* CUBIC SPLINE TABLE ELECTROSTATICS */
1069 F
= vftab
[vfitab
+1];
1070 Geps
= vfeps
*vftab
[vfitab
+2];
1071 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
1073 FF
= Fp
+Geps
+2.0*Heps2
;
1074 felec
= -qq12
*FF
*vftabscale
*rinv12
;
1078 /* Calculate temporary vectorial force */
1083 /* Update vectorial force */
1087 f
[j_coord_offset
+DIM
*2+XX
] -= tx
;
1088 f
[j_coord_offset
+DIM
*2+YY
] -= ty
;
1089 f
[j_coord_offset
+DIM
*2+ZZ
] -= tz
;
1091 /**************************
1092 * CALCULATE INTERACTIONS *
1093 **************************/
1097 /* Calculate table index by multiplying r with table scale and truncate to integer */
1098 rt
= r13
*vftabscale
;
1101 vfitab
= 3*4*vfitab
;
1103 /* CUBIC SPLINE TABLE ELECTROSTATICS */
1104 F
= vftab
[vfitab
+1];
1105 Geps
= vfeps
*vftab
[vfitab
+2];
1106 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
1108 FF
= Fp
+Geps
+2.0*Heps2
;
1109 felec
= -qq13
*FF
*vftabscale
*rinv13
;
1113 /* Calculate temporary vectorial force */
1118 /* Update vectorial force */
1122 f
[j_coord_offset
+DIM
*3+XX
] -= tx
;
1123 f
[j_coord_offset
+DIM
*3+YY
] -= ty
;
1124 f
[j_coord_offset
+DIM
*3+ZZ
] -= tz
;
1126 /**************************
1127 * CALCULATE INTERACTIONS *
1128 **************************/
1132 /* Calculate table index by multiplying r with table scale and truncate to integer */
1133 rt
= r21
*vftabscale
;
1136 vfitab
= 3*4*vfitab
;
1138 /* CUBIC SPLINE TABLE ELECTROSTATICS */
1139 F
= vftab
[vfitab
+1];
1140 Geps
= vfeps
*vftab
[vfitab
+2];
1141 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
1143 FF
= Fp
+Geps
+2.0*Heps2
;
1144 felec
= -qq21
*FF
*vftabscale
*rinv21
;
1148 /* Calculate temporary vectorial force */
1153 /* Update vectorial force */
1157 f
[j_coord_offset
+DIM
*1+XX
] -= tx
;
1158 f
[j_coord_offset
+DIM
*1+YY
] -= ty
;
1159 f
[j_coord_offset
+DIM
*1+ZZ
] -= tz
;
1161 /**************************
1162 * CALCULATE INTERACTIONS *
1163 **************************/
1167 /* Calculate table index by multiplying r with table scale and truncate to integer */
1168 rt
= r22
*vftabscale
;
1171 vfitab
= 3*4*vfitab
;
1173 /* CUBIC SPLINE TABLE ELECTROSTATICS */
1174 F
= vftab
[vfitab
+1];
1175 Geps
= vfeps
*vftab
[vfitab
+2];
1176 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
1178 FF
= Fp
+Geps
+2.0*Heps2
;
1179 felec
= -qq22
*FF
*vftabscale
*rinv22
;
1183 /* Calculate temporary vectorial force */
1188 /* Update vectorial force */
1192 f
[j_coord_offset
+DIM
*2+XX
] -= tx
;
1193 f
[j_coord_offset
+DIM
*2+YY
] -= ty
;
1194 f
[j_coord_offset
+DIM
*2+ZZ
] -= tz
;
1196 /**************************
1197 * CALCULATE INTERACTIONS *
1198 **************************/
1202 /* Calculate table index by multiplying r with table scale and truncate to integer */
1203 rt
= r23
*vftabscale
;
1206 vfitab
= 3*4*vfitab
;
1208 /* CUBIC SPLINE TABLE ELECTROSTATICS */
1209 F
= vftab
[vfitab
+1];
1210 Geps
= vfeps
*vftab
[vfitab
+2];
1211 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
1213 FF
= Fp
+Geps
+2.0*Heps2
;
1214 felec
= -qq23
*FF
*vftabscale
*rinv23
;
1218 /* Calculate temporary vectorial force */
1223 /* Update vectorial force */
1227 f
[j_coord_offset
+DIM
*3+XX
] -= tx
;
1228 f
[j_coord_offset
+DIM
*3+YY
] -= ty
;
1229 f
[j_coord_offset
+DIM
*3+ZZ
] -= tz
;
1231 /**************************
1232 * CALCULATE INTERACTIONS *
1233 **************************/
1237 /* Calculate table index by multiplying r with table scale and truncate to integer */
1238 rt
= r31
*vftabscale
;
1241 vfitab
= 3*4*vfitab
;
1243 /* CUBIC SPLINE TABLE ELECTROSTATICS */
1244 F
= vftab
[vfitab
+1];
1245 Geps
= vfeps
*vftab
[vfitab
+2];
1246 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
1248 FF
= Fp
+Geps
+2.0*Heps2
;
1249 felec
= -qq31
*FF
*vftabscale
*rinv31
;
1253 /* Calculate temporary vectorial force */
1258 /* Update vectorial force */
1262 f
[j_coord_offset
+DIM
*1+XX
] -= tx
;
1263 f
[j_coord_offset
+DIM
*1+YY
] -= ty
;
1264 f
[j_coord_offset
+DIM
*1+ZZ
] -= tz
;
1266 /**************************
1267 * CALCULATE INTERACTIONS *
1268 **************************/
1272 /* Calculate table index by multiplying r with table scale and truncate to integer */
1273 rt
= r32
*vftabscale
;
1276 vfitab
= 3*4*vfitab
;
1278 /* CUBIC SPLINE TABLE ELECTROSTATICS */
1279 F
= vftab
[vfitab
+1];
1280 Geps
= vfeps
*vftab
[vfitab
+2];
1281 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
1283 FF
= Fp
+Geps
+2.0*Heps2
;
1284 felec
= -qq32
*FF
*vftabscale
*rinv32
;
1288 /* Calculate temporary vectorial force */
1293 /* Update vectorial force */
1297 f
[j_coord_offset
+DIM
*2+XX
] -= tx
;
1298 f
[j_coord_offset
+DIM
*2+YY
] -= ty
;
1299 f
[j_coord_offset
+DIM
*2+ZZ
] -= tz
;
1301 /**************************
1302 * CALCULATE INTERACTIONS *
1303 **************************/
1307 /* Calculate table index by multiplying r with table scale and truncate to integer */
1308 rt
= r33
*vftabscale
;
1311 vfitab
= 3*4*vfitab
;
1313 /* CUBIC SPLINE TABLE ELECTROSTATICS */
1314 F
= vftab
[vfitab
+1];
1315 Geps
= vfeps
*vftab
[vfitab
+2];
1316 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
1318 FF
= Fp
+Geps
+2.0*Heps2
;
1319 felec
= -qq33
*FF
*vftabscale
*rinv33
;
1323 /* Calculate temporary vectorial force */
1328 /* Update vectorial force */
1332 f
[j_coord_offset
+DIM
*3+XX
] -= tx
;
1333 f
[j_coord_offset
+DIM
*3+YY
] -= ty
;
1334 f
[j_coord_offset
+DIM
*3+ZZ
] -= tz
;
1336 /* Inner loop uses 380 flops */
1338 /* End of innermost loop */
1341 f
[i_coord_offset
+DIM
*0+XX
] += fix0
;
1342 f
[i_coord_offset
+DIM
*0+YY
] += fiy0
;
1343 f
[i_coord_offset
+DIM
*0+ZZ
] += fiz0
;
1347 f
[i_coord_offset
+DIM
*1+XX
] += fix1
;
1348 f
[i_coord_offset
+DIM
*1+YY
] += fiy1
;
1349 f
[i_coord_offset
+DIM
*1+ZZ
] += fiz1
;
1353 f
[i_coord_offset
+DIM
*2+XX
] += fix2
;
1354 f
[i_coord_offset
+DIM
*2+YY
] += fiy2
;
1355 f
[i_coord_offset
+DIM
*2+ZZ
] += fiz2
;
1359 f
[i_coord_offset
+DIM
*3+XX
] += fix3
;
1360 f
[i_coord_offset
+DIM
*3+YY
] += fiy3
;
1361 f
[i_coord_offset
+DIM
*3+ZZ
] += fiz3
;
1365 fshift
[i_shift_offset
+XX
] += tx
;
1366 fshift
[i_shift_offset
+YY
] += ty
;
1367 fshift
[i_shift_offset
+ZZ
] += tz
;
1369 /* Increment number of inner iterations */
1370 inneriter
+= j_index_end
- j_index_start
;
1372 /* Outer loop uses 39 flops */
1375 /* Increment number of outer iterations */
1378 /* Update outer/inner flops */
1380 inc_nrnb(nrnb
,eNR_NBKERNEL_ELEC_VDW_W4W4_F
,outeriter
*39 + inneriter
*380);