src/gromacs/nbnxm/nbnxm_setup.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   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 <memory>
  46
  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"
  64
  65 #include "freeenergydispatch.h"
  66 #include "grid.h"
  67 #include "nbnxm_geometry.h"
  68 #include "nbnxm_simd.h"
  69 #include "pairlist.h"
  70 #include "pairlistset.h"
  71 #include "pairlistsets.h"
  72 #include "pairsearch.h"
  73
  74 namespace Nbnxm
  75 {
  76
  77 /*! \brief Resources that can be used to execute non-bonded kernels on */
  78 enum class NonbondedResource : int
  79 {
  80     Cpu,
  81     Gpu,
  82     EmulateGpu
  83 };
  84
  85 /*! \brief Returns whether CPU SIMD support exists for the given inputrec
  86  *
  87  * If the return value is FALSE and fplog/cr != NULL, prints a fallback
  88  * message to fplog/stderr.
  89  */
  90 static bool nbnxn_simd_supported(const gmx::MDLogger& mdlog, const t_inputrec& inputrec)
  91 {
  92     if (inputrec.vdwtype == VanDerWaalsType::Pme && inputrec.ljpme_combination_rule == LongRangeVdW::LB)
  93     {
  94         /* LJ PME with LB combination rule does 7 mesh operations.
  95          * This so slow that we don't compile SIMD non-bonded kernels
  96          * for that. */
  97         GMX_LOG(mdlog.warning)
  98                 .asParagraph()
  99                 .appendText(
 100                         "LJ-PME with Lorentz-Berthelot is not supported with SIMD kernels, falling "
 101                         "back to plain C kernels");
 102         return FALSE;
 103     }
 104
 105     return TRUE;
 106 }
 107
 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)
 111 {
 112     KernelSetup kernelSetup;
 113
 114     if (!GMX_SIMD)
 115     {
 116         kernelSetup.kernelType         = KernelType::Cpu4x4_PlainC;
 117         kernelSetup.ewaldExclusionType = EwaldExclusionType::Table;
 118     }
 119     else
 120     {
 121 #ifdef GMX_NBNXN_SIMD_4XN
 122         kernelSetup.kernelType = KernelType::Cpu4xN_Simd_4xN;
 123 #endif
 124 #ifdef GMX_NBNXN_SIMD_2XNN
 125         kernelSetup.kernelType = KernelType::Cpu4xN_Simd_2xNN;
 126 #endif
 127
 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.
 132          *
 133          * 4xN calculates more (zero) interactions, but has less pair-search
 134          * work and much better kernel instruction scheduling.
 135          *
 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.
 144          */
 145         kernelSetup.kernelType = KernelType::Cpu4xN_Simd_4xN;
 146
 147         if (!GMX_SIMD_HAVE_FMA && (EEL_PME_EWALD(inputrec.coulombtype) || EVDW_PME(inputrec.vdwtype)))
 148         {
 149             /* We have Ewald kernels without FMA (Intel Sandy/Ivy Bridge).
 150              * There are enough instructions to make 2x(4+4) efficient.
 151              */
 152             kernelSetup.kernelType = KernelType::Cpu4xN_Simd_2xNN;
 153         }
 154
 155         if (hardwareInfo.haveAmdZen1Cpu)
 156         {
 157             /* One 256-bit FMA per cycle makes 2xNN faster */
 158             kernelSetup.kernelType = KernelType::Cpu4xN_Simd_2xNN;
 159         }
 160 #endif /* GMX_NBNXN_SIMD_2XNN && GMX_NBNXN_SIMD_4XN */
 161
 162
 163         if (getenv("GMX_NBNXN_SIMD_4XN") != nullptr)
 164         {
 165 #ifdef GMX_NBNXN_SIMD_4XN
 166             kernelSetup.kernelType = KernelType::Cpu4xN_Simd_4xN;
 167 #else
 168             gmx_fatal(FARGS,
 169                       "SIMD 4xN kernels requested, but GROMACS has been compiled without support "
 170                       "for these kernels");
 171 #endif
 172         }
 173         if (getenv("GMX_NBNXN_SIMD_2XNN") != nullptr)
 174         {
 175 #ifdef GMX_NBNXN_SIMD_2XNN
 176             kernelSetup.kernelType = KernelType::Cpu4xN_Simd_2xNN;
 177 #else
 178             gmx_fatal(FARGS,
 179                       "SIMD 2x(N+N) kernels requested, but GROMACS has been compiled without "
 180                       "support for these kernels");
 181 #endif
 182         }
 183
 184         /* Analytical Ewald exclusion correction is only an option in
 185          * the SIMD kernel.
 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.
 192          */
 193         MSVC_DIAGNOSTIC_IGNORE(6285) // Always zero because compile time constant
 194         if (
 195 #if GMX_SIMD
 196                 (GMX_SIMD_REAL_WIDTH >= 8 || (GMX_SIMD_REAL_WIDTH >= 4 && GMX_SIMD_HAVE_FMA && !GMX_DOUBLE)) &&
 197 #endif
 198                 !hardwareInfo.haveAmdZen1Cpu)
 199         {
 200             kernelSetup.ewaldExclusionType = EwaldExclusionType::Analytical;
 201         }
 202         MSVC_DIAGNOSTIC_RESET
 203         else { kernelSetup.ewaldExclusionType = EwaldExclusionType::Table; }
 204         if (getenv("GMX_NBNXN_EWALD_TABLE") != nullptr)
 205         {
 206             kernelSetup.ewaldExclusionType = EwaldExclusionType::Table;
 207         }
 208         if (getenv("GMX_NBNXN_EWALD_ANALYTICAL") != nullptr)
 209         {
 210             kernelSetup.ewaldExclusionType = EwaldExclusionType::Analytical;
 211         }
 212     }
 213
 214     return kernelSetup;
 215 }
 216
 217 const char* lookup_kernel_name(const KernelType kernelType)
 218 {
 219     switch (kernelType)
 220     {
 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:
 225 #if GMX_SIMD
 226             return "SIMD";
 227 #else  // GMX_SIMD
 228             return "not available";
 229 #endif // GMX_SIMD
 230         case KernelType::Gpu8x8x8: return "GPU";
 231         case KernelType::Cpu8x8x8_PlainC: return "plain C";
 232
 233         default: gmx_fatal(FARGS, "Illegal kernel type selected");
 234     }
 235 };
 236
 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)
 243 {
 244     KernelSetup kernelSetup;
 245
 246     if (nonbondedResource == NonbondedResource::EmulateGpu)
 247     {
 248         kernelSetup.kernelType         = KernelType::Cpu8x8x8_PlainC;
 249         kernelSetup.ewaldExclusionType = EwaldExclusionType::DecidedByGpuModule;
 250
 251         GMX_LOG(mdlog.warning).asParagraph().appendText("Emulating a GPU run on the CPU (slow)");
 252     }
 253     else if (nonbondedResource == NonbondedResource::Gpu)
 254     {
 255         kernelSetup.kernelType         = KernelType::Gpu8x8x8;
 256         kernelSetup.ewaldExclusionType = EwaldExclusionType::DecidedByGpuModule;
 257     }
 258     else
 259     {
 260         if (use_simd_kernels && nbnxn_simd_supported(mdlog, inputrec))
 261         {
 262             kernelSetup = pick_nbnxn_kernel_cpu(inputrec, hardwareInfo);
 263         }
 264         else
 265         {
 266             kernelSetup.kernelType         = KernelType::Cpu4x4_PlainC;
 267             kernelSetup.ewaldExclusionType = EwaldExclusionType::Analytical;
 268         }
 269     }
 270
 271     GMX_LOG(mdlog.info)
 272             .asParagraph()
 273             .appendTextFormatted("Using %s %dx%d nonbonded short-range kernels",
 274                                  lookup_kernel_name(kernelSetup.kernelType),
 275                                  IClusterSizePerKernelType[kernelSetup.kernelType],
 276                                  JClusterSizePerKernelType[kernelSetup.kernelType]);
 277
 278     if (KernelType::Cpu4x4_PlainC == kernelSetup.kernelType
 279         || KernelType::Cpu8x8x8_PlainC == kernelSetup.kernelType)
 280     {
 281         GMX_LOG(mdlog.warning)
 282                 .asParagraph()
 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));
 287     }
 288
 289     GMX_RELEASE_ASSERT(kernelSetup.kernelType != KernelType::NotSet
 290                                && kernelSetup.ewaldExclusionType != EwaldExclusionType::NotSet,
 291                        "All kernel setup parameters should be set here");
 292
 293     return kernelSetup;
 294 }
 295
 296 } // namespace Nbnxm
 297
 298 PairlistSets::PairlistSets(const PairlistParams& pairlistParams,
 299                            const bool            haveMultipleDomains,
 300                            const int             minimumIlistCountForGpuBalancing) :
 301     params_(pairlistParams), minimumIlistCountForGpuBalancing_(minimumIlistCountForGpuBalancing)
 302 {
 303     localSet_ = std::make_unique<PairlistSet>(params_);
 304
 305     if (haveMultipleDomains)
 306     {
 307         nonlocalSet_ = std::make_unique<PairlistSet>(params_);
 308     }
 309 }
 310
 311 namespace Nbnxm
 312 {
 313
 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)
 316 {
 317     if (const char* env = getenv("GMX_NB_MIN_CI"))
 318     {
 319         char* end = nullptr;
 320
 321         int minimumIlistCount = strtol(env, &end, 10);
 322         if (!end || (*end != 0) || minimumIlistCount < 0)
 323         {
 324             gmx_fatal(
 325                     FARGS, "Invalid value passed in GMX_NB_MIN_CI=%s, non-negative integer required", env);
 326         }
 327
 328         if (debug)
 329         {
 330             fprintf(debug, "Neighbor-list balancing parameter: %d (passed as env. var.)\n", minimumIlistCount);
 331         }
 332         return minimumIlistCount;
 333     }
 334     else
 335     {
 336         int minimumIlistCount = gpu_min_ci_balanced(nbnxmGpu);
 337         if (debug)
 338         {
 339             fprintf(debug,
 340                     "Neighbor-list balancing parameter: %d (auto-adjusted to the number of GPU "
 341                     "multi-processors)\n",
 342                     minimumIlistCount);
 343         }
 344         return minimumIlistCount;
 345     }
 346 }
 347
 348 static int getENbnxnInitCombRule(const t_forcerec& forcerec)
 349 {
 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)
 354     {
 355         /* Plain LJ cut-off: we can optimize with combination rules */
 356         return enbnxninitcombruleDETECT;
 357     }
 358     else if (forcerec.ic->vdwtype == VanDerWaalsType::Pme)
 359     {
 360         /* LJ-PME: we need to use a combination rule for the grid */
 361         if (forcerec.ljpme_combination_rule == LongRangeVdW::Geom)
 362         {
 363             return enbnxninitcombruleGEOM;
 364         }
 365         else
 366         {
 367             return enbnxninitcombruleLB;
 368         }
 369     }
 370     else
 371     {
 372         /* We use a full combination matrix: no rule required */
 373         return enbnxninitcombruleNONE;
 374     }
 375 }
 376
 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,
 385                                                    matrix                          box,
 386                                                    gmx_wallcycle*                  wcycle)
 387 {
 388     const bool emulateGpu = (getenv("GMX_EMULATE_GPU") != nullptr);
 389
 390     GMX_RELEASE_ASSERT(!(emulateGpu && useGpuForNonbonded),
 391                        "When GPU emulation is active, there cannot be a GPU assignment");
 392
 393     NonbondedResource nonbondedResource;
 394     if (useGpuForNonbonded)
 395     {
 396         nonbondedResource = NonbondedResource::Gpu;
 397     }
 398     else if (emulateGpu)
 399     {
 400         nonbondedResource = NonbondedResource::EmulateGpu;
 401     }
 402     else
 403     {
 404         nonbondedResource = NonbondedResource::Cpu;
 405     }
 406
 407     Nbnxm::KernelSetup kernelSetup = pick_nbnxn_kernel(
 408             mdlog, forcerec.use_simd_kernels, hardwareInfo, nonbondedResource, inputrec);
 409
 410     const bool haveMultipleDomains = havePPDomainDecomposition(commrec);
 411
 412     bool bFEP_NonBonded = (forcerec.efep != FreeEnergyPerturbationType::No)
 413                           && haveFepPerturbedNBInteractions(mtop);
 414     PairlistParams pairlistParams(
 415             kernelSetup.kernelType, bFEP_NonBonded, inputrec.rlist, haveMultipleDomains);
 416
 417     setupDynamicPairlistPruning(mdlog, inputrec, mtop, box, *forcerec.ic, &pairlistParams);
 418
 419     const int enbnxninitcombrule = getENbnxnInitCombRule(forcerec);
 420
 421     auto pinPolicy = (useGpuForNonbonded ? gmx::PinningPolicy::PinnedIfSupported
 422                                          : gmx::PinningPolicy::CannotBePinned);
 423
 424     int mimimumNumEnergyGroupNonbonded = inputrec.opts.ngener;
 425     if (inputrec.opts.ngener - inputrec.nwall == 1)
 426     {
 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.
 430          */
 431         mimimumNumEnergyGroupNonbonded = 1;
 432     }
 433
 434     auto nbat = std::make_unique<nbnxn_atomdata_t>(
 435             pinPolicy,
 436             mdlog,
 437             kernelSetup.kernelType,
 438             enbnxninitcombrule,
 439             forcerec.ntype,
 440             forcerec.nbfp,
 441             mimimumNumEnergyGroupNonbonded,
 442             (useGpuForNonbonded || emulateGpu) ? 1 : gmx_omp_nthreads_get(ModuleMultiThread::Nonbonded));
 443
 444     NbnxmGpu* gpu_nbv                          = nullptr;
 445     int       minimumIlistCountForGpuBalancing = 0;
 446     if (useGpuForNonbonded)
 447     {
 448         /* init the NxN GPU data; the last argument tells whether we'll have
 449          * both local and non-local NB calculation on GPU */
 450         GMX_RELEASE_ASSERT(
 451                 (deviceStreamManager != nullptr),
 452                 "Device stream manager should be initialized in order to use GPU for non-bonded.");
 453         gpu_nbv = gpu_init(
 454                 *deviceStreamManager, forcerec.ic.get(), pairlistParams, nbat.get(), haveMultipleDomains);
 455
 456         minimumIlistCountForGpuBalancing = getMinimumIlistCountForGpuBalancing(gpu_nbv);
 457     }
 458
 459     auto pairlistSets = std::make_unique<PairlistSets>(
 460             pairlistParams, haveMultipleDomains, minimumIlistCountForGpuBalancing);
 461
 462     auto pairSearch = std::make_unique<PairSearch>(
 463             inputrec.pbcType,
 464             EI_TPI(inputrec.eI),
 465             haveDDAtomOrdering(*commrec) ? &commrec->dd->numCells : nullptr,
 466             haveDDAtomOrdering(*commrec) ? domdec_zones(commrec->dd) : nullptr,
 467             pairlistParams.pairlistType,
 468             bFEP_NonBonded,
 469             gmx_omp_nthreads_get(ModuleMultiThread::Pairsearch),
 470             pinPolicy);
 471
 472     return std::make_unique<nonbonded_verlet_t>(
 473             std::move(pairlistSets), std::move(pairSearch), std::move(nbat), kernelSetup, gpu_nbv, wcycle);
 474 }
 475
 476 } // namespace Nbnxm
 477
 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),
 488     wcycle_(wcycle),
 489     gpu_nbv(gpu_nbv_ptr)
 490 {
 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");
 494
 495     if (pairlistSets_->params().haveFep)
 496     {
 497         freeEnergyDispatch_ = std::make_unique<FreeEnergyDispatch>(nbat->params().nenergrp);
 498     }
 499 }
 500
 501 nonbonded_verlet_t::~nonbonded_verlet_t()
 502 {
 503     Nbnxm::gpu_free(gpu_nbv);
 504 }