src/gromacs/listed_forces/listed_forces.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2014,2015,2016,2017,2018 by the GROMACS development team.
   5  * Copyright (c) 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 /*! \internal \file
  37  *
  38  * \brief This file defines high-level functions for mdrun to compute
  39  * energies and forces for listed interactions.
  40  *
  41  * \author Mark Abraham <mark.j.abraham@gmail.com>
  42  *
  43  * \ingroup module_listed_forces
  44  */
  45 #include "gmxpre.h"
  46
  47 #include "gromacs/utility/arrayref.h"
  48 #include "gromacs/utility/enumerationhelpers.h"
  49 #include "listed_forces.h"
  50
  51 #include <cassert>
  52
  53 #include <algorithm>
  54 #include <array>
  55 #include <numeric>
  56
  57 #include "gromacs/gmxlib/network.h"
  58 #include "gromacs/gmxlib/nrnb.h"
  59 #include "gromacs/listed_forces/bonded.h"
  60 #include "gromacs/listed_forces/disre.h"
  61 #include "gromacs/listed_forces/orires.h"
  62 #include "gromacs/listed_forces/pairs.h"
  63 #include "gromacs/listed_forces/position_restraints.h"
  64 #include "gromacs/math/vec.h"
  65 #include "gromacs/mdlib/enerdata_utils.h"
  66 #include "gromacs/mdlib/force.h"
  67 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  68 #include "gromacs/mdtypes/commrec.h"
  69 #include "gromacs/mdtypes/mdatom.h"
  70 #include "gromacs/mdtypes/fcdata.h"
  71 #include "gromacs/mdtypes/forceoutput.h"
  72 #include "gromacs/mdtypes/forcerec.h"
  73 #include "gromacs/mdtypes/inputrec.h"
  74 #include "gromacs/mdtypes/md_enums.h"
  75 #include "gromacs/mdtypes/simulation_workload.h"
  76 #include "gromacs/pbcutil/ishift.h"
  77 #include "gromacs/pbcutil/pbc.h"
  78 #include "gromacs/timing/wallcycle.h"
  79 #include "gromacs/topology/topology.h"
  80 #include "gromacs/utility/exceptions.h"
  81 #include "gromacs/utility/fatalerror.h"
  82
  83 #include "listed_internal.h"
  84 #include "manage_threading.h"
  85 #include "utilities.h"
  86
  87 ListedForces::ListedForces(const gmx_ffparams_t&      ffparams,
  88                            const int                  numEnergyGroups,
  89                            const int                  numThreads,
  90                            const InteractionSelection interactionSelection,
  91                            FILE*                      fplog) :
  92     idefSelection_(ffparams),
  93     threading_(std::make_unique<bonded_threading_t>(numThreads, numEnergyGroups, fplog)),
  94     interactionSelection_(interactionSelection),
  95     foreignEnergyGroups_(std::make_unique<gmx_grppairener_t>(numEnergyGroups))
  96 {
  97 }
  98
  99 ListedForces::ListedForces(ListedForces&& o) noexcept = default;
 100
 101 ListedForces::~ListedForces() = default;
 102
 103 //! Copies the selection interactions from \p idefSrc to \p idef
 104 static void selectInteractions(InteractionDefinitions*                  idef,
 105                                const InteractionDefinitions&            idefSrc,
 106                                const ListedForces::InteractionSelection interactionSelection)
 107 {
 108     const bool selectPairs =
 109             interactionSelection.test(static_cast<int>(ListedForces::InteractionGroup::Pairs));
 110     const bool selectDihedrals =
 111             interactionSelection.test(static_cast<int>(ListedForces::InteractionGroup::Dihedrals));
 112     const bool selectAngles =
 113             interactionSelection.test(static_cast<int>(ListedForces::InteractionGroup::Angles));
 114     const bool selectRest =
 115             interactionSelection.test(static_cast<int>(ListedForces::InteractionGroup::Rest));
 116
 117     for (int ftype = 0; ftype < F_NRE; ftype++)
 118     {
 119         const t_interaction_function& ifunc = interaction_function[ftype];
 120         if (ifunc.flags & IF_BOND)
 121         {
 122             bool assign = false;
 123             if (ifunc.flags & IF_PAIR)
 124             {
 125                 assign = selectPairs;
 126             }
 127             else if (ifunc.flags & IF_DIHEDRAL)
 128             {
 129                 assign = selectDihedrals;
 130             }
 131             else if (ifunc.flags & IF_ATYPE)
 132             {
 133                 assign = selectAngles;
 134             }
 135             else
 136             {
 137                 assign = selectRest;
 138             }
 139             if (assign)
 140             {
 141                 idef->il[ftype] = idefSrc.il[ftype];
 142             }
 143             else
 144             {
 145                 idef->il[ftype].clear();
 146             }
 147         }
 148     }
 149 }
 150
 151 void ListedForces::setup(const InteractionDefinitions& domainIdef, const int numAtomsForce, const bool useGpu)
 152 {
 153     if (interactionSelection_.all())
 154     {
 155         // Avoid the overhead of copying all interaction lists by simply setting the reference to the domain idef
 156         idef_ = &domainIdef;
 157     }
 158     else
 159     {
 160         idef_ = &idefSelection_;
 161
 162         selectInteractions(&idefSelection_, domainIdef, interactionSelection_);
 163
 164         idefSelection_.ilsort = domainIdef.ilsort;
 165
 166         if (interactionSelection_.test(static_cast<int>(ListedForces::InteractionGroup::Rest)))
 167         {
 168             idefSelection_.iparams_posres   = domainIdef.iparams_posres;
 169             idefSelection_.iparams_fbposres = domainIdef.iparams_fbposres;
 170         }
 171         else
 172         {
 173             idefSelection_.iparams_posres.clear();
 174             idefSelection_.iparams_fbposres.clear();
 175         }
 176     }
 177
 178     setup_bonded_threading(threading_.get(), numAtomsForce, useGpu, *idef_);
 179
 180     if (idef_->ilsort == ilsortFE_SORTED)
 181     {
 182         forceBufferLambda_.resize(numAtomsForce * sizeof(rvec4) / sizeof(real));
 183         shiftForceBufferLambda_.resize(gmx::c_numShiftVectors);
 184     }
 185 }
 186
 187 namespace
 188 {
 189
 190 using gmx::ArrayRef;
 191
 192 /*! \brief Return true if ftype is an explicit pair-listed LJ or
 193  * COULOMB interaction type: bonded LJ (usually 1-4), or special
 194  * listed non-bonded for FEP. */
 195 bool isPairInteraction(int ftype)
 196 {
 197     return ((ftype) >= F_LJ14 && (ftype) <= F_LJC_PAIRS_NB);
 198 }
 199
 200 /*! \brief Returns the bonded kernel flavor
 201  *
 202  * Note that energies are always requested when the virial
 203  * is requested (performance gain would be small).
 204  * Note that currently we do not have bonded kernels that
 205  * do not compute forces.
 206  */
 207 BondedKernelFlavor selectBondedKernelFlavor(const gmx::StepWorkload& stepWork,
 208                                             const bool               useSimdKernels,
 209                                             const bool               havePerturbedInteractions)
 210 {
 211     BondedKernelFlavor flavor;
 212     if (stepWork.computeEnergy || stepWork.computeVirial)
 213     {
 214         if (stepWork.computeVirial)
 215         {
 216             flavor = BondedKernelFlavor::ForcesAndVirialAndEnergy;
 217         }
 218         else
 219         {
 220             flavor = BondedKernelFlavor::ForcesAndEnergy;
 221         }
 222     }
 223     else
 224     {
 225         if (useSimdKernels && !havePerturbedInteractions)
 226         {
 227             flavor = BondedKernelFlavor::ForcesSimdWhenAvailable;
 228         }
 229         else
 230         {
 231             flavor = BondedKernelFlavor::ForcesNoSimd;
 232         }
 233     }
 234
 235     return flavor;
 236 }
 237
 238 /*! \brief Calculate one element of the list of bonded interactions
 239     for this thread */
 240 real calc_one_bond(int                           thread,
 241                    int                           ftype,
 242                    const InteractionDefinitions& idef,
 243                    ArrayRef<const int>           iatoms,
 244                    const int                     numNonperturbedInteractions,
 245                    const WorkDivision&           workDivision,
 246                    const rvec                    x[],
 247                    rvec4                         f[],
 248                    rvec                          fshift[],
 249                    const t_forcerec*             fr,
 250                    const t_pbc*                  pbc,
 251                    gmx_grppairener_t*            grpp,
 252                    t_nrnb*                       nrnb,
 253                    gmx::ArrayRef<const real>     lambda,
 254                    gmx::ArrayRef<real>           dvdl,
 255                    const t_mdatoms*              md,
 256                    t_fcdata*                     fcd,
 257                    const gmx::StepWorkload&      stepWork,
 258                    int*                          global_atom_index)
 259 {
 260     GMX_ASSERT(idef.ilsort == ilsortNO_FE || idef.ilsort == ilsortFE_SORTED,
 261                "The topology should be marked either as no FE or sorted on FE");
 262
 263     const bool havePerturbedInteractions =
 264             (idef.ilsort == ilsortFE_SORTED && numNonperturbedInteractions < iatoms.ssize());
 265     BondedKernelFlavor flavor =
 266             selectBondedKernelFlavor(stepWork, fr->use_simd_kernels, havePerturbedInteractions);
 267     FreeEnergyPerturbationCouplingType efptFTYPE;
 268     if (IS_RESTRAINT_TYPE(ftype))
 269     {
 270         efptFTYPE = FreeEnergyPerturbationCouplingType::Restraint;
 271     }
 272     else
 273     {
 274         efptFTYPE = FreeEnergyPerturbationCouplingType::Bonded;
 275     }
 276
 277     const int nat1   = interaction_function[ftype].nratoms + 1;
 278     const int nbonds = iatoms.ssize() / nat1;
 279
 280     GMX_ASSERT(fr->listedForcesGpu != nullptr || workDivision.end(ftype) == iatoms.ssize(),
 281                "The thread division should match the topology");
 282
 283     const int nb0 = workDivision.bound(ftype, thread);
 284     const int nbn = workDivision.bound(ftype, thread + 1) - nb0;
 285
 286     ArrayRef<const t_iparams> iparams = idef.iparams;
 287
 288     real v = 0;
 289     if (!isPairInteraction(ftype))
 290     {
 291         if (ftype == F_CMAP)
 292         {
 293             /* TODO The execution time for CMAP dihedrals might be
 294                nice to account to its own subtimer, but first
 295                wallcycle needs to be extended to support calling from
 296                multiple threads. */
 297             v = cmap_dihs(nbn,
 298                           iatoms.data() + nb0,
 299                           iparams.data(),
 300                           &idef.cmap_grid,
 301                           x,
 302                           f,
 303                           fshift,
 304                           pbc,
 305                           lambda[static_cast<int>(efptFTYPE)],
 306                           &(dvdl[static_cast<int>(efptFTYPE)]),
 307                           gmx::arrayRefFromArray(md->chargeA, md->nr),
 308                           fcd,
 309                           nullptr,
 310                           nullptr,
 311                           global_atom_index);
 312         }
 313         else
 314         {
 315             v = calculateSimpleBond(ftype,
 316                                     nbn,
 317                                     iatoms.data() + nb0,
 318                                     iparams.data(),
 319                                     x,
 320                                     f,
 321                                     fshift,
 322                                     pbc,
 323                                     lambda[static_cast<int>(efptFTYPE)],
 324                                     &(dvdl[static_cast<int>(efptFTYPE)]),
 325                                     gmx::arrayRefFromArray(md->chargeA, md->nr),
 326                                     fcd,
 327                                     fcd->disres,
 328                                     fcd->orires.get(),
 329                                     global_atom_index,
 330                                     flavor);
 331         }
 332     }
 333     else
 334     {
 335         /* TODO The execution time for pairs might be nice to account
 336            to its own subtimer, but first wallcycle needs to be
 337            extended to support calling from multiple threads. */
 338         do_pairs(ftype,
 339                  nbn,
 340                  iatoms.data() + nb0,
 341                  iparams.data(),
 342                  x,
 343                  f,
 344                  fshift,
 345                  pbc,
 346                  lambda.data(),
 347                  dvdl.data(),
 348                  md->chargeA ? gmx::arrayRefFromArray(md->chargeA, md->nr) : gmx::ArrayRef<real>{},
 349                  md->chargeB ? gmx::arrayRefFromArray(md->chargeB, md->nr) : gmx::ArrayRef<real>{},
 350                  md->bPerturbed ? gmx::arrayRefFromArray(md->bPerturbed, md->nr) : gmx::ArrayRef<bool>(),
 351                  gmx::arrayRefFromArray(md->cENER, md->nr),
 352                  md->nPerturbed,
 353                  fr,
 354                  havePerturbedInteractions,
 355                  stepWork,
 356                  grpp,
 357                  global_atom_index);
 358     }
 359
 360     if (thread == 0)
 361     {
 362         inc_nrnb(nrnb, nrnbIndex(ftype), nbonds);
 363     }
 364
 365     return v;
 366 }
 367
 368 } // namespace
 369
 370 /*! \brief Compute the bonded part of the listed forces, parallelized over threads
 371  */
 372 static void calcBondedForces(const InteractionDefinitions& idef,
 373                              bonded_threading_t*           bt,
 374                              const rvec                    x[],
 375                              const t_forcerec*             fr,
 376                              const t_pbc*                  pbc_null,
 377                              rvec*                         fshiftMasterBuffer,
 378                              gmx_enerdata_t*               enerd,
 379                              t_nrnb*                       nrnb,
 380                              gmx::ArrayRef<const real>     lambda,
 381                              gmx::ArrayRef<real>           dvdl,
 382                              const t_mdatoms*              md,
 383                              t_fcdata*                     fcd,
 384                              const gmx::StepWorkload&      stepWork,
 385                              int*                          global_atom_index)
 386 {
 387 #pragma omp parallel for num_threads(bt->nthreads) schedule(static)
 388     for (int thread = 0; thread < bt->nthreads; thread++)
 389     {
 390         try
 391         {
 392             auto& threadBuffer = bt->threadedForceBuffer.threadForceBuffer(thread);
 393             /* thread stuff */
 394             rvec*               fshift;
 395             gmx::ArrayRef<real> dvdlt;
 396             gmx::ArrayRef<real> epot;
 397             gmx_grppairener_t*  grpp;
 398
 399             threadBuffer.clearForcesAndEnergies();
 400
 401             rvec4* ft = threadBuffer.forceBuffer();
 402
 403             /* Thread 0 writes directly to the main output buffers.
 404              * We might want to reconsider this.
 405              */
 406             if (thread == 0)
 407             {
 408                 fshift = fshiftMasterBuffer;
 409                 epot   = enerd->term;
 410                 grpp   = &enerd->grpp;
 411                 dvdlt  = dvdl;
 412             }
 413             else
 414             {
 415                 fshift = as_rvec_array(threadBuffer.shiftForces().data());
 416                 epot   = threadBuffer.energyTerms();
 417                 grpp   = &threadBuffer.groupPairEnergies();
 418                 dvdlt  = threadBuffer.dvdl();
 419             }
 420             /* Loop over all bonded force types to calculate the bonded forces */
 421             for (int ftype = 0; (ftype < F_NRE); ftype++)
 422             {
 423                 const InteractionList& ilist = idef.il[ftype];
 424                 if (!ilist.empty() && ftype_is_bonded_potential(ftype))
 425                 {
 426                     ArrayRef<const int> iatoms = gmx::makeConstArrayRef(ilist.iatoms);
 427
 428                     real v = calc_one_bond(thread,
 429                                            ftype,
 430                                            idef,
 431                                            iatoms,
 432                                            idef.numNonperturbedInteractions[ftype],
 433                                            bt->workDivision,
 434                                            x,
 435                                            ft,
 436                                            fshift,
 437                                            fr,
 438                                            pbc_null,
 439                                            grpp,
 440                                            nrnb,
 441                                            lambda,
 442                                            dvdlt,
 443                                            md,
 444                                            fcd,
 445                                            stepWork,
 446                                            global_atom_index);
 447                     epot[ftype] += v;
 448                 }
 449             }
 450         }
 451         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 452     }
 453 }
 454
 455 bool ListedForces::haveRestraints(const t_fcdata& fcdata) const
 456 {
 457     GMX_ASSERT(fcdata.disres, "NMR distance restraint object should be set up");
 458
 459     return (!idef_->il[F_POSRES].empty() || !idef_->il[F_FBPOSRES].empty() || fcdata.orires
 460             || fcdata.disres->nres > 0);
 461 }
 462
 463 bool ListedForces::haveCpuBondeds() const
 464 {
 465     return threading_->haveBondeds;
 466 }
 467
 468 bool ListedForces::haveCpuListedForces(const t_fcdata& fcdata) const
 469 {
 470     return haveCpuBondeds() || haveRestraints(fcdata);
 471 }
 472
 473 namespace
 474 {
 475
 476 /*! \brief Calculates all listed force interactions. */
 477 void calc_listed(struct gmx_wallcycle*         wcycle,
 478                  const InteractionDefinitions& idef,
 479                  bonded_threading_t*           bt,
 480                  const rvec                    x[],
 481                  gmx::ForceOutputs*            forceOutputs,
 482                  const t_forcerec*             fr,
 483                  const t_pbc*                  pbc,
 484                  gmx_enerdata_t*               enerd,
 485                  t_nrnb*                       nrnb,
 486                  gmx::ArrayRef<const real>     lambda,
 487                  const t_mdatoms*              md,
 488                  t_fcdata*                     fcd,
 489                  int*                          global_atom_index,
 490                  const gmx::StepWorkload&      stepWork)
 491 {
 492     if (bt->haveBondeds)
 493     {
 494         gmx::ForceWithShiftForces& forceWithShiftForces = forceOutputs->forceWithShiftForces();
 495
 496         wallcycle_sub_start(wcycle, WallCycleSubCounter::Listed);
 497         /* The dummy array is to have a place to store the dhdl at other values
 498            of lambda, which will be thrown away in the end */
 499         gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, real> dvdl = { 0 };
 500         calcBondedForces(idef,
 501                          bt,
 502                          x,
 503                          fr,
 504                          fr->bMolPBC ? pbc : nullptr,
 505                          as_rvec_array(forceWithShiftForces.shiftForces().data()),
 506                          enerd,
 507                          nrnb,
 508                          lambda,
 509                          dvdl,
 510                          md,
 511                          fcd,
 512                          stepWork,
 513                          global_atom_index);
 514         wallcycle_sub_stop(wcycle, WallCycleSubCounter::Listed);
 515
 516         wallcycle_sub_start(wcycle, WallCycleSubCounter::ListedBufOps);
 517         bt->threadedForceBuffer.reduce(
 518                 &forceWithShiftForces, enerd->term.data(), &enerd->grpp, dvdl, stepWork, 1);
 519
 520         if (stepWork.computeDhdl)
 521         {
 522             for (auto i : keysOf(enerd->dvdl_lin))
 523             {
 524                 enerd->dvdl_nonlin[i] += dvdl[i];
 525             }
 526         }
 527         wallcycle_sub_stop(wcycle, WallCycleSubCounter::ListedBufOps);
 528     }
 529
 530     /* Copy the sum of violations for the distance restraints from fcd */
 531     if (fcd)
 532     {
 533         enerd->term[F_DISRESVIOL] = fcd->disres->sumviol;
 534     }
 535 }
 536
 537 /*! \brief As calc_listed(), but only determines the potential energy
 538  * for the perturbed interactions.
 539  *
 540  * The shift forces in fr are not affected.
 541  */
 542 void calc_listed_lambda(const InteractionDefinitions& idef,
 543                         bonded_threading_t*           bt,
 544                         const rvec                    x[],
 545                         const t_forcerec*             fr,
 546                         const struct t_pbc*           pbc,
 547                         gmx::ArrayRef<real>           forceBufferLambda,
 548                         gmx::ArrayRef<gmx::RVec>      shiftForceBufferLambda,
 549                         gmx_grppairener_t*            grpp,
 550                         gmx::ArrayRef<real>           epot,
 551                         gmx::ArrayRef<real>           dvdl,
 552                         t_nrnb*                       nrnb,
 553                         gmx::ArrayRef<const real>     lambda,
 554                         const t_mdatoms*              md,
 555                         t_fcdata*                     fcd,
 556                         int*                          global_atom_index)
 557 {
 558     WorkDivision& workDivision = bt->foreignLambdaWorkDivision;
 559
 560     const t_pbc* pbc_null;
 561     if (fr->bMolPBC)
 562     {
 563         pbc_null = pbc;
 564     }
 565     else
 566     {
 567         pbc_null = nullptr;
 568     }
 569
 570     /* We already have the forces, so we use temp buffers here */
 571     std::fill(forceBufferLambda.begin(), forceBufferLambda.end(), 0.0_real);
 572     std::fill(shiftForceBufferLambda.begin(),
 573               shiftForceBufferLambda.end(),
 574               gmx::RVec{ 0.0_real, 0.0_real, 0.0_real });
 575     rvec4* f      = reinterpret_cast<rvec4*>(forceBufferLambda.data());
 576     rvec*  fshift = as_rvec_array(shiftForceBufferLambda.data());
 577
 578     /* Loop over all bonded force types to calculate the bonded energies */
 579     for (int ftype = 0; (ftype < F_NRE); ftype++)
 580     {
 581         if (ftype_is_bonded_potential(ftype))
 582         {
 583             const InteractionList& ilist = idef.il[ftype];
 584             /* Create a temporary iatom list with only perturbed interactions */
 585             const int           numNonperturbed = idef.numNonperturbedInteractions[ftype];
 586             ArrayRef<const int> iatomsPerturbed = gmx::constArrayRefFromArray(
 587                     ilist.iatoms.data() + numNonperturbed, ilist.size() - numNonperturbed);
 588             if (!iatomsPerturbed.empty())
 589             {
 590                 /* Set the work range of thread 0 to the perturbed bondeds */
 591                 workDivision.setBound(ftype, 0, 0);
 592                 workDivision.setBound(ftype, 1, iatomsPerturbed.ssize());
 593
 594                 gmx::StepWorkload tempFlags;
 595                 tempFlags.computeEnergy = true;
 596                 real v                  = calc_one_bond(0,
 597                                        ftype,
 598                                        idef,
 599                                        iatomsPerturbed,
 600                                        iatomsPerturbed.ssize(),
 601                                        workDivision,
 602                                        x,
 603                                        f,
 604                                        fshift,
 605                                        fr,
 606                                        pbc_null,
 607                                        grpp,
 608                                        nrnb,
 609                                        lambda,
 610                                        dvdl,
 611                                        md,
 612                                        fcd,
 613                                        tempFlags,
 614                                        global_atom_index);
 615                 epot[ftype] += v;
 616             }
 617         }
 618     }
 619 }
 620
 621 } // namespace
 622
 623 void ListedForces::calculate(struct gmx_wallcycle*                     wcycle,
 624                              const matrix                              box,
 625                              const t_lambda*                           fepvals,
 626                              const t_commrec*                          cr,
 627                              const gmx_multisim_t*                     ms,
 628                              gmx::ArrayRefWithPadding<const gmx::RVec> coordinates,
 629                              gmx::ArrayRef<const gmx::RVec>            xWholeMolecules,
 630                              t_fcdata*                                 fcdata,
 631                              const history_t*                          hist,
 632                              gmx::ForceOutputs*                        forceOutputs,
 633                              const t_forcerec*                         fr,
 634                              const struct t_pbc*                       pbc,
 635                              gmx_enerdata_t*                           enerd,
 636                              t_nrnb*                                   nrnb,
 637                              gmx::ArrayRef<const real>                 lambda,
 638                              const t_mdatoms*                          md,
 639                              int*                                      global_atom_index,
 640                              const gmx::StepWorkload&                  stepWork)
 641 {
 642     if (interactionSelection_.none() || !stepWork.computeListedForces)
 643     {
 644         return;
 645     }
 646
 647     const InteractionDefinitions& idef = *idef_;
 648
 649     // Todo: replace all rvec use here with ArrayRefWithPadding
 650     const rvec* x = as_rvec_array(coordinates.paddedArrayRef().data());
 651
 652     const bool calculateRestInteractions =
 653             interactionSelection_.test(static_cast<int>(ListedForces::InteractionGroup::Rest));
 654
 655     t_pbc pbc_full; /* Full PBC is needed for position restraints */
 656     if (calculateRestInteractions && haveRestraints(*fcdata))
 657     {
 658         if (!idef.il[F_POSRES].empty() || !idef.il[F_FBPOSRES].empty())
 659         {
 660             /* Not enough flops to bother counting */
 661             set_pbc(&pbc_full, fr->pbcType, box);
 662         }
 663
 664         /* TODO Use of restraints triggers further function calls
 665            inside the loop over calc_one_bond(), but those are too
 666            awkward to account to this subtimer properly in the present
 667            code. We don't test / care much about performance with
 668            restraints, anyway. */
 669         wallcycle_sub_start(wcycle, WallCycleSubCounter::Restraints);
 670
 671         if (!idef.il[F_POSRES].empty())
 672         {
 673             posres_wrapper(nrnb, idef, &pbc_full, x, enerd, lambda, fr, &forceOutputs->forceWithVirial());
 674         }
 675
 676         if (!idef.il[F_FBPOSRES].empty())
 677         {
 678             fbposres_wrapper(nrnb, idef, &pbc_full, x, enerd, fr, &forceOutputs->forceWithVirial());
 679         }
 680
 681         /* Do pre force calculation stuff which might require communication */
 682         if (fcdata->orires)
 683         {
 684             GMX_ASSERT(!xWholeMolecules.empty(), "Need whole molecules for orienation restraints");
 685             enerd->term[F_ORIRESDEV] = calc_orires_dev(ms,
 686                                                        idef.il[F_ORIRES].size(),
 687                                                        idef.il[F_ORIRES].iatoms.data(),
 688                                                        idef.iparams.data(),
 689                                                        xWholeMolecules,
 690                                                        x,
 691                                                        fr->bMolPBC ? pbc : nullptr,
 692                                                        fcdata->orires.get());
 693         }
 694         if (fcdata->disres->nres > 0)
 695         {
 696             calc_disres_R_6(cr,
 697                             ms,
 698                             idef.il[F_DISRES].size(),
 699                             idef.il[F_DISRES].iatoms.data(),
 700                             x,
 701                             fr->bMolPBC ? pbc : nullptr,
 702                             fcdata->disres,
 703                             hist);
 704         }
 705
 706         wallcycle_sub_stop(wcycle, WallCycleSubCounter::Restraints);
 707     }
 708
 709     calc_listed(wcycle, idef, threading_.get(), x, forceOutputs, fr, pbc, enerd, nrnb, lambda, md, fcdata, global_atom_index, stepWork);
 710
 711     /* Check if we have to determine energy differences
 712      * at foreign lambda's.
 713      */
 714     if (fepvals->n_lambda > 0 && stepWork.computeDhdl)
 715     {
 716         gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, real> dvdl = { 0 };
 717         if (!idef.il[F_POSRES].empty())
 718         {
 719             posres_wrapper_lambda(wcycle, fepvals, idef, &pbc_full, x, enerd, lambda, fr);
 720         }
 721         if (idef.ilsort != ilsortNO_FE)
 722         {
 723             wallcycle_sub_start(wcycle, WallCycleSubCounter::ListedFep);
 724             if (idef.ilsort != ilsortFE_SORTED)
 725             {
 726                 gmx_incons("The bonded interactions are not sorted for free energy");
 727             }
 728             for (int i = 0; i < 1 + enerd->foreignLambdaTerms.numLambdas(); i++)
 729             {
 730                 gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, real> lam_i;
 731                 foreignEnergyGroups_->clear();
 732                 std::array<real, F_NRE> foreign_term = { 0 };
 733                 for (auto j : keysOf(lam_i))
 734                 {
 735                     lam_i[j] = (i == 0 ? lambda[static_cast<int>(j)] : fepvals->all_lambda[j][i - 1]);
 736                 }
 737                 calc_listed_lambda(idef,
 738                                    threading_.get(),
 739                                    x,
 740                                    fr,
 741                                    pbc,
 742                                    forceBufferLambda_,
 743                                    shiftForceBufferLambda_,
 744                                    foreignEnergyGroups_.get(),
 745                                    foreign_term,
 746                                    dvdl,
 747                                    nrnb,
 748                                    lam_i,
 749                                    md,
 750                                    fcdata,
 751                                    global_atom_index);
 752                 sum_epot(*foreignEnergyGroups_, foreign_term.data());
 753                 const double dvdlSum = std::accumulate(std::begin(dvdl), std::end(dvdl), 0.);
 754                 std::fill(std::begin(dvdl), std::end(dvdl), 0.0);
 755                 enerd->foreignLambdaTerms.accumulate(i, foreign_term[F_EPOT], dvdlSum);
 756             }
 757             wallcycle_sub_stop(wcycle, WallCycleSubCounter::ListedFep);
 758         }
 759     }
 760 }