src/gromacs/nbnxm/kerneldispatch.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014,2015,2016 by the GROMACS development team.
   5  * Copyright (c) 2017,2018,2019,2020,2021, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36
  37 #include "gmxpre.h"
  38
  39 #include "gromacs/gmxlib/nrnb.h"
  40 #include "gromacs/gmxlib/nonbonded/nb_free_energy.h"
  41 #include "gromacs/gmxlib/nonbonded/nonbonded.h"
  42 #include "gromacs/math/vectypes.h"
  43 #include "gromacs/mdlib/enerdata_utils.h"
  44 #include "gromacs/mdlib/force.h"
  45 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  46 #include "gromacs/mdtypes/enerdata.h"
  47 #include "gromacs/mdtypes/forceoutput.h"
  48 #include "gromacs/mdtypes/forcerec.h"
  49 #include "gromacs/mdtypes/inputrec.h"
  50 #include "gromacs/mdtypes/interaction_const.h"
  51 #include "gromacs/mdtypes/md_enums.h"
  52 #include "gromacs/mdtypes/mdatom.h"
  53 #include "gromacs/mdtypes/nblist.h"
  54 #include "gromacs/mdtypes/simulation_workload.h"
  55 #include "gromacs/nbnxm/gpu_data_mgmt.h"
  56 #include "gromacs/nbnxm/nbnxm.h"
  57 #include "gromacs/simd/simd.h"
  58 #include "gromacs/timing/wallcycle.h"
  59 #include "gromacs/utility/enumerationhelpers.h"
  60 #include "gromacs/utility/fatalerror.h"
  61 #include "gromacs/utility/gmxassert.h"
  62 #include "gromacs/utility/real.h"
  63
  64 #include "kernel_common.h"
  65 #include "nbnxm_gpu.h"
  66 #include "nbnxm_simd.h"
  67 #include "pairlistset.h"
  68 #include "pairlistsets.h"
  69 #include "kernels_reference/kernel_gpu_ref.h"
  70 #define INCLUDE_KERNELFUNCTION_TABLES
  71 #include "kernels_reference/kernel_ref.h"
  72 #ifdef GMX_NBNXN_SIMD_2XNN
  73 #    include "kernels_simd_2xmm/kernels.h"
  74 #endif
  75 #ifdef GMX_NBNXN_SIMD_4XN
  76 #    include "kernels_simd_4xm/kernels.h"
  77 #endif
  78 #undef INCLUDE_FUNCTION_TABLES
  79
  80 /*! \brief Clears the energy group output buffers
  81  *
  82  * \param[in,out] out  nbnxn kernel output struct
  83  */
  84 static void clearGroupEnergies(nbnxn_atomdata_output_t* out)
  85 {
  86     std::fill(out->Vvdw.begin(), out->Vvdw.end(), 0.0_real);
  87     std::fill(out->Vc.begin(), out->Vc.end(), 0.0_real);
  88     std::fill(out->VSvdw.begin(), out->VSvdw.end(), 0.0_real);
  89     std::fill(out->VSc.begin(), out->VSc.end(), 0.0_real);
  90 }
  91
  92 /*! \brief Reduce the group-pair energy buffers produced by a SIMD kernel
  93  * to single terms in the output buffers.
  94  *
  95  * The SIMD kernels produce a large number of energy buffer in SIMD registers
  96  * to avoid scattered reads and writes.
  97  *
  98  * \tparam        unrollj         The unroll size for j-particles in the SIMD kernel
  99  * \param[in]     numGroups       The number of energy groups
 100  * \param[in]     numGroups_2log  Log2 of numGroups, rounded up
 101  * \param[in,out] out             Struct with energy buffers
 102  */
 103 template<int unrollj>
 104 static void reduceGroupEnergySimdBuffers(int numGroups, int numGroups_2log, nbnxn_atomdata_output_t* out)
 105 {
 106     const int unrollj_half = unrollj / 2;
 107     /* Energies are stored in SIMD registers with size 2^numGroups_2log */
 108     const int numGroupsStorage = (1 << numGroups_2log);
 109
 110     const real* gmx_restrict vVdwSimd     = out->VSvdw.data();
 111     const real* gmx_restrict vCoulombSimd = out->VSc.data();
 112     real* gmx_restrict vVdw               = out->Vvdw.data();
 113     real* gmx_restrict vCoulomb           = out->Vc.data();
 114
 115     /* The size of the SIMD energy group buffer array is:
 116      * numGroups*numGroups*numGroupsStorage*unrollj_half*simd_width
 117      */
 118     for (int i = 0; i < numGroups; i++)
 119     {
 120         for (int j1 = 0; j1 < numGroups; j1++)
 121         {
 122             for (int j0 = 0; j0 < numGroups; j0++)
 123             {
 124                 int c = ((i * numGroups + j1) * numGroupsStorage + j0) * unrollj_half * unrollj;
 125                 for (int s = 0; s < unrollj_half; s++)
 126                 {
 127                     vVdw[i * numGroups + j0] += vVdwSimd[c + 0];
 128                     vVdw[i * numGroups + j1] += vVdwSimd[c + 1];
 129                     vCoulomb[i * numGroups + j0] += vCoulombSimd[c + 0];
 130                     vCoulomb[i * numGroups + j1] += vCoulombSimd[c + 1];
 131                     c += unrollj + 2;
 132                 }
 133             }
 134         }
 135     }
 136 }
 137
 138 static int getCoulombKernelType(const Nbnxm::KernelSetup& kernelSetup, const interaction_const_t& ic)
 139 {
 140
 141     if (EEL_RF(ic.eeltype) || ic.eeltype == CoulombInteractionType::Cut)
 142     {
 143         return coulktRF;
 144     }
 145     else
 146     {
 147         if (kernelSetup.ewaldExclusionType == Nbnxm::EwaldExclusionType::Table)
 148         {
 149             if (ic.rcoulomb == ic.rvdw)
 150             {
 151                 return coulktTAB;
 152             }
 153             else
 154             {
 155                 return coulktTAB_TWIN;
 156             }
 157         }
 158         else
 159         {
 160             if (ic.rcoulomb == ic.rvdw)
 161             {
 162                 return coulktEWALD;
 163             }
 164             else
 165             {
 166                 return coulktEWALD_TWIN;
 167             }
 168         }
 169     }
 170 }
 171
 172 static int getVdwKernelType(const Nbnxm::KernelSetup&       kernelSetup,
 173                             const nbnxn_atomdata_t::Params& nbatParams,
 174                             const interaction_const_t&      ic)
 175 {
 176     if (ic.vdwtype == VanDerWaalsType::Cut)
 177     {
 178         switch (ic.vdw_modifier)
 179         {
 180             case InteractionModifiers::None:
 181             case InteractionModifiers::PotShift:
 182                 switch (nbatParams.ljCombinationRule)
 183                 {
 184                     case LJCombinationRule::Geometric: return vdwktLJCUT_COMBGEOM;
 185                     case LJCombinationRule::LorentzBerthelot: return vdwktLJCUT_COMBLB;
 186                     case LJCombinationRule::None: return vdwktLJCUT_COMBNONE;
 187                     default: gmx_incons("Unknown combination rule");
 188                 }
 189             case InteractionModifiers::ForceSwitch: return vdwktLJFORCESWITCH;
 190             case InteractionModifiers::PotSwitch: return vdwktLJPOTSWITCH;
 191             default:
 192                 std::string errorMsg =
 193                         gmx::formatString("Unsupported VdW interaction modifier %s (%d)",
 194                                           enumValueToString(ic.vdw_modifier),
 195                                           static_cast<int>(ic.vdw_modifier));
 196                 gmx_incons(errorMsg);
 197         }
 198     }
 199     else if (ic.vdwtype == VanDerWaalsType::Pme)
 200     {
 201         if (ic.ljpme_comb_rule == LongRangeVdW::Geom)
 202         {
 203             return vdwktLJEWALDCOMBGEOM;
 204         }
 205         else
 206         {
 207             /* At setup we (should have) selected the C reference kernel */
 208             GMX_RELEASE_ASSERT(kernelSetup.kernelType == Nbnxm::KernelType::Cpu4x4_PlainC,
 209                                "Only the C reference nbnxn SIMD kernel supports LJ-PME with LB "
 210                                "combination rules");
 211             return vdwktLJEWALDCOMBLB;
 212         }
 213     }
 214     else
 215     {
 216         std::string errorMsg = gmx::formatString("Unsupported VdW interaction type %s (%d)",
 217                                                  enumValueToString(ic.vdwtype),
 218                                                  static_cast<int>(ic.vdwtype));
 219         gmx_incons(errorMsg);
 220     }
 221 }
 222
 223 /*! \brief Dispatches the non-bonded N versus M atom cluster CPU kernels.
 224  *
 225  * OpenMP parallelization is performed within this function.
 226  * Energy reduction, but not force and shift force reduction, is performed
 227  * within this function.
 228  *
 229  * \param[in]     pairlistSet   Pairlists with local or non-local interactions to compute
 230  * \param[in]     kernelSetup   The non-bonded kernel setup
 231  * \param[in,out] nbat          The atomdata for the interactions
 232  * \param[in]     ic            Non-bonded interaction constants
 233  * \param[in]     shiftVectors  The PBC shift vectors
 234  * \param[in]     stepWork      Flags that tell what to compute
 235  * \param[in]     clearF        Enum that tells if to clear the force output buffer
 236  * \param[out]    vCoulomb      Output buffer for Coulomb energies
 237  * \param[out]    vVdw          Output buffer for Van der Waals energies
 238  * \param[in]     wcycle        Pointer to cycle counting data structure.
 239  */
 240 static void nbnxn_kernel_cpu(const PairlistSet&             pairlistSet,
 241                              const Nbnxm::KernelSetup&      kernelSetup,
 242                              nbnxn_atomdata_t*              nbat,
 243                              const interaction_const_t&     ic,
 244                              gmx::ArrayRef<const gmx::RVec> shiftVectors,
 245                              const gmx::StepWorkload&       stepWork,
 246                              int                            clearF,
 247                              real*                          vCoulomb,
 248                              real*                          vVdw,
 249                              gmx_wallcycle*                 wcycle)
 250 {
 251
 252     const nbnxn_atomdata_t::Params& nbatParams = nbat->params();
 253
 254     const int coulkt = getCoulombKernelType(kernelSetup, ic);
 255     const int vdwkt  = getVdwKernelType(kernelSetup, nbatParams, ic);
 256
 257     gmx::ArrayRef<const NbnxnPairlistCpu> pairlists = pairlistSet.cpuLists();
 258
 259     const auto* shiftVecPointer = as_rvec_array(shiftVectors.data());
 260
 261     int gmx_unused nthreads = gmx_omp_nthreads_get(ModuleMultiThread::Nonbonded);
 262     wallcycle_sub_start(wcycle, WallCycleSubCounter::NonbondedClear);
 263 #pragma omp parallel for schedule(static) num_threads(nthreads)
 264     for (gmx::index nb = 0; nb < pairlists.ssize(); nb++)
 265     {
 266         // Presently, the kernels do not call C++ code that can throw,
 267         // so no need for a try/catch pair in this OpenMP region.
 268         nbnxn_atomdata_output_t* out = &nbat->out[nb];
 269
 270         if (clearF == enbvClearFYes)
 271         {
 272             clearForceBuffer(nbat, nb);
 273
 274             clear_fshift(out->fshift.data());
 275         }
 276
 277         if (nb == 0)
 278         {
 279             wallcycle_sub_stop(wcycle, WallCycleSubCounter::NonbondedClear);
 280             wallcycle_sub_start(wcycle, WallCycleSubCounter::NonbondedKernel);
 281         }
 282
 283         // TODO: Change to reference
 284         const NbnxnPairlistCpu* pairlist = &pairlists[nb];
 285
 286         if (!stepWork.computeEnergy)
 287         {
 288             /* Don't calculate energies */
 289             switch (kernelSetup.kernelType)
 290             {
 291                 case Nbnxm::KernelType::Cpu4x4_PlainC:
 292                     nbnxn_kernel_noener_ref[coulkt][vdwkt](pairlist, nbat, &ic, shiftVecPointer, out);
 293                     break;
 294 #ifdef GMX_NBNXN_SIMD_2XNN
 295                 case Nbnxm::KernelType::Cpu4xN_Simd_2xNN:
 296                     nbnxm_kernel_noener_simd_2xmm[coulkt][vdwkt](pairlist, nbat, &ic, shiftVecPointer, out);
 297                     break;
 298 #endif
 299 #ifdef GMX_NBNXN_SIMD_4XN
 300                 case Nbnxm::KernelType::Cpu4xN_Simd_4xN:
 301                     nbnxm_kernel_noener_simd_4xm[coulkt][vdwkt](pairlist, nbat, &ic, shiftVecPointer, out);
 302                     break;
 303 #endif
 304                 default: GMX_RELEASE_ASSERT(false, "Unsupported kernel architecture");
 305             }
 306         }
 307         else if (out->Vvdw.size() == 1)
 308         {
 309             /* A single energy group (pair) */
 310             out->Vvdw[0] = 0;
 311             out->Vc[0]   = 0;
 312
 313             switch (kernelSetup.kernelType)
 314             {
 315                 case Nbnxm::KernelType::Cpu4x4_PlainC:
 316                     nbnxn_kernel_ener_ref[coulkt][vdwkt](pairlist, nbat, &ic, shiftVecPointer, out);
 317                     break;
 318 #ifdef GMX_NBNXN_SIMD_2XNN
 319                 case Nbnxm::KernelType::Cpu4xN_Simd_2xNN:
 320                     nbnxm_kernel_ener_simd_2xmm[coulkt][vdwkt](pairlist, nbat, &ic, shiftVecPointer, out);
 321                     break;
 322 #endif
 323 #ifdef GMX_NBNXN_SIMD_4XN
 324                 case Nbnxm::KernelType::Cpu4xN_Simd_4xN:
 325                     nbnxm_kernel_ener_simd_4xm[coulkt][vdwkt](pairlist, nbat, &ic, shiftVecPointer, out);
 326                     break;
 327 #endif
 328                 default: GMX_RELEASE_ASSERT(false, "Unsupported kernel architecture");
 329             }
 330         }
 331         else
 332         {
 333             /* Calculate energy group contributions */
 334             clearGroupEnergies(out);
 335
 336             int unrollj = 0;
 337
 338             switch (kernelSetup.kernelType)
 339             {
 340                 case Nbnxm::KernelType::Cpu4x4_PlainC:
 341                     unrollj = c_nbnxnCpuIClusterSize;
 342                     nbnxn_kernel_energrp_ref[coulkt][vdwkt](pairlist, nbat, &ic, shiftVecPointer, out);
 343                     break;
 344 #ifdef GMX_NBNXN_SIMD_2XNN
 345                 case Nbnxm::KernelType::Cpu4xN_Simd_2xNN:
 346                     unrollj = GMX_SIMD_REAL_WIDTH / 2;
 347                     nbnxm_kernel_energrp_simd_2xmm[coulkt][vdwkt](
 348                             pairlist, nbat, &ic, shiftVecPointer, out);
 349                     break;
 350 #endif
 351 #ifdef GMX_NBNXN_SIMD_4XN
 352                 case Nbnxm::KernelType::Cpu4xN_Simd_4xN:
 353                     unrollj = GMX_SIMD_REAL_WIDTH;
 354                     nbnxm_kernel_energrp_simd_4xm[coulkt][vdwkt](pairlist, nbat, &ic, shiftVecPointer, out);
 355                     break;
 356 #endif
 357                 default: GMX_RELEASE_ASSERT(false, "Unsupported kernel architecture");
 358             }
 359
 360             if (kernelSetup.kernelType != Nbnxm::KernelType::Cpu4x4_PlainC)
 361             {
 362                 switch (unrollj)
 363                 {
 364                     case 2:
 365                         reduceGroupEnergySimdBuffers<2>(nbatParams.nenergrp, nbatParams.neg_2log, out);
 366                         break;
 367                     case 4:
 368                         reduceGroupEnergySimdBuffers<4>(nbatParams.nenergrp, nbatParams.neg_2log, out);
 369                         break;
 370                     case 8:
 371                         reduceGroupEnergySimdBuffers<8>(nbatParams.nenergrp, nbatParams.neg_2log, out);
 372                         break;
 373                     default: GMX_RELEASE_ASSERT(false, "Unsupported j-unroll size");
 374                 }
 375             }
 376         }
 377     }
 378     wallcycle_sub_stop(wcycle, WallCycleSubCounter::NonbondedKernel);
 379
 380     if (stepWork.computeEnergy)
 381     {
 382         reduce_energies_over_lists(nbat, pairlists.ssize(), vVdw, vCoulomb);
 383     }
 384 }
 385
 386 static void accountFlops(t_nrnb*                    nrnb,
 387                          const PairlistSet&         pairlistSet,
 388                          const nonbonded_verlet_t&  nbv,
 389                          const interaction_const_t& ic,
 390                          const gmx::StepWorkload&   stepWork)
 391 {
 392     const bool usingGpuKernels = nbv.useGpu();
 393
 394     int enr_nbnxn_kernel_ljc = eNRNB;
 395     if (EEL_RF(ic.eeltype) || ic.eeltype == CoulombInteractionType::Cut)
 396     {
 397         enr_nbnxn_kernel_ljc = eNR_NBNXN_LJ_RF;
 398     }
 399     else if ((!usingGpuKernels && nbv.kernelSetup().ewaldExclusionType == Nbnxm::EwaldExclusionType::Analytical)
 400              || (usingGpuKernels && Nbnxm::gpu_is_kernel_ewald_analytical(nbv.gpu_nbv)))
 401     {
 402         enr_nbnxn_kernel_ljc = eNR_NBNXN_LJ_EWALD;
 403     }
 404     else
 405     {
 406         enr_nbnxn_kernel_ljc = eNR_NBNXN_LJ_TAB;
 407     }
 408     int enr_nbnxn_kernel_lj = eNR_NBNXN_LJ;
 409     if (stepWork.computeEnergy)
 410     {
 411         /* In eNR_??? the nbnxn F+E kernels are always the F kernel + 1 */
 412         enr_nbnxn_kernel_ljc += 1;
 413         enr_nbnxn_kernel_lj += 1;
 414     }
 415
 416     inc_nrnb(nrnb, enr_nbnxn_kernel_ljc, pairlistSet.natpair_ljq_);
 417     inc_nrnb(nrnb, enr_nbnxn_kernel_lj, pairlistSet.natpair_lj_);
 418     /* The Coulomb-only kernels are offset -eNR_NBNXN_LJ_RF+eNR_NBNXN_RF */
 419     inc_nrnb(nrnb, enr_nbnxn_kernel_ljc - eNR_NBNXN_LJ_RF + eNR_NBNXN_RF, pairlistSet.natpair_q_);
 420
 421     if (ic.vdw_modifier == InteractionModifiers::ForceSwitch)
 422     {
 423         /* We add up the switch cost separately */
 424         inc_nrnb(nrnb,
 425                  eNR_NBNXN_ADD_LJ_FSW + (stepWork.computeEnergy ? 1 : 0),
 426                  pairlistSet.natpair_ljq_ + pairlistSet.natpair_lj_);
 427     }
 428     if (ic.vdw_modifier == InteractionModifiers::PotSwitch)
 429     {
 430         /* We add up the switch cost separately */
 431         inc_nrnb(nrnb,
 432                  eNR_NBNXN_ADD_LJ_PSW + (stepWork.computeEnergy ? 1 : 0),
 433                  pairlistSet.natpair_ljq_ + pairlistSet.natpair_lj_);
 434     }
 435     if (ic.vdwtype == VanDerWaalsType::Pme)
 436     {
 437         /* We add up the LJ Ewald cost separately */
 438         inc_nrnb(nrnb,
 439                  eNR_NBNXN_ADD_LJ_EWALD + (stepWork.computeEnergy ? 1 : 0),
 440                  pairlistSet.natpair_ljq_ + pairlistSet.natpair_lj_);
 441     }
 442 }
 443
 444 void nonbonded_verlet_t::dispatchNonbondedKernel(gmx::InteractionLocality       iLocality,
 445                                                  const interaction_const_t&     ic,
 446                                                  const gmx::StepWorkload&       stepWork,
 447                                                  int                            clearF,
 448                                                  gmx::ArrayRef<const gmx::RVec> shiftvec,
 449                                                  gmx::ArrayRef<real> repulsionDispersionSR,
 450                                                  gmx::ArrayRef<real> CoulombSR,
 451                                                  t_nrnb*             nrnb) const
 452 {
 453     const PairlistSet& pairlistSet = pairlistSets().pairlistSet(iLocality);
 454
 455     switch (kernelSetup().kernelType)
 456     {
 457         case Nbnxm::KernelType::Cpu4x4_PlainC:
 458         case Nbnxm::KernelType::Cpu4xN_Simd_4xN:
 459         case Nbnxm::KernelType::Cpu4xN_Simd_2xNN:
 460             nbnxn_kernel_cpu(pairlistSet,
 461                              kernelSetup(),
 462                              nbat.get(),
 463                              ic,
 464                              shiftvec,
 465                              stepWork,
 466                              clearF,
 467                              CoulombSR.data(),
 468                              repulsionDispersionSR.data(),
 469                              wcycle_);
 470             break;
 471
 472         case Nbnxm::KernelType::Gpu8x8x8:
 473             Nbnxm::gpu_launch_kernel(gpu_nbv, stepWork, iLocality);
 474             break;
 475
 476         case Nbnxm::KernelType::Cpu8x8x8_PlainC:
 477             nbnxn_kernel_gpu_ref(pairlistSet.gpuList(),
 478                                  nbat.get(),
 479                                  &ic,
 480                                  shiftvec,
 481                                  stepWork,
 482                                  clearF,
 483                                  nbat->out[0].f,
 484                                  nbat->out[0].fshift.data(),
 485                                  CoulombSR.data(),
 486                                  repulsionDispersionSR.data());
 487             break;
 488
 489         default: GMX_RELEASE_ASSERT(false, "Invalid nonbonded kernel type passed!");
 490     }
 491
 492     accountFlops(nrnb, pairlistSet, *this, ic, stepWork);
 493 }
 494
 495 void nonbonded_verlet_t::dispatchFreeEnergyKernel(gmx::InteractionLocality       iLocality,
 496                                                   const t_forcerec&              fr,
 497                                                   gmx::ArrayRef<const gmx::RVec> coords,
 498                                                   gmx::ForceWithShiftForces* forceWithShiftForces,
 499                                                   gmx::ArrayRef<const real>  chargeA,
 500                                                   gmx::ArrayRef<const real>  chargeB,
 501                                                   gmx::ArrayRef<const int>   typeA,
 502                                                   gmx::ArrayRef<const int>   typeB,
 503                                                   t_lambda*                  fepvals,
 504                                                   gmx::ArrayRef<const real>  lambda,
 505                                                   gmx_enerdata_t*            enerd,
 506                                                   const gmx::StepWorkload&   stepWork,
 507                                                   t_nrnb*                    nrnb)
 508 {
 509     const auto nbl_fep = pairlistSets().pairlistSet(iLocality).fepLists();
 510
 511     /* When the first list is empty, all are empty and there is nothing to do */
 512     if (!pairlistSets().params().haveFep || nbl_fep[0]->nrj == 0)
 513     {
 514         return;
 515     }
 516
 517     int donb_flags = 0;
 518     /* Add short-range interactions */
 519     donb_flags |= GMX_NONBONDED_DO_SR;
 520
 521     if (stepWork.computeForces)
 522     {
 523         donb_flags |= GMX_NONBONDED_DO_FORCE;
 524     }
 525     if (stepWork.computeVirial)
 526     {
 527         donb_flags |= GMX_NONBONDED_DO_SHIFTFORCE;
 528     }
 529     if (stepWork.computeEnergy)
 530     {
 531         donb_flags |= GMX_NONBONDED_DO_POTENTIAL;
 532     }
 533
 534     gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, real> dvdl_nb      = { 0 };
 535     int                                                             kernelFlags  = donb_flags;
 536     gmx::ArrayRef<const real>                                       kernelLambda = lambda;
 537     gmx::ArrayRef<real>                                             kernelDvdl   = dvdl_nb;
 538
 539     gmx::ArrayRef<real> energygrp_elec = enerd->grpp.energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR];
 540     gmx::ArrayRef<real> energygrp_vdw = enerd->grpp.energyGroupPairTerms[NonBondedEnergyTerms::LJSR];
 541
 542     GMX_ASSERT(gmx_omp_nthreads_get(ModuleMultiThread::Nonbonded) == nbl_fep.ssize(),
 543                "Number of lists should be same as number of NB threads");
 544
 545     wallcycle_sub_start(wcycle_, WallCycleSubCounter::NonbondedFep);
 546 #pragma omp parallel for schedule(static) num_threads(nbl_fep.ssize())
 547     for (gmx::index th = 0; th < nbl_fep.ssize(); th++)
 548     {
 549         try
 550         {
 551             gmx_nb_free_energy_kernel(*nbl_fep[th],
 552                                       coords,
 553                                       forceWithShiftForces,
 554                                       fr.use_simd_kernels,
 555                                       fr.ntype,
 556                                       fr.rlist,
 557                                       *fr.ic,
 558                                       fr.shift_vec,
 559                                       fr.nbfp,
 560                                       fr.ljpme_c6grid,
 561                                       chargeA,
 562                                       chargeB,
 563                                       typeA,
 564                                       typeB,
 565                                       kernelFlags,
 566                                       kernelLambda,
 567                                       kernelDvdl,
 568                                       energygrp_elec,
 569                                       energygrp_vdw,
 570                                       nrnb);
 571         }
 572         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 573     }
 574
 575     if (fepvals->sc_alpha != 0)
 576     {
 577         enerd->dvdl_nonlin[FreeEnergyPerturbationCouplingType::Vdw] +=
 578                 dvdl_nb[FreeEnergyPerturbationCouplingType::Vdw];
 579         enerd->dvdl_nonlin[FreeEnergyPerturbationCouplingType::Coul] +=
 580                 dvdl_nb[FreeEnergyPerturbationCouplingType::Coul];
 581     }
 582     else
 583     {
 584         enerd->dvdl_lin[FreeEnergyPerturbationCouplingType::Vdw] +=
 585                 dvdl_nb[FreeEnergyPerturbationCouplingType::Vdw];
 586         enerd->dvdl_lin[FreeEnergyPerturbationCouplingType::Coul] +=
 587                 dvdl_nb[FreeEnergyPerturbationCouplingType::Coul];
 588     }
 589
 590     /* If we do foreign lambda and we have soft-core interactions
 591      * we have to recalculate the (non-linear) energies contributions.
 592      */
 593     if (fepvals->n_lambda > 0 && stepWork.computeDhdl && fepvals->sc_alpha != 0)
 594     {
 595         gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, real> lam_i;
 596         kernelFlags = (donb_flags & ~(GMX_NONBONDED_DO_FORCE | GMX_NONBONDED_DO_SHIFTFORCE))
 597                       | GMX_NONBONDED_DO_FOREIGNLAMBDA;
 598         kernelLambda = lam_i;
 599         kernelDvdl   = dvdl_nb;
 600         gmx::ArrayRef<real> energygrp_elec =
 601                 enerd->foreign_grpp.energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR];
 602         gmx::ArrayRef<real> energygrp_vdw =
 603                 enerd->foreign_grpp.energyGroupPairTerms[NonBondedEnergyTerms::LJSR];
 604
 605         for (gmx::index i = 0; i < 1 + enerd->foreignLambdaTerms.numLambdas(); i++)
 606         {
 607             std::fill(std::begin(dvdl_nb), std::end(dvdl_nb), 0);
 608             for (int j = 0; j < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); j++)
 609             {
 610                 lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i - 1]);
 611             }
 612             reset_foreign_enerdata(enerd);
 613 #pragma omp parallel for schedule(static) num_threads(nbl_fep.ssize())
 614             for (gmx::index th = 0; th < nbl_fep.ssize(); th++)
 615             {
 616                 try
 617                 {
 618                     gmx_nb_free_energy_kernel(*nbl_fep[th],
 619                                               coords,
 620                                               forceWithShiftForces,
 621                                               fr.use_simd_kernels,
 622                                               fr.ntype,
 623                                               fr.rlist,
 624                                               *fr.ic,
 625                                               fr.shift_vec,
 626                                               fr.nbfp,
 627                                               fr.ljpme_c6grid,
 628                                               chargeA,
 629                                               chargeB,
 630                                               typeA,
 631                                               typeB,
 632                                               kernelFlags,
 633                                               kernelLambda,
 634                                               kernelDvdl,
 635                                               energygrp_elec,
 636                                               energygrp_vdw,
 637                                               nrnb);
 638                 }
 639                 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 640             }
 641
 642             sum_epot(enerd->foreign_grpp, enerd->foreign_term.data());
 643             enerd->foreignLambdaTerms.accumulate(
 644                     i,
 645                     enerd->foreign_term[F_EPOT],
 646                     dvdl_nb[FreeEnergyPerturbationCouplingType::Vdw]
 647                             + dvdl_nb[FreeEnergyPerturbationCouplingType::Coul]);
 648         }
 649     }
 650     wallcycle_sub_stop(wcycle_, WallCycleSubCounter::NonbondedFep);
 651 }