0565d3faf0b11854e38d35974d713939415b39e3
[alexxy/gromacs.git] / src / gromacs / nbnxm / nbnxm_setup.cpp
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2019, 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 /*! \internal \file
36  * \brief Common functions for the different NBNXN GPU implementations.
37  *
38  * \author Berk Hess <hess@kth.se>
39  *
40  * \ingroup module_nbnxm
41  */
42
43 #include "gmxpre.h"
44
45 #include "gromacs/domdec/domdec.h"
46 #include "gromacs/domdec/domdec_struct.h"
47 #include "gromacs/hardware/hw_info.h"
48 #include "gromacs/mdlib/gmx_omp_nthreads.h"
49 #include "gromacs/mdtypes/commrec.h"
50 #include "gromacs/mdtypes/forcerec.h"
51 #include "gromacs/mdtypes/inputrec.h"
52 #include "gromacs/nbnxm/atomdata.h"
53 #include "gromacs/nbnxm/gpu_data_mgmt.h"
54 #include "gromacs/nbnxm/nbnxm.h"
55 #include "gromacs/nbnxm/nbnxm_geometry.h"
56 #include "gromacs/nbnxm/nbnxm_simd.h"
57 #include "gromacs/nbnxm/pairlist_tuning.h"
58 #include "gromacs/simd/simd.h"
59 #include "gromacs/utility/fatalerror.h"
60 #include "gromacs/utility/logger.h"
61
62 #include "grid.h"
63 #include "internal.h"
64
65 namespace Nbnxm
66 {
67
68 /*! \brief Returns whether CPU SIMD support exists for the given inputrec
69  *
70  * If the return value is FALSE and fplog/cr != NULL, prints a fallback
71  * message to fplog/stderr.
72  */
73 static gmx_bool nbnxn_simd_supported(const gmx::MDLogger &mdlog,
74                                      const t_inputrec    *ir)
75 {
76     if (ir->vdwtype == evdwPME && ir->ljpme_combination_rule == eljpmeLB)
77     {
78         /* LJ PME with LB combination rule does 7 mesh operations.
79          * This so slow that we don't compile SIMD non-bonded kernels
80          * for that. */
81         GMX_LOG(mdlog.warning).asParagraph().appendText("LJ-PME with Lorentz-Berthelot is not supported with SIMD kernels, falling back to plain C kernels");
82         return FALSE;
83     }
84
85     return TRUE;
86 }
87
88 /*! \brief Returns the most suitable CPU kernel type and Ewald handling */
89 static void pick_nbnxn_kernel_cpu(const t_inputrec gmx_unused    *ir,
90                                   int                            *kernel_type,
91                                   int                            *ewald_excl,
92                                   const gmx_hw_info_t gmx_unused &hardwareInfo)
93 {
94     *kernel_type = nbnxnk4x4_PlainC;
95     *ewald_excl  = ewaldexclTable;
96
97 #if GMX_SIMD
98     {
99 #ifdef GMX_NBNXN_SIMD_4XN
100         *kernel_type = nbnxnk4xN_SIMD_4xN;
101 #endif
102 #ifdef GMX_NBNXN_SIMD_2XNN
103         *kernel_type = nbnxnk4xN_SIMD_2xNN;
104 #endif
105
106 #if defined GMX_NBNXN_SIMD_2XNN && defined GMX_NBNXN_SIMD_4XN
107         /* We need to choose if we want 2x(N+N) or 4xN kernels.
108          * This is based on the SIMD acceleration choice and CPU information
109          * detected at runtime.
110          *
111          * 4xN calculates more (zero) interactions, but has less pair-search
112          * work and much better kernel instruction scheduling.
113          *
114          * Up till now we have only seen that on Intel Sandy/Ivy Bridge,
115          * which doesn't have FMA, both the analytical and tabulated Ewald
116          * kernels have similar pair rates for 4x8 and 2x(4+4), so we choose
117          * 2x(4+4) because it results in significantly fewer pairs.
118          * For RF, the raw pair rate of the 4x8 kernel is higher than 2x(4+4),
119          * 10% with HT, 50% without HT. As we currently don't detect the actual
120          * use of HT, use 4x8 to avoid a potential performance hit.
121          * On Intel Haswell 4x8 is always faster.
122          */
123         *kernel_type = nbnxnk4xN_SIMD_4xN;
124
125 #if !GMX_SIMD_HAVE_FMA
126         if (EEL_PME_EWALD(ir->coulombtype) ||
127             EVDW_PME(ir->vdwtype))
128         {
129             /* We have Ewald kernels without FMA (Intel Sandy/Ivy Bridge).
130              * There are enough instructions to make 2x(4+4) efficient.
131              */
132             *kernel_type = nbnxnk4xN_SIMD_2xNN;
133         }
134 #endif
135         if (hardwareInfo.haveAmdZenCpu)
136         {
137             /* One 256-bit FMA per cycle makes 2xNN faster */
138             *kernel_type = nbnxnk4xN_SIMD_2xNN;
139         }
140 #endif      /* GMX_NBNXN_SIMD_2XNN && GMX_NBNXN_SIMD_4XN */
141
142
143         if (getenv("GMX_NBNXN_SIMD_4XN") != nullptr)
144         {
145 #ifdef GMX_NBNXN_SIMD_4XN
146             *kernel_type = nbnxnk4xN_SIMD_4xN;
147 #else
148             gmx_fatal(FARGS, "SIMD 4xN kernels requested, but GROMACS has been compiled without support for these kernels");
149 #endif
150         }
151         if (getenv("GMX_NBNXN_SIMD_2XNN") != nullptr)
152         {
153 #ifdef GMX_NBNXN_SIMD_2XNN
154             *kernel_type = nbnxnk4xN_SIMD_2xNN;
155 #else
156             gmx_fatal(FARGS, "SIMD 2x(N+N) kernels requested, but GROMACS has been compiled without support for these kernels");
157 #endif
158         }
159
160         /* Analytical Ewald exclusion correction is only an option in
161          * the SIMD kernel.
162          * Since table lookup's don't parallelize with SIMD, analytical
163          * will probably always be faster for a SIMD width of 8 or more.
164          * With FMA analytical is sometimes faster for a width if 4 as well.
165          * In single precision, this is faster on Bulldozer.
166          */
167 #if GMX_SIMD_REAL_WIDTH >= 8 || \
168         (GMX_SIMD_REAL_WIDTH >= 4 && GMX_SIMD_HAVE_FMA && !GMX_DOUBLE)
169         /* On AMD Zen, tabulated Ewald kernels are faster on all 4 combinations
170          * of single or double precision and 128 or 256-bit AVX2.
171          */
172         if (!hardwareInfo.haveAmdZenCpu)
173         {
174             *ewald_excl = ewaldexclAnalytical;
175         }
176 #endif
177         if (getenv("GMX_NBNXN_EWALD_TABLE") != nullptr)
178         {
179             *ewald_excl = ewaldexclTable;
180         }
181         if (getenv("GMX_NBNXN_EWALD_ANALYTICAL") != nullptr)
182         {
183             *ewald_excl = ewaldexclAnalytical;
184         }
185
186     }
187 #endif  // GMX_SIMD
188 }
189
190 const char *lookup_kernel_name(int kernel_type)
191 {
192     const char *returnvalue = nullptr;
193     switch (kernel_type)
194     {
195         case nbnxnkNotSet:
196             returnvalue = "not set";
197             break;
198         case nbnxnk4x4_PlainC:
199             returnvalue = "plain C";
200             break;
201         case nbnxnk4xN_SIMD_4xN:
202         case nbnxnk4xN_SIMD_2xNN:
203 #if GMX_SIMD
204             returnvalue = "SIMD";
205 #else  // GMX_SIMD
206             returnvalue = "not available";
207 #endif // GMX_SIMD
208             break;
209         case nbnxnk8x8x8_GPU: returnvalue    = "GPU"; break;
210         case nbnxnk8x8x8_PlainC: returnvalue = "plain C"; break;
211
212         case nbnxnkNR:
213         default:
214             gmx_fatal(FARGS, "Illegal kernel type selected");
215     }
216     return returnvalue;
217 };
218
219 /*! \brief Returns the most suitable kernel type and Ewald handling */
220 static void pick_nbnxn_kernel(const gmx::MDLogger &mdlog,
221                               gmx_bool             use_simd_kernels,
222                               const gmx_hw_info_t &hardwareInfo,
223                               gmx_bool             bUseGPU,
224                               EmulateGpuNonbonded  emulateGpu,
225                               const t_inputrec    *ir,
226                               int                 *kernel_type,
227                               int                 *ewald_excl,
228                               gmx_bool             bDoNonbonded)
229 {
230     GMX_RELEASE_ASSERT(kernel_type, "Need a valid kernel_type pointer");
231
232     *kernel_type = nbnxnkNotSet;
233     *ewald_excl  = ewaldexclTable;
234
235     if (emulateGpu == EmulateGpuNonbonded::Yes)
236     {
237         *kernel_type = nbnxnk8x8x8_PlainC;
238
239         if (bDoNonbonded)
240         {
241             GMX_LOG(mdlog.warning).asParagraph().appendText("Emulating a GPU run on the CPU (slow)");
242         }
243     }
244     else if (bUseGPU)
245     {
246         *kernel_type = nbnxnk8x8x8_GPU;
247     }
248
249     if (*kernel_type == nbnxnkNotSet)
250     {
251         if (use_simd_kernels &&
252             nbnxn_simd_supported(mdlog, ir))
253         {
254             pick_nbnxn_kernel_cpu(ir, kernel_type, ewald_excl, hardwareInfo);
255         }
256         else
257         {
258             *kernel_type = nbnxnk4x4_PlainC;
259         }
260     }
261
262     if (bDoNonbonded)
263     {
264         GMX_LOG(mdlog.info).asParagraph().appendTextFormatted(
265                 "Using %s %dx%d nonbonded short-range kernels",
266                 lookup_kernel_name(*kernel_type),
267                 nbnxn_kernel_to_cluster_i_size(*kernel_type),
268                 nbnxn_kernel_to_cluster_j_size(*kernel_type));
269
270         if (nbnxnk4x4_PlainC == *kernel_type ||
271             nbnxnk8x8x8_PlainC == *kernel_type)
272         {
273             GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted(
274                     "WARNING: Using the slow %s kernels. This should\n"
275                     "not happen during routine usage on supported platforms.",
276                     lookup_kernel_name(*kernel_type));
277         }
278     }
279 }
280
281 void init_nb_verlet(const gmx::MDLogger     &mdlog,
282                     nonbonded_verlet_t     **nb_verlet,
283                     gmx_bool                 bFEP_NonBonded,
284                     const t_inputrec        *ir,
285                     const t_forcerec        *fr,
286                     const t_commrec         *cr,
287                     const gmx_hw_info_t     &hardwareInfo,
288                     const gmx_device_info_t *deviceInfo,
289                     const gmx_mtop_t        *mtop,
290                     matrix                   box)
291 {
292     nonbonded_verlet_t *nbv;
293     char               *env;
294
295     nbv = new nonbonded_verlet_t();
296
297     nbv->emulateGpu = ((getenv("GMX_EMULATE_GPU") != nullptr) ? EmulateGpuNonbonded::Yes : EmulateGpuNonbonded::No);
298     nbv->bUseGPU    = deviceInfo != nullptr;
299
300     GMX_RELEASE_ASSERT(!(nbv->emulateGpu == EmulateGpuNonbonded::Yes && nbv->bUseGPU), "When GPU emulation is active, there cannot be a GPU assignment");
301
302     nbv->nbs             = nullptr;
303     nbv->min_ci_balanced = 0;
304
305     nbv->ngrp = (DOMAINDECOMP(cr) ? 2 : 1);
306     for (int i = 0; i < nbv->ngrp; i++)
307     {
308         nbv->grp[i].nbl_lists.nnbl = 0;
309         nbv->grp[i].kernel_type    = nbnxnkNotSet;
310
311         if (i == 0) /* local */
312         {
313             pick_nbnxn_kernel(mdlog, fr->use_simd_kernels, hardwareInfo,
314                               nbv->bUseGPU, nbv->emulateGpu, ir,
315                               &nbv->grp[i].kernel_type,
316                               &nbv->grp[i].ewald_excl,
317                               fr->bNonbonded);
318         }
319         else /* non-local */
320         {
321             /* Use the same kernel for local and non-local interactions */
322             nbv->grp[i].kernel_type = nbv->grp[0].kernel_type;
323             nbv->grp[i].ewald_excl  = nbv->grp[0].ewald_excl;
324         }
325     }
326
327     nbv->listParams = std::make_unique<NbnxnListParameters>(ir->rlist);
328     setupDynamicPairlistPruning(mdlog, ir, mtop, box, nbv->grp[0].kernel_type, fr->ic,
329                                 nbv->listParams.get());
330
331     nbv->nbs = std::make_unique<nbnxn_search>(DOMAINDECOMP(cr) ? &cr->dd->nc : nullptr,
332                                               DOMAINDECOMP(cr) ? domdec_zones(cr->dd) : nullptr,
333                                               bFEP_NonBonded,
334                                               gmx_omp_nthreads_get(emntPairsearch));
335
336     for (int i = 0; i < nbv->ngrp; i++)
337     {
338         nbnxn_init_pairlist_set(&nbv->grp[i].nbl_lists,
339                                 nbnxn_kernel_pairlist_simple(nbv->grp[i].kernel_type),
340                                 /* 8x8x8 "non-simple" lists are ATM always combined */
341                                 !nbnxn_kernel_pairlist_simple(nbv->grp[i].kernel_type));
342     }
343
344     int      enbnxninitcombrule;
345     if (fr->ic->vdwtype == evdwCUT &&
346         (fr->ic->vdw_modifier == eintmodNONE ||
347          fr->ic->vdw_modifier == eintmodPOTSHIFT) &&
348         getenv("GMX_NO_LJ_COMB_RULE") == nullptr)
349     {
350         /* Plain LJ cut-off: we can optimize with combination rules */
351         enbnxninitcombrule = enbnxninitcombruleDETECT;
352     }
353     else if (fr->ic->vdwtype == evdwPME)
354     {
355         /* LJ-PME: we need to use a combination rule for the grid */
356         if (fr->ljpme_combination_rule == eljpmeGEOM)
357         {
358             enbnxninitcombrule = enbnxninitcombruleGEOM;
359         }
360         else
361         {
362             enbnxninitcombrule = enbnxninitcombruleLB;
363         }
364     }
365     else
366     {
367         /* We use a full combination matrix: no rule required */
368         enbnxninitcombrule = enbnxninitcombruleNONE;
369     }
370
371     nbv->nbat = new nbnxn_atomdata_t(nbv->bUseGPU ? gmx::PinningPolicy::PinnedIfSupported : gmx::PinningPolicy::CannotBePinned);
372     int mimimumNumEnergyGroupNonbonded = ir->opts.ngener;
373     if (ir->opts.ngener - ir->nwall == 1)
374     {
375         /* We have only one non-wall energy group, we do not need energy group
376          * support in the non-bondeds kernels, since all non-bonded energy
377          * contributions go to the first element of the energy group matrix.
378          */
379         mimimumNumEnergyGroupNonbonded = 1;
380     }
381     bool bSimpleList = nbnxn_kernel_pairlist_simple(nbv->grp[0].kernel_type);
382     nbnxn_atomdata_init(mdlog,
383                         nbv->nbat,
384                         nbv->grp[0].kernel_type,
385                         enbnxninitcombrule,
386                         fr->ntype, fr->nbfp,
387                         mimimumNumEnergyGroupNonbonded,
388                         bSimpleList ? gmx_omp_nthreads_get(emntNonbonded) : 1);
389
390     if (nbv->bUseGPU)
391     {
392         /* init the NxN GPU data; the last argument tells whether we'll have
393          * both local and non-local NB calculation on GPU */
394         gpu_init(&nbv->gpu_nbv,
395                  deviceInfo,
396                  fr->ic,
397                  nbv->listParams.get(),
398                  nbv->nbat,
399                  cr->nodeid,
400                  (nbv->ngrp > 1));
401
402         if ((env = getenv("GMX_NB_MIN_CI")) != nullptr)
403         {
404             char *end;
405
406             nbv->min_ci_balanced = strtol(env, &end, 10);
407             if (!end || (*end != 0) || nbv->min_ci_balanced < 0)
408             {
409                 gmx_fatal(FARGS, "Invalid value passed in GMX_NB_MIN_CI=%s, non-negative integer required", env);
410             }
411
412             if (debug)
413             {
414                 fprintf(debug, "Neighbor-list balancing parameter: %d (passed as env. var.)\n",
415                         nbv->min_ci_balanced);
416             }
417         }
418         else
419         {
420             nbv->min_ci_balanced = gpu_min_ci_balanced(nbv->gpu_nbv);
421             if (debug)
422             {
423                 fprintf(debug, "Neighbor-list balancing parameter: %d (auto-adjusted to the number of GPU multi-processors)\n",
424                         nbv->min_ci_balanced);
425             }
426         }
427
428     }
429
430     *nb_verlet = nbv;
431 }
432
433 } // namespace Nbnxm