2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013,2014.2015,2017,2018, 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/gmxlib/nrnb.h"
48 * Gromacs nonbonded kernel: nb_kernel_ElecEw_VdwCSTab_GeomW4W4_VF_c
49 * Electrostatics interaction: Ewald
50 * VdW interaction: CubicSplineTable
51 * Geometry: Water4-Water4
52 * Calculate force/pot: PotentialAndForce
55 nb_kernel_ElecEw_VdwCSTab_GeomW4W4_VF_c
56 (t_nblist
* gmx_restrict nlist
,
57 rvec
* gmx_restrict xx
,
58 rvec
* gmx_restrict ff
,
59 struct t_forcerec
* gmx_restrict fr
,
60 t_mdatoms
* gmx_restrict mdatoms
,
61 nb_kernel_data_t gmx_unused
* gmx_restrict kernel_data
,
62 t_nrnb
* gmx_restrict nrnb
)
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
;
106 real ewtabscale
,eweps
,sh_ewald
,ewrt
,ewtabhalfspace
;
114 jindex
= nlist
->jindex
;
116 shiftidx
= nlist
->shift
;
118 shiftvec
= fr
->shift_vec
[0];
119 fshift
= fr
->fshift
[0];
120 facel
= fr
->ic
->epsfac
;
121 charge
= mdatoms
->chargeA
;
122 nvdwtype
= fr
->ntype
;
124 vdwtype
= mdatoms
->typeA
;
126 vftab
= kernel_data
->table_vdw
->data
;
127 vftabscale
= kernel_data
->table_vdw
->scale
;
129 sh_ewald
= fr
->ic
->sh_ewald
;
130 ewtab
= fr
->ic
->tabq_coul_FDV0
;
131 ewtabscale
= fr
->ic
->tabq_scale
;
132 ewtabhalfspace
= 0.5/ewtabscale
;
134 /* Setup water-specific parameters */
135 inr
= nlist
->iinr
[0];
136 iq1
= facel
*charge
[inr
+1];
137 iq2
= facel
*charge
[inr
+2];
138 iq3
= facel
*charge
[inr
+3];
139 vdwioffset0
= 2*nvdwtype
*vdwtype
[inr
+0];
144 vdwjidx0
= 2*vdwtype
[inr
+0];
145 c6_00
= vdwparam
[vdwioffset0
+vdwjidx0
];
146 c12_00
= vdwparam
[vdwioffset0
+vdwjidx0
+1];
160 /* Start outer loop over neighborlists */
161 for(iidx
=0; iidx
<nri
; iidx
++)
163 /* Load shift vector for this list */
164 i_shift_offset
= DIM
*shiftidx
[iidx
];
165 shX
= shiftvec
[i_shift_offset
+XX
];
166 shY
= shiftvec
[i_shift_offset
+YY
];
167 shZ
= shiftvec
[i_shift_offset
+ZZ
];
169 /* Load limits for loop over neighbors */
170 j_index_start
= jindex
[iidx
];
171 j_index_end
= jindex
[iidx
+1];
173 /* Get outer coordinate index */
175 i_coord_offset
= DIM
*inr
;
177 /* Load i particle coords and add shift vector */
178 ix0
= shX
+ x
[i_coord_offset
+DIM
*0+XX
];
179 iy0
= shY
+ x
[i_coord_offset
+DIM
*0+YY
];
180 iz0
= shZ
+ x
[i_coord_offset
+DIM
*0+ZZ
];
181 ix1
= shX
+ x
[i_coord_offset
+DIM
*1+XX
];
182 iy1
= shY
+ x
[i_coord_offset
+DIM
*1+YY
];
183 iz1
= shZ
+ x
[i_coord_offset
+DIM
*1+ZZ
];
184 ix2
= shX
+ x
[i_coord_offset
+DIM
*2+XX
];
185 iy2
= shY
+ x
[i_coord_offset
+DIM
*2+YY
];
186 iz2
= shZ
+ x
[i_coord_offset
+DIM
*2+ZZ
];
187 ix3
= shX
+ x
[i_coord_offset
+DIM
*3+XX
];
188 iy3
= shY
+ x
[i_coord_offset
+DIM
*3+YY
];
189 iz3
= shZ
+ x
[i_coord_offset
+DIM
*3+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
];
225 jx3
= x
[j_coord_offset
+DIM
*3+XX
];
226 jy3
= x
[j_coord_offset
+DIM
*3+YY
];
227 jz3
= x
[j_coord_offset
+DIM
*3+ZZ
];
229 /* Calculate displacement vector */
261 /* Calculate squared distance and things based on it */
262 rsq00
= dx00
*dx00
+dy00
*dy00
+dz00
*dz00
;
263 rsq11
= dx11
*dx11
+dy11
*dy11
+dz11
*dz11
;
264 rsq12
= dx12
*dx12
+dy12
*dy12
+dz12
*dz12
;
265 rsq13
= dx13
*dx13
+dy13
*dy13
+dz13
*dz13
;
266 rsq21
= dx21
*dx21
+dy21
*dy21
+dz21
*dz21
;
267 rsq22
= dx22
*dx22
+dy22
*dy22
+dz22
*dz22
;
268 rsq23
= dx23
*dx23
+dy23
*dy23
+dz23
*dz23
;
269 rsq31
= dx31
*dx31
+dy31
*dy31
+dz31
*dz31
;
270 rsq32
= dx32
*dx32
+dy32
*dy32
+dz32
*dz32
;
271 rsq33
= dx33
*dx33
+dy33
*dy33
+dz33
*dz33
;
273 rinv00
= 1.0/sqrt(rsq00
);
274 rinv11
= 1.0/sqrt(rsq11
);
275 rinv12
= 1.0/sqrt(rsq12
);
276 rinv13
= 1.0/sqrt(rsq13
);
277 rinv21
= 1.0/sqrt(rsq21
);
278 rinv22
= 1.0/sqrt(rsq22
);
279 rinv23
= 1.0/sqrt(rsq23
);
280 rinv31
= 1.0/sqrt(rsq31
);
281 rinv32
= 1.0/sqrt(rsq32
);
282 rinv33
= 1.0/sqrt(rsq33
);
284 rinvsq11
= rinv11
*rinv11
;
285 rinvsq12
= rinv12
*rinv12
;
286 rinvsq13
= rinv13
*rinv13
;
287 rinvsq21
= rinv21
*rinv21
;
288 rinvsq22
= rinv22
*rinv22
;
289 rinvsq23
= rinv23
*rinv23
;
290 rinvsq31
= rinv31
*rinv31
;
291 rinvsq32
= rinv32
*rinv32
;
292 rinvsq33
= rinv33
*rinv33
;
294 /**************************
295 * CALCULATE INTERACTIONS *
296 **************************/
300 /* Calculate table index by multiplying r with table scale and truncate to integer */
306 /* CUBIC SPLINE TABLE DISPERSION */
310 Geps
= vfeps
*vftab
[vfitab
+2];
311 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
315 FF
= Fp
+Geps
+2.0*Heps2
;
318 /* CUBIC SPLINE TABLE REPULSION */
321 Geps
= vfeps
*vftab
[vfitab
+6];
322 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+7];
326 FF
= Fp
+Geps
+2.0*Heps2
;
329 fvdw
= -(fvdw6
+fvdw12
)*vftabscale
*rinv00
;
331 /* Update potential sums from outer loop */
336 /* Calculate temporary vectorial force */
341 /* Update vectorial force */
345 f
[j_coord_offset
+DIM
*0+XX
] -= tx
;
346 f
[j_coord_offset
+DIM
*0+YY
] -= ty
;
347 f
[j_coord_offset
+DIM
*0+ZZ
] -= tz
;
349 /**************************
350 * CALCULATE INTERACTIONS *
351 **************************/
355 /* EWALD ELECTROSTATICS */
357 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
358 ewrt
= r11
*ewtabscale
;
362 felec
= ewtab
[ewitab
]+eweps
*ewtab
[ewitab
+1];
363 velec
= qq11
*(rinv11
-(ewtab
[ewitab
+2]-ewtabhalfspace
*eweps
*(ewtab
[ewitab
]+felec
)));
364 felec
= qq11
*rinv11
*(rinvsq11
-felec
);
366 /* Update potential sums from outer loop */
371 /* Calculate temporary vectorial force */
376 /* Update vectorial force */
380 f
[j_coord_offset
+DIM
*1+XX
] -= tx
;
381 f
[j_coord_offset
+DIM
*1+YY
] -= ty
;
382 f
[j_coord_offset
+DIM
*1+ZZ
] -= tz
;
384 /**************************
385 * CALCULATE INTERACTIONS *
386 **************************/
390 /* EWALD ELECTROSTATICS */
392 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
393 ewrt
= r12
*ewtabscale
;
397 felec
= ewtab
[ewitab
]+eweps
*ewtab
[ewitab
+1];
398 velec
= qq12
*(rinv12
-(ewtab
[ewitab
+2]-ewtabhalfspace
*eweps
*(ewtab
[ewitab
]+felec
)));
399 felec
= qq12
*rinv12
*(rinvsq12
-felec
);
401 /* Update potential sums from outer loop */
406 /* Calculate temporary vectorial force */
411 /* Update vectorial force */
415 f
[j_coord_offset
+DIM
*2+XX
] -= tx
;
416 f
[j_coord_offset
+DIM
*2+YY
] -= ty
;
417 f
[j_coord_offset
+DIM
*2+ZZ
] -= tz
;
419 /**************************
420 * CALCULATE INTERACTIONS *
421 **************************/
425 /* EWALD ELECTROSTATICS */
427 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
428 ewrt
= r13
*ewtabscale
;
432 felec
= ewtab
[ewitab
]+eweps
*ewtab
[ewitab
+1];
433 velec
= qq13
*(rinv13
-(ewtab
[ewitab
+2]-ewtabhalfspace
*eweps
*(ewtab
[ewitab
]+felec
)));
434 felec
= qq13
*rinv13
*(rinvsq13
-felec
);
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 /* EWALD ELECTROSTATICS */
462 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
463 ewrt
= r21
*ewtabscale
;
467 felec
= ewtab
[ewitab
]+eweps
*ewtab
[ewitab
+1];
468 velec
= qq21
*(rinv21
-(ewtab
[ewitab
+2]-ewtabhalfspace
*eweps
*(ewtab
[ewitab
]+felec
)));
469 felec
= qq21
*rinv21
*(rinvsq21
-felec
);
471 /* Update potential sums from outer loop */
476 /* Calculate temporary vectorial force */
481 /* Update vectorial force */
485 f
[j_coord_offset
+DIM
*1+XX
] -= tx
;
486 f
[j_coord_offset
+DIM
*1+YY
] -= ty
;
487 f
[j_coord_offset
+DIM
*1+ZZ
] -= tz
;
489 /**************************
490 * CALCULATE INTERACTIONS *
491 **************************/
495 /* EWALD ELECTROSTATICS */
497 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
498 ewrt
= r22
*ewtabscale
;
502 felec
= ewtab
[ewitab
]+eweps
*ewtab
[ewitab
+1];
503 velec
= qq22
*(rinv22
-(ewtab
[ewitab
+2]-ewtabhalfspace
*eweps
*(ewtab
[ewitab
]+felec
)));
504 felec
= qq22
*rinv22
*(rinvsq22
-felec
);
506 /* Update potential sums from outer loop */
511 /* Calculate temporary vectorial force */
516 /* Update vectorial force */
520 f
[j_coord_offset
+DIM
*2+XX
] -= tx
;
521 f
[j_coord_offset
+DIM
*2+YY
] -= ty
;
522 f
[j_coord_offset
+DIM
*2+ZZ
] -= tz
;
524 /**************************
525 * CALCULATE INTERACTIONS *
526 **************************/
530 /* EWALD ELECTROSTATICS */
532 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
533 ewrt
= r23
*ewtabscale
;
537 felec
= ewtab
[ewitab
]+eweps
*ewtab
[ewitab
+1];
538 velec
= qq23
*(rinv23
-(ewtab
[ewitab
+2]-ewtabhalfspace
*eweps
*(ewtab
[ewitab
]+felec
)));
539 felec
= qq23
*rinv23
*(rinvsq23
-felec
);
541 /* Update potential sums from outer loop */
546 /* Calculate temporary vectorial force */
551 /* Update vectorial force */
555 f
[j_coord_offset
+DIM
*3+XX
] -= tx
;
556 f
[j_coord_offset
+DIM
*3+YY
] -= ty
;
557 f
[j_coord_offset
+DIM
*3+ZZ
] -= tz
;
559 /**************************
560 * CALCULATE INTERACTIONS *
561 **************************/
565 /* EWALD ELECTROSTATICS */
567 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
568 ewrt
= r31
*ewtabscale
;
572 felec
= ewtab
[ewitab
]+eweps
*ewtab
[ewitab
+1];
573 velec
= qq31
*(rinv31
-(ewtab
[ewitab
+2]-ewtabhalfspace
*eweps
*(ewtab
[ewitab
]+felec
)));
574 felec
= qq31
*rinv31
*(rinvsq31
-felec
);
576 /* Update potential sums from outer loop */
581 /* Calculate temporary vectorial force */
586 /* Update vectorial force */
590 f
[j_coord_offset
+DIM
*1+XX
] -= tx
;
591 f
[j_coord_offset
+DIM
*1+YY
] -= ty
;
592 f
[j_coord_offset
+DIM
*1+ZZ
] -= tz
;
594 /**************************
595 * CALCULATE INTERACTIONS *
596 **************************/
600 /* EWALD ELECTROSTATICS */
602 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
603 ewrt
= r32
*ewtabscale
;
607 felec
= ewtab
[ewitab
]+eweps
*ewtab
[ewitab
+1];
608 velec
= qq32
*(rinv32
-(ewtab
[ewitab
+2]-ewtabhalfspace
*eweps
*(ewtab
[ewitab
]+felec
)));
609 felec
= qq32
*rinv32
*(rinvsq32
-felec
);
611 /* Update potential sums from outer loop */
616 /* Calculate temporary vectorial force */
621 /* Update vectorial force */
625 f
[j_coord_offset
+DIM
*2+XX
] -= tx
;
626 f
[j_coord_offset
+DIM
*2+YY
] -= ty
;
627 f
[j_coord_offset
+DIM
*2+ZZ
] -= tz
;
629 /**************************
630 * CALCULATE INTERACTIONS *
631 **************************/
635 /* EWALD ELECTROSTATICS */
637 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
638 ewrt
= r33
*ewtabscale
;
642 felec
= ewtab
[ewitab
]+eweps
*ewtab
[ewitab
+1];
643 velec
= qq33
*(rinv33
-(ewtab
[ewitab
+2]-ewtabhalfspace
*eweps
*(ewtab
[ewitab
]+felec
)));
644 felec
= qq33
*rinv33
*(rinvsq33
-felec
);
646 /* Update potential sums from outer loop */
651 /* Calculate temporary vectorial force */
656 /* Update vectorial force */
660 f
[j_coord_offset
+DIM
*3+XX
] -= tx
;
661 f
[j_coord_offset
+DIM
*3+YY
] -= ty
;
662 f
[j_coord_offset
+DIM
*3+ZZ
] -= tz
;
664 /* Inner loop uses 415 flops */
666 /* End of innermost loop */
669 f
[i_coord_offset
+DIM
*0+XX
] += fix0
;
670 f
[i_coord_offset
+DIM
*0+YY
] += fiy0
;
671 f
[i_coord_offset
+DIM
*0+ZZ
] += fiz0
;
675 f
[i_coord_offset
+DIM
*1+XX
] += fix1
;
676 f
[i_coord_offset
+DIM
*1+YY
] += fiy1
;
677 f
[i_coord_offset
+DIM
*1+ZZ
] += fiz1
;
681 f
[i_coord_offset
+DIM
*2+XX
] += fix2
;
682 f
[i_coord_offset
+DIM
*2+YY
] += fiy2
;
683 f
[i_coord_offset
+DIM
*2+ZZ
] += fiz2
;
687 f
[i_coord_offset
+DIM
*3+XX
] += fix3
;
688 f
[i_coord_offset
+DIM
*3+YY
] += fiy3
;
689 f
[i_coord_offset
+DIM
*3+ZZ
] += fiz3
;
693 fshift
[i_shift_offset
+XX
] += tx
;
694 fshift
[i_shift_offset
+YY
] += ty
;
695 fshift
[i_shift_offset
+ZZ
] += tz
;
698 /* Update potential energies */
699 kernel_data
->energygrp_elec
[ggid
] += velecsum
;
700 kernel_data
->energygrp_vdw
[ggid
] += vvdwsum
;
702 /* Increment number of inner iterations */
703 inneriter
+= j_index_end
- j_index_start
;
705 /* Outer loop uses 41 flops */
708 /* Increment number of outer iterations */
711 /* Update outer/inner flops */
713 inc_nrnb(nrnb
,eNR_NBKERNEL_ELEC_VDW_W4W4_VF
,outeriter
*41 + inneriter
*415);
716 * Gromacs nonbonded kernel: nb_kernel_ElecEw_VdwCSTab_GeomW4W4_F_c
717 * Electrostatics interaction: Ewald
718 * VdW interaction: CubicSplineTable
719 * Geometry: Water4-Water4
720 * Calculate force/pot: Force
723 nb_kernel_ElecEw_VdwCSTab_GeomW4W4_F_c
724 (t_nblist
* gmx_restrict nlist
,
725 rvec
* gmx_restrict xx
,
726 rvec
* gmx_restrict ff
,
727 struct t_forcerec
* gmx_restrict fr
,
728 t_mdatoms
* gmx_restrict mdatoms
,
729 nb_kernel_data_t gmx_unused
* gmx_restrict kernel_data
,
730 t_nrnb
* gmx_restrict nrnb
)
732 int i_shift_offset
,i_coord_offset
,j_coord_offset
;
733 int j_index_start
,j_index_end
;
734 int nri
,inr
,ggid
,iidx
,jidx
,jnr
,outeriter
,inneriter
;
735 real shX
,shY
,shZ
,tx
,ty
,tz
,fscal
,rcutoff
,rcutoff2
;
736 int *iinr
,*jindex
,*jjnr
,*shiftidx
,*gid
;
737 real
*shiftvec
,*fshift
,*x
,*f
;
739 real ix0
,iy0
,iz0
,fix0
,fiy0
,fiz0
,iq0
,isai0
;
741 real ix1
,iy1
,iz1
,fix1
,fiy1
,fiz1
,iq1
,isai1
;
743 real ix2
,iy2
,iz2
,fix2
,fiy2
,fiz2
,iq2
,isai2
;
745 real ix3
,iy3
,iz3
,fix3
,fiy3
,fiz3
,iq3
,isai3
;
747 real jx0
,jy0
,jz0
,fjx0
,fjy0
,fjz0
,jq0
,isaj0
;
749 real jx1
,jy1
,jz1
,fjx1
,fjy1
,fjz1
,jq1
,isaj1
;
751 real jx2
,jy2
,jz2
,fjx2
,fjy2
,fjz2
,jq2
,isaj2
;
753 real jx3
,jy3
,jz3
,fjx3
,fjy3
,fjz3
,jq3
,isaj3
;
754 real dx00
,dy00
,dz00
,rsq00
,rinv00
,rinvsq00
,r00
,qq00
,c6_00
,c12_00
,cexp1_00
,cexp2_00
;
755 real dx11
,dy11
,dz11
,rsq11
,rinv11
,rinvsq11
,r11
,qq11
,c6_11
,c12_11
,cexp1_11
,cexp2_11
;
756 real dx12
,dy12
,dz12
,rsq12
,rinv12
,rinvsq12
,r12
,qq12
,c6_12
,c12_12
,cexp1_12
,cexp2_12
;
757 real dx13
,dy13
,dz13
,rsq13
,rinv13
,rinvsq13
,r13
,qq13
,c6_13
,c12_13
,cexp1_13
,cexp2_13
;
758 real dx21
,dy21
,dz21
,rsq21
,rinv21
,rinvsq21
,r21
,qq21
,c6_21
,c12_21
,cexp1_21
,cexp2_21
;
759 real dx22
,dy22
,dz22
,rsq22
,rinv22
,rinvsq22
,r22
,qq22
,c6_22
,c12_22
,cexp1_22
,cexp2_22
;
760 real dx23
,dy23
,dz23
,rsq23
,rinv23
,rinvsq23
,r23
,qq23
,c6_23
,c12_23
,cexp1_23
,cexp2_23
;
761 real dx31
,dy31
,dz31
,rsq31
,rinv31
,rinvsq31
,r31
,qq31
,c6_31
,c12_31
,cexp1_31
,cexp2_31
;
762 real dx32
,dy32
,dz32
,rsq32
,rinv32
,rinvsq32
,r32
,qq32
,c6_32
,c12_32
,cexp1_32
,cexp2_32
;
763 real dx33
,dy33
,dz33
,rsq33
,rinv33
,rinvsq33
,r33
,qq33
,c6_33
,c12_33
,cexp1_33
,cexp2_33
;
764 real velec
,felec
,velecsum
,facel
,crf
,krf
,krf2
;
767 real rinvsix
,rvdw
,vvdw
,vvdw6
,vvdw12
,fvdw
,fvdw6
,fvdw12
,vvdwsum
,br
,vvdwexp
,sh_vdw_invrcut6
;
771 real rt
,vfeps
,vftabscale
,Y
,F
,Geps
,Heps2
,Fp
,VV
,FF
;
774 real ewtabscale
,eweps
,sh_ewald
,ewrt
,ewtabhalfspace
;
782 jindex
= nlist
->jindex
;
784 shiftidx
= nlist
->shift
;
786 shiftvec
= fr
->shift_vec
[0];
787 fshift
= fr
->fshift
[0];
788 facel
= fr
->ic
->epsfac
;
789 charge
= mdatoms
->chargeA
;
790 nvdwtype
= fr
->ntype
;
792 vdwtype
= mdatoms
->typeA
;
794 vftab
= kernel_data
->table_vdw
->data
;
795 vftabscale
= kernel_data
->table_vdw
->scale
;
797 sh_ewald
= fr
->ic
->sh_ewald
;
798 ewtab
= fr
->ic
->tabq_coul_F
;
799 ewtabscale
= fr
->ic
->tabq_scale
;
800 ewtabhalfspace
= 0.5/ewtabscale
;
802 /* Setup water-specific parameters */
803 inr
= nlist
->iinr
[0];
804 iq1
= facel
*charge
[inr
+1];
805 iq2
= facel
*charge
[inr
+2];
806 iq3
= facel
*charge
[inr
+3];
807 vdwioffset0
= 2*nvdwtype
*vdwtype
[inr
+0];
812 vdwjidx0
= 2*vdwtype
[inr
+0];
813 c6_00
= vdwparam
[vdwioffset0
+vdwjidx0
];
814 c12_00
= vdwparam
[vdwioffset0
+vdwjidx0
+1];
828 /* Start outer loop over neighborlists */
829 for(iidx
=0; iidx
<nri
; iidx
++)
831 /* Load shift vector for this list */
832 i_shift_offset
= DIM
*shiftidx
[iidx
];
833 shX
= shiftvec
[i_shift_offset
+XX
];
834 shY
= shiftvec
[i_shift_offset
+YY
];
835 shZ
= shiftvec
[i_shift_offset
+ZZ
];
837 /* Load limits for loop over neighbors */
838 j_index_start
= jindex
[iidx
];
839 j_index_end
= jindex
[iidx
+1];
841 /* Get outer coordinate index */
843 i_coord_offset
= DIM
*inr
;
845 /* Load i particle coords and add shift vector */
846 ix0
= shX
+ x
[i_coord_offset
+DIM
*0+XX
];
847 iy0
= shY
+ x
[i_coord_offset
+DIM
*0+YY
];
848 iz0
= shZ
+ x
[i_coord_offset
+DIM
*0+ZZ
];
849 ix1
= shX
+ x
[i_coord_offset
+DIM
*1+XX
];
850 iy1
= shY
+ x
[i_coord_offset
+DIM
*1+YY
];
851 iz1
= shZ
+ x
[i_coord_offset
+DIM
*1+ZZ
];
852 ix2
= shX
+ x
[i_coord_offset
+DIM
*2+XX
];
853 iy2
= shY
+ x
[i_coord_offset
+DIM
*2+YY
];
854 iz2
= shZ
+ x
[i_coord_offset
+DIM
*2+ZZ
];
855 ix3
= shX
+ x
[i_coord_offset
+DIM
*3+XX
];
856 iy3
= shY
+ x
[i_coord_offset
+DIM
*3+YY
];
857 iz3
= shZ
+ x
[i_coord_offset
+DIM
*3+ZZ
];
872 /* Start inner kernel loop */
873 for(jidx
=j_index_start
; jidx
<j_index_end
; jidx
++)
875 /* Get j neighbor index, and coordinate index */
877 j_coord_offset
= DIM
*jnr
;
879 /* load j atom coordinates */
880 jx0
= x
[j_coord_offset
+DIM
*0+XX
];
881 jy0
= x
[j_coord_offset
+DIM
*0+YY
];
882 jz0
= x
[j_coord_offset
+DIM
*0+ZZ
];
883 jx1
= x
[j_coord_offset
+DIM
*1+XX
];
884 jy1
= x
[j_coord_offset
+DIM
*1+YY
];
885 jz1
= x
[j_coord_offset
+DIM
*1+ZZ
];
886 jx2
= x
[j_coord_offset
+DIM
*2+XX
];
887 jy2
= x
[j_coord_offset
+DIM
*2+YY
];
888 jz2
= x
[j_coord_offset
+DIM
*2+ZZ
];
889 jx3
= x
[j_coord_offset
+DIM
*3+XX
];
890 jy3
= x
[j_coord_offset
+DIM
*3+YY
];
891 jz3
= x
[j_coord_offset
+DIM
*3+ZZ
];
893 /* Calculate displacement vector */
925 /* Calculate squared distance and things based on it */
926 rsq00
= dx00
*dx00
+dy00
*dy00
+dz00
*dz00
;
927 rsq11
= dx11
*dx11
+dy11
*dy11
+dz11
*dz11
;
928 rsq12
= dx12
*dx12
+dy12
*dy12
+dz12
*dz12
;
929 rsq13
= dx13
*dx13
+dy13
*dy13
+dz13
*dz13
;
930 rsq21
= dx21
*dx21
+dy21
*dy21
+dz21
*dz21
;
931 rsq22
= dx22
*dx22
+dy22
*dy22
+dz22
*dz22
;
932 rsq23
= dx23
*dx23
+dy23
*dy23
+dz23
*dz23
;
933 rsq31
= dx31
*dx31
+dy31
*dy31
+dz31
*dz31
;
934 rsq32
= dx32
*dx32
+dy32
*dy32
+dz32
*dz32
;
935 rsq33
= dx33
*dx33
+dy33
*dy33
+dz33
*dz33
;
937 rinv00
= 1.0/sqrt(rsq00
);
938 rinv11
= 1.0/sqrt(rsq11
);
939 rinv12
= 1.0/sqrt(rsq12
);
940 rinv13
= 1.0/sqrt(rsq13
);
941 rinv21
= 1.0/sqrt(rsq21
);
942 rinv22
= 1.0/sqrt(rsq22
);
943 rinv23
= 1.0/sqrt(rsq23
);
944 rinv31
= 1.0/sqrt(rsq31
);
945 rinv32
= 1.0/sqrt(rsq32
);
946 rinv33
= 1.0/sqrt(rsq33
);
948 rinvsq11
= rinv11
*rinv11
;
949 rinvsq12
= rinv12
*rinv12
;
950 rinvsq13
= rinv13
*rinv13
;
951 rinvsq21
= rinv21
*rinv21
;
952 rinvsq22
= rinv22
*rinv22
;
953 rinvsq23
= rinv23
*rinv23
;
954 rinvsq31
= rinv31
*rinv31
;
955 rinvsq32
= rinv32
*rinv32
;
956 rinvsq33
= rinv33
*rinv33
;
958 /**************************
959 * CALCULATE INTERACTIONS *
960 **************************/
964 /* Calculate table index by multiplying r with table scale and truncate to integer */
970 /* CUBIC SPLINE TABLE DISPERSION */
973 Geps
= vfeps
*vftab
[vfitab
+2];
974 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+3];
976 FF
= Fp
+Geps
+2.0*Heps2
;
979 /* CUBIC SPLINE TABLE REPULSION */
981 Geps
= vfeps
*vftab
[vfitab
+6];
982 Heps2
= vfeps
*vfeps
*vftab
[vfitab
+7];
984 FF
= Fp
+Geps
+2.0*Heps2
;
986 fvdw
= -(fvdw6
+fvdw12
)*vftabscale
*rinv00
;
990 /* Calculate temporary vectorial force */
995 /* Update vectorial force */
999 f
[j_coord_offset
+DIM
*0+XX
] -= tx
;
1000 f
[j_coord_offset
+DIM
*0+YY
] -= ty
;
1001 f
[j_coord_offset
+DIM
*0+ZZ
] -= tz
;
1003 /**************************
1004 * CALCULATE INTERACTIONS *
1005 **************************/
1009 /* EWALD ELECTROSTATICS */
1011 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1012 ewrt
= r11
*ewtabscale
;
1014 eweps
= ewrt
-ewitab
;
1015 felec
= (1.0-eweps
)*ewtab
[ewitab
]+eweps
*ewtab
[ewitab
+1];
1016 felec
= qq11
*rinv11
*(rinvsq11
-felec
);
1020 /* Calculate temporary vectorial force */
1025 /* Update vectorial force */
1029 f
[j_coord_offset
+DIM
*1+XX
] -= tx
;
1030 f
[j_coord_offset
+DIM
*1+YY
] -= ty
;
1031 f
[j_coord_offset
+DIM
*1+ZZ
] -= tz
;
1033 /**************************
1034 * CALCULATE INTERACTIONS *
1035 **************************/
1039 /* EWALD ELECTROSTATICS */
1041 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1042 ewrt
= r12
*ewtabscale
;
1044 eweps
= ewrt
-ewitab
;
1045 felec
= (1.0-eweps
)*ewtab
[ewitab
]+eweps
*ewtab
[ewitab
+1];
1046 felec
= qq12
*rinv12
*(rinvsq12
-felec
);
1050 /* Calculate temporary vectorial force */
1055 /* Update vectorial force */
1059 f
[j_coord_offset
+DIM
*2+XX
] -= tx
;
1060 f
[j_coord_offset
+DIM
*2+YY
] -= ty
;
1061 f
[j_coord_offset
+DIM
*2+ZZ
] -= tz
;
1063 /**************************
1064 * CALCULATE INTERACTIONS *
1065 **************************/
1069 /* EWALD ELECTROSTATICS */
1071 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1072 ewrt
= r13
*ewtabscale
;
1074 eweps
= ewrt
-ewitab
;
1075 felec
= (1.0-eweps
)*ewtab
[ewitab
]+eweps
*ewtab
[ewitab
+1];
1076 felec
= qq13
*rinv13
*(rinvsq13
-felec
);
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
;
1093 /**************************
1094 * CALCULATE INTERACTIONS *
1095 **************************/
1099 /* EWALD ELECTROSTATICS */
1101 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1102 ewrt
= r21
*ewtabscale
;
1104 eweps
= ewrt
-ewitab
;
1105 felec
= (1.0-eweps
)*ewtab
[ewitab
]+eweps
*ewtab
[ewitab
+1];
1106 felec
= qq21
*rinv21
*(rinvsq21
-felec
);
1110 /* Calculate temporary vectorial force */
1115 /* Update vectorial force */
1119 f
[j_coord_offset
+DIM
*1+XX
] -= tx
;
1120 f
[j_coord_offset
+DIM
*1+YY
] -= ty
;
1121 f
[j_coord_offset
+DIM
*1+ZZ
] -= tz
;
1123 /**************************
1124 * CALCULATE INTERACTIONS *
1125 **************************/
1129 /* EWALD ELECTROSTATICS */
1131 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1132 ewrt
= r22
*ewtabscale
;
1134 eweps
= ewrt
-ewitab
;
1135 felec
= (1.0-eweps
)*ewtab
[ewitab
]+eweps
*ewtab
[ewitab
+1];
1136 felec
= qq22
*rinv22
*(rinvsq22
-felec
);
1140 /* Calculate temporary vectorial force */
1145 /* Update vectorial force */
1149 f
[j_coord_offset
+DIM
*2+XX
] -= tx
;
1150 f
[j_coord_offset
+DIM
*2+YY
] -= ty
;
1151 f
[j_coord_offset
+DIM
*2+ZZ
] -= tz
;
1153 /**************************
1154 * CALCULATE INTERACTIONS *
1155 **************************/
1159 /* EWALD ELECTROSTATICS */
1161 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1162 ewrt
= r23
*ewtabscale
;
1164 eweps
= ewrt
-ewitab
;
1165 felec
= (1.0-eweps
)*ewtab
[ewitab
]+eweps
*ewtab
[ewitab
+1];
1166 felec
= qq23
*rinv23
*(rinvsq23
-felec
);
1170 /* Calculate temporary vectorial force */
1175 /* Update vectorial force */
1179 f
[j_coord_offset
+DIM
*3+XX
] -= tx
;
1180 f
[j_coord_offset
+DIM
*3+YY
] -= ty
;
1181 f
[j_coord_offset
+DIM
*3+ZZ
] -= tz
;
1183 /**************************
1184 * CALCULATE INTERACTIONS *
1185 **************************/
1189 /* EWALD ELECTROSTATICS */
1191 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1192 ewrt
= r31
*ewtabscale
;
1194 eweps
= ewrt
-ewitab
;
1195 felec
= (1.0-eweps
)*ewtab
[ewitab
]+eweps
*ewtab
[ewitab
+1];
1196 felec
= qq31
*rinv31
*(rinvsq31
-felec
);
1200 /* Calculate temporary vectorial force */
1205 /* Update vectorial force */
1209 f
[j_coord_offset
+DIM
*1+XX
] -= tx
;
1210 f
[j_coord_offset
+DIM
*1+YY
] -= ty
;
1211 f
[j_coord_offset
+DIM
*1+ZZ
] -= tz
;
1213 /**************************
1214 * CALCULATE INTERACTIONS *
1215 **************************/
1219 /* EWALD ELECTROSTATICS */
1221 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1222 ewrt
= r32
*ewtabscale
;
1224 eweps
= ewrt
-ewitab
;
1225 felec
= (1.0-eweps
)*ewtab
[ewitab
]+eweps
*ewtab
[ewitab
+1];
1226 felec
= qq32
*rinv32
*(rinvsq32
-felec
);
1230 /* Calculate temporary vectorial force */
1235 /* Update vectorial force */
1239 f
[j_coord_offset
+DIM
*2+XX
] -= tx
;
1240 f
[j_coord_offset
+DIM
*2+YY
] -= ty
;
1241 f
[j_coord_offset
+DIM
*2+ZZ
] -= tz
;
1243 /**************************
1244 * CALCULATE INTERACTIONS *
1245 **************************/
1249 /* EWALD ELECTROSTATICS */
1251 /* Calculate Ewald table index by multiplying r with scale and truncate to integer */
1252 ewrt
= r33
*ewtabscale
;
1254 eweps
= ewrt
-ewitab
;
1255 felec
= (1.0-eweps
)*ewtab
[ewitab
]+eweps
*ewtab
[ewitab
+1];
1256 felec
= qq33
*rinv33
*(rinvsq33
-felec
);
1260 /* Calculate temporary vectorial force */
1265 /* Update vectorial force */
1269 f
[j_coord_offset
+DIM
*3+XX
] -= tx
;
1270 f
[j_coord_offset
+DIM
*3+YY
] -= ty
;
1271 f
[j_coord_offset
+DIM
*3+ZZ
] -= tz
;
1273 /* Inner loop uses 344 flops */
1275 /* End of innermost loop */
1278 f
[i_coord_offset
+DIM
*0+XX
] += fix0
;
1279 f
[i_coord_offset
+DIM
*0+YY
] += fiy0
;
1280 f
[i_coord_offset
+DIM
*0+ZZ
] += fiz0
;
1284 f
[i_coord_offset
+DIM
*1+XX
] += fix1
;
1285 f
[i_coord_offset
+DIM
*1+YY
] += fiy1
;
1286 f
[i_coord_offset
+DIM
*1+ZZ
] += fiz1
;
1290 f
[i_coord_offset
+DIM
*2+XX
] += fix2
;
1291 f
[i_coord_offset
+DIM
*2+YY
] += fiy2
;
1292 f
[i_coord_offset
+DIM
*2+ZZ
] += fiz2
;
1296 f
[i_coord_offset
+DIM
*3+XX
] += fix3
;
1297 f
[i_coord_offset
+DIM
*3+YY
] += fiy3
;
1298 f
[i_coord_offset
+DIM
*3+ZZ
] += fiz3
;
1302 fshift
[i_shift_offset
+XX
] += tx
;
1303 fshift
[i_shift_offset
+YY
] += ty
;
1304 fshift
[i_shift_offset
+ZZ
] += tz
;
1306 /* Increment number of inner iterations */
1307 inneriter
+= j_index_end
- j_index_start
;
1309 /* Outer loop uses 39 flops */
1312 /* Increment number of outer iterations */
1315 /* Update outer/inner flops */
1317 inc_nrnb(nrnb
,eNR_NBKERNEL_ELEC_VDW_W4W4_F
,outeriter
*39 + inneriter
*344);