Avoid using function calls in OpenMP directives
[alexxy/gromacs.git] / src / gromacs / mdlib / nbnxn_kernels / simd_4xn / nbnxn_kernel_simd_4xn.c
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
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.
8  *
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.
13  *
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.
18  *
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.
23  *
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.
31  *
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.
34  */
35 /*
36  * Note: this file was generated by the Verlet kernel generator for
37  * kernel type 4xn.
38  */
39
40 #ifdef HAVE_CONFIG_H
41 #include <config.h>
42 #endif
43
44 #include "typedefs.h"
45
46 #include "gromacs/mdlib/nbnxn_simd.h"
47
48 #ifdef GMX_NBNXN_SIMD_4XN
49
50 #include "gromacs/simd/vector_operations.h"
51
52 #if !(GMX_SIMD_REAL_WIDTH == 2 || GMX_SIMD_REAL_WIDTH == 4 || GMX_SIMD_REAL_WIDTH == 8)
53 #error "unsupported SIMD width"
54 #endif
55
56 #define GMX_SIMD_J_UNROLL_SIZE 1
57 #include "nbnxn_kernel_simd_4xn.h"
58 #include "../nbnxn_kernel_common.h"
59 #include "gmx_omp_nthreads.h"
60 #include "types/force_flags.h"
61 #include "gmx_fatal.h"
62
63 /*! \brief Kinds of electrostatic treatments in SIMD Verlet kernels
64  */
65 enum {
66     coulktRF, coulktTAB, coulktTAB_TWIN, coulktEWALD, coulktEWALD_TWIN, coulktNR
67 };
68
69 /*! \brief Kinds of Van der Waals treatments in SIMD Verlet kernels
70  */
71 enum {
72     vdwktLJCUT_COMBGEOM, vdwktLJCUT_COMBLB, vdwktLJCUT_COMBNONE, vdwktLJFORCESWITCH, vdwktLJPOTSWITCH, vdwktLJEWALDCOMBGEOM, vdwktNR
73 };
74
75 /* Declare and define the kernel function pointer lookup tables.
76  * The minor index of the array goes over both the LJ combination rules,
77  * which is only supported by plain cut-off, and the LJ switch/PME functions.
78  */
79 static p_nbk_func_noener p_nbk_noener[coulktNR][vdwktNR] =
80 {
81     {
82         nbnxn_kernel_ElecRF_VdwLJCombGeom_F_4xn,
83         nbnxn_kernel_ElecRF_VdwLJCombLB_F_4xn,
84         nbnxn_kernel_ElecRF_VdwLJ_F_4xn,
85         nbnxn_kernel_ElecRF_VdwLJFSw_F_4xn,
86         nbnxn_kernel_ElecRF_VdwLJPSw_F_4xn,
87         nbnxn_kernel_ElecRF_VdwLJEwCombGeom_F_4xn,
88     },
89     {
90         nbnxn_kernel_ElecQSTab_VdwLJCombGeom_F_4xn,
91         nbnxn_kernel_ElecQSTab_VdwLJCombLB_F_4xn,
92         nbnxn_kernel_ElecQSTab_VdwLJ_F_4xn,
93         nbnxn_kernel_ElecQSTab_VdwLJFSw_F_4xn,
94         nbnxn_kernel_ElecQSTab_VdwLJPSw_F_4xn,
95         nbnxn_kernel_ElecQSTab_VdwLJEwCombGeom_F_4xn,
96     },
97     {
98         nbnxn_kernel_ElecQSTabTwinCut_VdwLJCombGeom_F_4xn,
99         nbnxn_kernel_ElecQSTabTwinCut_VdwLJCombLB_F_4xn,
100         nbnxn_kernel_ElecQSTabTwinCut_VdwLJ_F_4xn,
101         nbnxn_kernel_ElecQSTabTwinCut_VdwLJFSw_F_4xn,
102         nbnxn_kernel_ElecQSTabTwinCut_VdwLJPSw_F_4xn,
103         nbnxn_kernel_ElecQSTabTwinCut_VdwLJEwCombGeom_F_4xn,
104     },
105     {
106         nbnxn_kernel_ElecEw_VdwLJCombGeom_F_4xn,
107         nbnxn_kernel_ElecEw_VdwLJCombLB_F_4xn,
108         nbnxn_kernel_ElecEw_VdwLJ_F_4xn,
109         nbnxn_kernel_ElecEw_VdwLJFSw_F_4xn,
110         nbnxn_kernel_ElecEw_VdwLJPSw_F_4xn,
111         nbnxn_kernel_ElecEw_VdwLJEwCombGeom_F_4xn,
112     },
113     {
114         nbnxn_kernel_ElecEwTwinCut_VdwLJCombGeom_F_4xn,
115         nbnxn_kernel_ElecEwTwinCut_VdwLJCombLB_F_4xn,
116         nbnxn_kernel_ElecEwTwinCut_VdwLJ_F_4xn,
117         nbnxn_kernel_ElecEwTwinCut_VdwLJFSw_F_4xn,
118         nbnxn_kernel_ElecEwTwinCut_VdwLJPSw_F_4xn,
119         nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombGeom_F_4xn,
120     },
121 };
122
123 static p_nbk_func_ener p_nbk_ener[coulktNR][vdwktNR] =
124 {
125     {
126         nbnxn_kernel_ElecRF_VdwLJCombGeom_VF_4xn,
127         nbnxn_kernel_ElecRF_VdwLJCombLB_VF_4xn,
128         nbnxn_kernel_ElecRF_VdwLJ_VF_4xn,
129         nbnxn_kernel_ElecRF_VdwLJFSw_VF_4xn,
130         nbnxn_kernel_ElecRF_VdwLJPSw_VF_4xn,
131         nbnxn_kernel_ElecRF_VdwLJEwCombGeom_VF_4xn,
132     },
133     {
134         nbnxn_kernel_ElecQSTab_VdwLJCombGeom_VF_4xn,
135         nbnxn_kernel_ElecQSTab_VdwLJCombLB_VF_4xn,
136         nbnxn_kernel_ElecQSTab_VdwLJ_VF_4xn,
137         nbnxn_kernel_ElecQSTab_VdwLJFSw_VF_4xn,
138         nbnxn_kernel_ElecQSTab_VdwLJPSw_VF_4xn,
139         nbnxn_kernel_ElecQSTab_VdwLJEwCombGeom_VF_4xn,
140     },
141     {
142         nbnxn_kernel_ElecQSTabTwinCut_VdwLJCombGeom_VF_4xn,
143         nbnxn_kernel_ElecQSTabTwinCut_VdwLJCombLB_VF_4xn,
144         nbnxn_kernel_ElecQSTabTwinCut_VdwLJ_VF_4xn,
145         nbnxn_kernel_ElecQSTabTwinCut_VdwLJFSw_VF_4xn,
146         nbnxn_kernel_ElecQSTabTwinCut_VdwLJPSw_VF_4xn,
147         nbnxn_kernel_ElecQSTabTwinCut_VdwLJEwCombGeom_VF_4xn,
148     },
149     {
150         nbnxn_kernel_ElecEw_VdwLJCombGeom_VF_4xn,
151         nbnxn_kernel_ElecEw_VdwLJCombLB_VF_4xn,
152         nbnxn_kernel_ElecEw_VdwLJ_VF_4xn,
153         nbnxn_kernel_ElecEw_VdwLJFSw_VF_4xn,
154         nbnxn_kernel_ElecEw_VdwLJPSw_VF_4xn,
155         nbnxn_kernel_ElecEw_VdwLJEwCombGeom_VF_4xn,
156     },
157     {
158         nbnxn_kernel_ElecEwTwinCut_VdwLJCombGeom_VF_4xn,
159         nbnxn_kernel_ElecEwTwinCut_VdwLJCombLB_VF_4xn,
160         nbnxn_kernel_ElecEwTwinCut_VdwLJ_VF_4xn,
161         nbnxn_kernel_ElecEwTwinCut_VdwLJFSw_VF_4xn,
162         nbnxn_kernel_ElecEwTwinCut_VdwLJPSw_VF_4xn,
163         nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombGeom_VF_4xn,
164     },
165 };
166
167 static p_nbk_func_ener p_nbk_energrp[coulktNR][vdwktNR] =
168 {
169     {
170         nbnxn_kernel_ElecRF_VdwLJCombGeom_VgrpF_4xn,
171         nbnxn_kernel_ElecRF_VdwLJCombLB_VgrpF_4xn,
172         nbnxn_kernel_ElecRF_VdwLJ_VgrpF_4xn,
173         nbnxn_kernel_ElecRF_VdwLJFSw_VgrpF_4xn,
174         nbnxn_kernel_ElecRF_VdwLJPSw_VgrpF_4xn,
175         nbnxn_kernel_ElecRF_VdwLJEwCombGeom_VgrpF_4xn,
176     },
177     {
178         nbnxn_kernel_ElecQSTab_VdwLJCombGeom_VgrpF_4xn,
179         nbnxn_kernel_ElecQSTab_VdwLJCombLB_VgrpF_4xn,
180         nbnxn_kernel_ElecQSTab_VdwLJ_VgrpF_4xn,
181         nbnxn_kernel_ElecQSTab_VdwLJFSw_VgrpF_4xn,
182         nbnxn_kernel_ElecQSTab_VdwLJPSw_VgrpF_4xn,
183         nbnxn_kernel_ElecQSTab_VdwLJEwCombGeom_VgrpF_4xn,
184     },
185     {
186         nbnxn_kernel_ElecQSTabTwinCut_VdwLJCombGeom_VgrpF_4xn,
187         nbnxn_kernel_ElecQSTabTwinCut_VdwLJCombLB_VgrpF_4xn,
188         nbnxn_kernel_ElecQSTabTwinCut_VdwLJ_VgrpF_4xn,
189         nbnxn_kernel_ElecQSTabTwinCut_VdwLJFSw_VgrpF_4xn,
190         nbnxn_kernel_ElecQSTabTwinCut_VdwLJPSw_VgrpF_4xn,
191         nbnxn_kernel_ElecQSTabTwinCut_VdwLJEwCombGeom_VgrpF_4xn,
192     },
193     {
194         nbnxn_kernel_ElecEw_VdwLJCombGeom_VgrpF_4xn,
195         nbnxn_kernel_ElecEw_VdwLJCombLB_VgrpF_4xn,
196         nbnxn_kernel_ElecEw_VdwLJ_VgrpF_4xn,
197         nbnxn_kernel_ElecEw_VdwLJFSw_VgrpF_4xn,
198         nbnxn_kernel_ElecEw_VdwLJPSw_VgrpF_4xn,
199         nbnxn_kernel_ElecEw_VdwLJEwCombGeom_VgrpF_4xn,
200     },
201     {
202         nbnxn_kernel_ElecEwTwinCut_VdwLJCombGeom_VgrpF_4xn,
203         nbnxn_kernel_ElecEwTwinCut_VdwLJCombLB_VgrpF_4xn,
204         nbnxn_kernel_ElecEwTwinCut_VdwLJ_VgrpF_4xn,
205         nbnxn_kernel_ElecEwTwinCut_VdwLJFSw_VgrpF_4xn,
206         nbnxn_kernel_ElecEwTwinCut_VdwLJPSw_VgrpF_4xn,
207         nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombGeom_VgrpF_4xn,
208     },
209 };
210
211
212 static void
213 reduce_group_energies(int ng, int ng_2log,
214                       const real *VSvdw, const real *VSc,
215                       real *Vvdw, real *Vc)
216 {
217     const int unrollj      = GMX_SIMD_REAL_WIDTH/GMX_SIMD_J_UNROLL_SIZE;
218     const int unrollj_half = unrollj/2;
219     int       ng_p2, i, j, j0, j1, c, s;
220
221     ng_p2 = (1<<ng_2log);
222
223     /* The size of the x86 SIMD energy group buffer array is:
224      * ng*ng*ng_p2*unrollj_half*simd_width
225      */
226     for (i = 0; i < ng; i++)
227     {
228         for (j = 0; j < ng; j++)
229         {
230             Vvdw[i*ng+j] = 0;
231             Vc[i*ng+j]   = 0;
232         }
233
234         for (j1 = 0; j1 < ng; j1++)
235         {
236             for (j0 = 0; j0 < ng; j0++)
237             {
238                 c = ((i*ng + j1)*ng_p2 + j0)*unrollj_half*unrollj;
239                 for (s = 0; s < unrollj_half; s++)
240                 {
241                     Vvdw[i*ng+j0] += VSvdw[c+0];
242                     Vvdw[i*ng+j1] += VSvdw[c+1];
243                     Vc  [i*ng+j0] += VSc  [c+0];
244                     Vc  [i*ng+j1] += VSc  [c+1];
245                     c             += unrollj + 2;
246                 }
247             }
248         }
249     }
250 }
251
252 #else /* GMX_NBNXN_SIMD_4XN */
253
254 #include "gmx_fatal.h"
255
256 #endif /* GMX_NBNXN_SIMD_4XN */
257
258 void
259 nbnxn_kernel_simd_4xn(nbnxn_pairlist_set_t      gmx_unused *nbl_list,
260                       const nbnxn_atomdata_t    gmx_unused *nbat,
261                       const interaction_const_t gmx_unused *ic,
262                       int                       gmx_unused  ewald_excl,
263                       rvec                      gmx_unused *shift_vec,
264                       int                       gmx_unused  force_flags,
265                       int                       gmx_unused  clearF,
266                       real                      gmx_unused *fshift,
267                       real                      gmx_unused *Vc,
268                       real                      gmx_unused *Vvdw)
269 #ifdef GMX_NBNXN_SIMD_4XN
270 {
271     int                nnbl;
272     nbnxn_pairlist_t **nbl;
273     int                coulkt, vdwkt = 0;
274     int                nb;
275     int                nthreads gmx_unused;
276
277     nnbl = nbl_list->nnbl;
278     nbl  = nbl_list->nbl;
279
280     if (EEL_RF(ic->eeltype) || ic->eeltype == eelCUT)
281     {
282         coulkt = coulktRF;
283     }
284     else
285     {
286         if (ewald_excl == ewaldexclTable)
287         {
288             if (ic->rcoulomb == ic->rvdw)
289             {
290                 coulkt = coulktTAB;
291             }
292             else
293             {
294                 coulkt = coulktTAB_TWIN;
295             }
296         }
297         else
298         {
299             if (ic->rcoulomb == ic->rvdw)
300             {
301                 coulkt = coulktEWALD;
302             }
303             else
304             {
305                 coulkt = coulktEWALD_TWIN;
306             }
307         }
308     }
309
310     if (ic->vdwtype == evdwCUT)
311     {
312         switch (ic->vdw_modifier)
313         {
314             case eintmodNONE:
315             case eintmodPOTSHIFT:
316                 switch (nbat->comb_rule)
317                 {
318                     case ljcrGEOM: vdwkt = vdwktLJCUT_COMBGEOM; break;
319                     case ljcrLB:   vdwkt = vdwktLJCUT_COMBLB;   break;
320                     case ljcrNONE: vdwkt = vdwktLJCUT_COMBNONE; break;
321                     default:       gmx_incons("Unknown combination rule");
322                 }
323                 break;
324             case eintmodFORCESWITCH:
325                 vdwkt = vdwktLJFORCESWITCH;
326                 break;
327             case eintmodPOTSWITCH:
328                 vdwkt = vdwktLJPOTSWITCH;
329                 break;
330             default:
331                 gmx_incons("Unsupported VdW interaction modifier");
332         }
333     }
334     else if (ic->vdwtype == evdwPME)
335     {
336         if (ic->ljpme_comb_rule == eljpmeLB)
337         {
338             gmx_incons("The nbnxn SIMD kernels don't suport LJ-PME with LB");
339         }
340         vdwkt = vdwktLJEWALDCOMBGEOM;
341     }
342     else
343     {
344         gmx_incons("Unsupported VdW interaction type");
345     }
346
347     nthreads = gmx_omp_nthreads_get(emntNonbonded);
348 #pragma omp parallel for schedule(static) num_threads(nthreads)
349     for (nb = 0; nb < nnbl; nb++)
350     {
351         nbnxn_atomdata_output_t *out;
352         real                    *fshift_p;
353
354         out = &nbat->out[nb];
355
356         if (clearF == enbvClearFYes)
357         {
358             clear_f(nbat, nb, out->f);
359         }
360
361         if ((force_flags & GMX_FORCE_VIRIAL) && nnbl == 1)
362         {
363             fshift_p = fshift;
364         }
365         else
366         {
367             fshift_p = out->fshift;
368
369             if (clearF == enbvClearFYes)
370             {
371                 clear_fshift(fshift_p);
372             }
373         }
374
375         if (!(force_flags & GMX_FORCE_ENERGY))
376         {
377             /* Don't calculate energies */
378             p_nbk_noener[coulkt][vdwkt](nbl[nb], nbat,
379                                         ic,
380                                         shift_vec,
381                                         out->f,
382                                         fshift_p);
383         }
384         else if (out->nV == 1)
385         {
386             /* No energy groups */
387             out->Vvdw[0] = 0;
388             out->Vc[0]   = 0;
389
390             p_nbk_ener[coulkt][vdwkt](nbl[nb], nbat,
391                                       ic,
392                                       shift_vec,
393                                       out->f,
394                                       fshift_p,
395                                       out->Vvdw,
396                                       out->Vc);
397         }
398         else
399         {
400             /* Calculate energy group contributions */
401             int i;
402
403             for (i = 0; i < out->nVS; i++)
404             {
405                 out->VSvdw[i] = 0;
406             }
407             for (i = 0; i < out->nVS; i++)
408             {
409                 out->VSc[i] = 0;
410             }
411
412             p_nbk_energrp[coulkt][vdwkt](nbl[nb], nbat,
413                                          ic,
414                                          shift_vec,
415                                          out->f,
416                                          fshift_p,
417                                          out->VSvdw,
418                                          out->VSc);
419
420             reduce_group_energies(nbat->nenergrp, nbat->neg_2log,
421                                   out->VSvdw, out->VSc,
422                                   out->Vvdw, out->Vc);
423         }
424     }
425
426     if (force_flags & GMX_FORCE_ENERGY)
427     {
428         reduce_energies_over_lists(nbat, nnbl, Vvdw, Vc);
429     }
430 }
431 #else
432 {
433     gmx_incons("nbnxn_kernel_simd_4xn called when such kernels "
434                " are not enabled.");
435 }
436 #endif
437 #undef GMX_SIMD_J_UNROLL_SIZE