2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2019,2020,2021, 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 * \brief Common functions for the different NBNXN GPU implementations.
38 * \author Berk Hess <hess@kth.se>
40 * \ingroup module_nbnxm
47 #include "gromacs/domdec/domdec.h"
48 #include "gromacs/domdec/domdec_struct.h"
49 #include "gromacs/hardware/hw_info.h"
50 #include "gromacs/mdlib/gmx_omp_nthreads.h"
51 #include "gromacs/mdtypes/commrec.h"
52 #include "gromacs/mdtypes/enerdata.h"
53 #include "gromacs/mdtypes/forcerec.h"
54 #include "gromacs/mdtypes/inputrec.h"
55 #include "gromacs/mdtypes/interaction_const.h"
56 #include "gromacs/nbnxm/atomdata.h"
57 #include "gromacs/nbnxm/gpu_data_mgmt.h"
58 #include "gromacs/nbnxm/nbnxm.h"
59 #include "gromacs/nbnxm/pairlist_tuning.h"
60 #include "gromacs/simd/simd.h"
61 #include "gromacs/topology/mtop_util.h"
62 #include "gromacs/utility/fatalerror.h"
63 #include "gromacs/utility/logger.h"
65 #include "freeenergydispatch.h"
67 #include "nbnxm_geometry.h"
68 #include "nbnxm_simd.h"
70 #include "pairlistset.h"
71 #include "pairlistsets.h"
72 #include "pairsearch.h"
77 /*! \brief Resources that can be used to execute non-bonded kernels on */
78 enum class NonbondedResource : int
85 /*! \brief Returns whether CPU SIMD support exists for the given inputrec
87 * If the return value is FALSE and fplog/cr != NULL, prints a fallback
88 * message to fplog/stderr.
90 static bool nbnxn_simd_supported(const gmx::MDLogger& mdlog, const t_inputrec& inputrec)
92 if (inputrec.vdwtype == VanDerWaalsType::Pme && inputrec.ljpme_combination_rule == LongRangeVdW::LB)
94 /* LJ PME with LB combination rule does 7 mesh operations.
95 * This so slow that we don't compile SIMD non-bonded kernels
97 GMX_LOG(mdlog.warning)
100 "LJ-PME with Lorentz-Berthelot is not supported with SIMD kernels, falling "
101 "back to plain C kernels");
108 /*! \brief Returns the most suitable CPU kernel type and Ewald handling */
109 static KernelSetup pick_nbnxn_kernel_cpu(const t_inputrec gmx_unused& inputrec,
110 const gmx_hw_info_t gmx_unused& hardwareInfo)
112 KernelSetup kernelSetup;
116 kernelSetup.kernelType = KernelType::Cpu4x4_PlainC;
117 kernelSetup.ewaldExclusionType = EwaldExclusionType::Table;
121 #ifdef GMX_NBNXN_SIMD_4XN
122 kernelSetup.kernelType = KernelType::Cpu4xN_Simd_4xN;
124 #ifdef GMX_NBNXN_SIMD_2XNN
125 kernelSetup.kernelType = KernelType::Cpu4xN_Simd_2xNN;
128 #if defined GMX_NBNXN_SIMD_2XNN && defined GMX_NBNXN_SIMD_4XN
129 /* We need to choose if we want 2x(N+N) or 4xN kernels.
130 * This is based on the SIMD acceleration choice and CPU information
131 * detected at runtime.
133 * 4xN calculates more (zero) interactions, but has less pair-search
134 * work and much better kernel instruction scheduling.
136 * Up till now we have only seen that on Intel Sandy/Ivy Bridge,
137 * which doesn't have FMA, both the analytical and tabulated Ewald
138 * kernels have similar pair rates for 4x8 and 2x(4+4), so we choose
139 * 2x(4+4) because it results in significantly fewer pairs.
140 * For RF, the raw pair rate of the 4x8 kernel is higher than 2x(4+4),
141 * 10% with HT, 50% without HT. As we currently don't detect the actual
142 * use of HT, use 4x8 to avoid a potential performance hit.
143 * On Intel Haswell 4x8 is always faster.
145 kernelSetup.kernelType = KernelType::Cpu4xN_Simd_4xN;
147 if (!GMX_SIMD_HAVE_FMA && (EEL_PME_EWALD(inputrec.coulombtype) || EVDW_PME(inputrec.vdwtype)))
149 /* We have Ewald kernels without FMA (Intel Sandy/Ivy Bridge).
150 * There are enough instructions to make 2x(4+4) efficient.
152 kernelSetup.kernelType = KernelType::Cpu4xN_Simd_2xNN;
155 if (hardwareInfo.haveAmdZen1Cpu)
157 /* One 256-bit FMA per cycle makes 2xNN faster */
158 kernelSetup.kernelType = KernelType::Cpu4xN_Simd_2xNN;
160 #endif /* GMX_NBNXN_SIMD_2XNN && GMX_NBNXN_SIMD_4XN */
163 if (getenv("GMX_NBNXN_SIMD_4XN") != nullptr)
165 #ifdef GMX_NBNXN_SIMD_4XN
166 kernelSetup.kernelType = KernelType::Cpu4xN_Simd_4xN;
169 "SIMD 4xN kernels requested, but GROMACS has been compiled without support "
170 "for these kernels");
173 if (getenv("GMX_NBNXN_SIMD_2XNN") != nullptr)
175 #ifdef GMX_NBNXN_SIMD_2XNN
176 kernelSetup.kernelType = KernelType::Cpu4xN_Simd_2xNN;
179 "SIMD 2x(N+N) kernels requested, but GROMACS has been compiled without "
180 "support for these kernels");
184 /* Analytical Ewald exclusion correction is only an option in
186 * Since table lookup's don't parallelize with SIMD, analytical
187 * will probably always be faster for a SIMD width of 8 or more.
188 * With FMA analytical is sometimes faster for a width if 4 as well.
189 * In single precision, this is faster on Bulldozer.
190 * On AMD Zen, tabulated Ewald kernels are faster on all 4 combinations
191 * of single or double precision and 128 or 256-bit AVX2.
193 MSVC_DIAGNOSTIC_IGNORE(6285) // Always zero because compile time constant
196 (GMX_SIMD_REAL_WIDTH >= 8 || (GMX_SIMD_REAL_WIDTH >= 4 && GMX_SIMD_HAVE_FMA && !GMX_DOUBLE)) &&
198 !hardwareInfo.haveAmdZen1Cpu)
200 kernelSetup.ewaldExclusionType = EwaldExclusionType::Analytical;
202 MSVC_DIAGNOSTIC_RESET
203 else { kernelSetup.ewaldExclusionType = EwaldExclusionType::Table; }
204 if (getenv("GMX_NBNXN_EWALD_TABLE") != nullptr)
206 kernelSetup.ewaldExclusionType = EwaldExclusionType::Table;
208 if (getenv("GMX_NBNXN_EWALD_ANALYTICAL") != nullptr)
210 kernelSetup.ewaldExclusionType = EwaldExclusionType::Analytical;
217 const char* lookup_kernel_name(const KernelType kernelType)
221 case KernelType::NotSet: return "not set";
222 case KernelType::Cpu4x4_PlainC: return "plain C";
223 case KernelType::Cpu4xN_Simd_4xN:
224 case KernelType::Cpu4xN_Simd_2xNN:
228 return "not available";
230 case KernelType::Gpu8x8x8: return "GPU";
231 case KernelType::Cpu8x8x8_PlainC: return "plain C";
233 default: gmx_fatal(FARGS, "Illegal kernel type selected");
237 /*! \brief Returns the most suitable kernel type and Ewald handling */
238 static KernelSetup pick_nbnxn_kernel(const gmx::MDLogger& mdlog,
239 gmx_bool use_simd_kernels,
240 const gmx_hw_info_t& hardwareInfo,
241 const NonbondedResource& nonbondedResource,
242 const t_inputrec& inputrec)
244 KernelSetup kernelSetup;
246 if (nonbondedResource == NonbondedResource::EmulateGpu)
248 kernelSetup.kernelType = KernelType::Cpu8x8x8_PlainC;
249 kernelSetup.ewaldExclusionType = EwaldExclusionType::DecidedByGpuModule;
251 GMX_LOG(mdlog.warning).asParagraph().appendText("Emulating a GPU run on the CPU (slow)");
253 else if (nonbondedResource == NonbondedResource::Gpu)
255 kernelSetup.kernelType = KernelType::Gpu8x8x8;
256 kernelSetup.ewaldExclusionType = EwaldExclusionType::DecidedByGpuModule;
260 if (use_simd_kernels && nbnxn_simd_supported(mdlog, inputrec))
262 kernelSetup = pick_nbnxn_kernel_cpu(inputrec, hardwareInfo);
266 kernelSetup.kernelType = KernelType::Cpu4x4_PlainC;
267 kernelSetup.ewaldExclusionType = EwaldExclusionType::Analytical;
273 .appendTextFormatted("Using %s %dx%d nonbonded short-range kernels",
274 lookup_kernel_name(kernelSetup.kernelType),
275 IClusterSizePerKernelType[kernelSetup.kernelType],
276 JClusterSizePerKernelType[kernelSetup.kernelType]);
278 if (KernelType::Cpu4x4_PlainC == kernelSetup.kernelType
279 || KernelType::Cpu8x8x8_PlainC == kernelSetup.kernelType)
281 GMX_LOG(mdlog.warning)
283 .appendTextFormatted(
284 "WARNING: Using the slow %s kernels. This should\n"
285 "not happen during routine usage on supported platforms.",
286 lookup_kernel_name(kernelSetup.kernelType));
289 GMX_RELEASE_ASSERT(kernelSetup.kernelType != KernelType::NotSet
290 && kernelSetup.ewaldExclusionType != EwaldExclusionType::NotSet,
291 "All kernel setup parameters should be set here");
298 PairlistSets::PairlistSets(const PairlistParams& pairlistParams,
299 const bool haveMultipleDomains,
300 const int minimumIlistCountForGpuBalancing) :
301 params_(pairlistParams), minimumIlistCountForGpuBalancing_(minimumIlistCountForGpuBalancing)
303 localSet_ = std::make_unique<PairlistSet>(params_);
305 if (haveMultipleDomains)
307 nonlocalSet_ = std::make_unique<PairlistSet>(params_);
314 /*! \brief Gets and returns the minimum i-list count for balacing based on the GPU used or env.var. when set */
315 static int getMinimumIlistCountForGpuBalancing(NbnxmGpu* nbnxmGpu)
317 if (const char* env = getenv("GMX_NB_MIN_CI"))
321 int minimumIlistCount = strtol(env, &end, 10);
322 if (!end || (*end != 0) || minimumIlistCount < 0)
325 FARGS, "Invalid value passed in GMX_NB_MIN_CI=%s, non-negative integer required", env);
330 fprintf(debug, "Neighbor-list balancing parameter: %d (passed as env. var.)\n", minimumIlistCount);
332 return minimumIlistCount;
336 int minimumIlistCount = gpu_min_ci_balanced(nbnxmGpu);
340 "Neighbor-list balancing parameter: %d (auto-adjusted to the number of GPU "
341 "multi-processors)\n",
344 return minimumIlistCount;
348 static int getENbnxnInitCombRule(const t_forcerec& forcerec)
350 if (forcerec.ic->vdwtype == VanDerWaalsType::Cut
351 && (forcerec.ic->vdw_modifier == InteractionModifiers::None
352 || forcerec.ic->vdw_modifier == InteractionModifiers::PotShift)
353 && getenv("GMX_NO_LJ_COMB_RULE") == nullptr)
355 /* Plain LJ cut-off: we can optimize with combination rules */
356 return enbnxninitcombruleDETECT;
358 else if (forcerec.ic->vdwtype == VanDerWaalsType::Pme)
360 /* LJ-PME: we need to use a combination rule for the grid */
361 if (forcerec.ljpme_combination_rule == LongRangeVdW::Geom)
363 return enbnxninitcombruleGEOM;
367 return enbnxninitcombruleLB;
372 /* We use a full combination matrix: no rule required */
373 return enbnxninitcombruleNONE;
377 std::unique_ptr<nonbonded_verlet_t> init_nb_verlet(const gmx::MDLogger& mdlog,
378 const t_inputrec& inputrec,
379 const t_forcerec& forcerec,
380 const t_commrec* commrec,
381 const gmx_hw_info_t& hardwareInfo,
382 bool useGpuForNonbonded,
383 const gmx::DeviceStreamManager* deviceStreamManager,
384 const gmx_mtop_t& mtop,
386 gmx_wallcycle* wcycle)
388 const bool emulateGpu = (getenv("GMX_EMULATE_GPU") != nullptr);
390 GMX_RELEASE_ASSERT(!(emulateGpu && useGpuForNonbonded),
391 "When GPU emulation is active, there cannot be a GPU assignment");
393 NonbondedResource nonbondedResource;
394 if (useGpuForNonbonded)
396 nonbondedResource = NonbondedResource::Gpu;
400 nonbondedResource = NonbondedResource::EmulateGpu;
404 nonbondedResource = NonbondedResource::Cpu;
407 Nbnxm::KernelSetup kernelSetup = pick_nbnxn_kernel(
408 mdlog, forcerec.use_simd_kernels, hardwareInfo, nonbondedResource, inputrec);
410 const bool haveMultipleDomains = havePPDomainDecomposition(commrec);
412 bool bFEP_NonBonded = (forcerec.efep != FreeEnergyPerturbationType::No)
413 && haveFepPerturbedNBInteractions(mtop);
414 PairlistParams pairlistParams(
415 kernelSetup.kernelType, bFEP_NonBonded, inputrec.rlist, haveMultipleDomains);
417 setupDynamicPairlistPruning(mdlog, inputrec, mtop, box, *forcerec.ic, &pairlistParams);
419 const int enbnxninitcombrule = getENbnxnInitCombRule(forcerec);
421 auto pinPolicy = (useGpuForNonbonded ? gmx::PinningPolicy::PinnedIfSupported
422 : gmx::PinningPolicy::CannotBePinned);
424 int mimimumNumEnergyGroupNonbonded = inputrec.opts.ngener;
425 if (inputrec.opts.ngener - inputrec.nwall == 1)
427 /* We have only one non-wall energy group, we do not need energy group
428 * support in the non-bondeds kernels, since all non-bonded energy
429 * contributions go to the first element of the energy group matrix.
431 mimimumNumEnergyGroupNonbonded = 1;
434 auto nbat = std::make_unique<nbnxn_atomdata_t>(
437 kernelSetup.kernelType,
441 mimimumNumEnergyGroupNonbonded,
442 (useGpuForNonbonded || emulateGpu) ? 1 : gmx_omp_nthreads_get(ModuleMultiThread::Nonbonded));
444 NbnxmGpu* gpu_nbv = nullptr;
445 int minimumIlistCountForGpuBalancing = 0;
446 if (useGpuForNonbonded)
448 /* init the NxN GPU data; the last argument tells whether we'll have
449 * both local and non-local NB calculation on GPU */
451 (deviceStreamManager != nullptr),
452 "Device stream manager should be initialized in order to use GPU for non-bonded.");
454 *deviceStreamManager, forcerec.ic.get(), pairlistParams, nbat.get(), haveMultipleDomains);
456 minimumIlistCountForGpuBalancing = getMinimumIlistCountForGpuBalancing(gpu_nbv);
459 auto pairlistSets = std::make_unique<PairlistSets>(
460 pairlistParams, haveMultipleDomains, minimumIlistCountForGpuBalancing);
462 auto pairSearch = std::make_unique<PairSearch>(
465 haveDDAtomOrdering(*commrec) ? &commrec->dd->numCells : nullptr,
466 haveDDAtomOrdering(*commrec) ? domdec_zones(commrec->dd) : nullptr,
467 pairlistParams.pairlistType,
469 gmx_omp_nthreads_get(ModuleMultiThread::Pairsearch),
472 return std::make_unique<nonbonded_verlet_t>(
473 std::move(pairlistSets), std::move(pairSearch), std::move(nbat), kernelSetup, gpu_nbv, wcycle);
478 nonbonded_verlet_t::nonbonded_verlet_t(std::unique_ptr<PairlistSets> pairlistSets,
479 std::unique_ptr<PairSearch> pairSearch,
480 std::unique_ptr<nbnxn_atomdata_t> nbat_in,
481 const Nbnxm::KernelSetup& kernelSetup,
482 NbnxmGpu* gpu_nbv_ptr,
483 gmx_wallcycle* wcycle) :
484 pairlistSets_(std::move(pairlistSets)),
485 pairSearch_(std::move(pairSearch)),
486 nbat(std::move(nbat_in)),
487 kernelSetup_(kernelSetup),
491 GMX_RELEASE_ASSERT(pairlistSets_, "Need valid pairlistSets");
492 GMX_RELEASE_ASSERT(pairSearch_, "Need valid search object");
493 GMX_RELEASE_ASSERT(nbat, "Need valid atomdata object");
495 if (pairlistSets_->params().haveFep)
497 freeEnergyDispatch_ = std::make_unique<FreeEnergyDispatch>(nbat->params().nenergrp);
501 nonbonded_verlet_t::~nonbonded_verlet_t()
503 Nbnxm::gpu_free(gpu_nbv);