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, 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 "listed_forces.h"
  48
  49 #include <cassert>
  50
  51 #include <algorithm>
  52 #include <array>
  53
  54 #include "gromacs/gmxlib/network.h"
  55 #include "gromacs/gmxlib/nrnb.h"
  56 #include "gromacs/listed_forces/bonded.h"
  57 #include "gromacs/listed_forces/disre.h"
  58 #include "gromacs/listed_forces/orires.h"
  59 #include "gromacs/listed_forces/pairs.h"
  60 #include "gromacs/listed_forces/position_restraints.h"
  61 #include "gromacs/math/vec.h"
  62 #include "gromacs/mdlib/enerdata_utils.h"
  63 #include "gromacs/mdlib/force.h"
  64 #include "gromacs/mdtypes/commrec.h"
  65 #include "gromacs/mdtypes/fcdata.h"
  66 #include "gromacs/mdtypes/forcerec.h"
  67 #include "gromacs/mdtypes/inputrec.h"
  68 #include "gromacs/mdtypes/md_enums.h"
  69 #include "gromacs/mdtypes/simulation_workload.h"
  70 #include "gromacs/pbcutil/ishift.h"
  71 #include "gromacs/pbcutil/pbc.h"
  72 #include "gromacs/timing/wallcycle.h"
  73 #include "gromacs/topology/topology.h"
  74 #include "gromacs/utility/exceptions.h"
  75 #include "gromacs/utility/fatalerror.h"
  76 #include "gromacs/utility/smalloc.h"
  77
  78 #include "listed_internal.h"
  79 #include "utilities.h"
  80
  81 namespace
  82 {
  83
  84 using gmx::ArrayRef;
  85
  86 /*! \brief Return true if ftype is an explicit pair-listed LJ or
  87  * COULOMB interaction type: bonded LJ (usually 1-4), or special
  88  * listed non-bonded for FEP. */
  89 bool isPairInteraction(int ftype)
  90 {
  91     return ((ftype) >= F_LJ14 && (ftype) <= F_LJC_PAIRS_NB);
  92 }
  93
  94 /*! \brief Zero thread-local output buffers */
  95 void zero_thread_output(f_thread_t* f_t)
  96 {
  97     constexpr int nelem_fa = sizeof(f_t->f[0]) / sizeof(real);
  98
  99     for (int i = 0; i < f_t->nblock_used; i++)
 100     {
 101         int a0 = f_t->block_index[i] * reduction_block_size;
 102         int a1 = a0 + reduction_block_size;
 103         for (int a = a0; a < a1; a++)
 104         {
 105             for (int d = 0; d < nelem_fa; d++)
 106             {
 107                 f_t->f[a][d] = 0;
 108             }
 109         }
 110     }
 111
 112     for (int i = 0; i < SHIFTS; i++)
 113     {
 114         clear_rvec(f_t->fshift[i]);
 115     }
 116     for (int i = 0; i < F_NRE; i++)
 117     {
 118         f_t->ener[i] = 0;
 119     }
 120     for (int i = 0; i < egNR; i++)
 121     {
 122         for (int j = 0; j < f_t->grpp.nener; j++)
 123         {
 124             f_t->grpp.ener[i][j] = 0;
 125         }
 126     }
 127     for (int i = 0; i < efptNR; i++)
 128     {
 129         f_t->dvdl[i] = 0;
 130     }
 131 }
 132
 133 /*! \brief The max thread number is arbitrary, we used a fixed number
 134  * to avoid memory management.  Using more than 16 threads is probably
 135  * never useful performance wise. */
 136 #define MAX_BONDED_THREADS 256
 137
 138 /*! \brief Reduce thread-local force buffers */
 139 void reduce_thread_forces(int n, gmx::ArrayRef<gmx::RVec> force, const bonded_threading_t* bt, int nthreads)
 140 {
 141     if (nthreads > MAX_BONDED_THREADS)
 142     {
 143         gmx_fatal(FARGS, "Can not reduce bonded forces on more than %d threads", MAX_BONDED_THREADS);
 144     }
 145
 146     rvec* gmx_restrict f = as_rvec_array(force.data());
 147
 148     /* This reduction can run on any number of threads,
 149      * independently of bt->nthreads.
 150      * But if nthreads matches bt->nthreads (which it currently does)
 151      * the uniform distribution of the touched blocks over nthreads will
 152      * match the distribution of bonded over threads well in most cases,
 153      * which means that threads mostly reduce their own data which increases
 154      * the number of cache hits.
 155      */
 156 #pragma omp parallel for num_threads(nthreads) schedule(static)
 157     for (int b = 0; b < bt->nblock_used; b++)
 158     {
 159         try
 160         {
 161             int    ind = bt->block_index[b];
 162             rvec4* fp[MAX_BONDED_THREADS];
 163
 164             /* Determine which threads contribute to this block */
 165             int nfb = 0;
 166             for (int ft = 0; ft < bt->nthreads; ft++)
 167             {
 168                 if (bitmask_is_set(bt->mask[ind], ft))
 169                 {
 170                     fp[nfb++] = bt->f_t[ft]->f;
 171                 }
 172             }
 173             if (nfb > 0)
 174             {
 175                 /* Reduce force buffers for threads that contribute */
 176                 int a0 = ind * reduction_block_size;
 177                 int a1 = (ind + 1) * reduction_block_size;
 178                 /* It would be nice if we could pad f to avoid this min */
 179                 a1 = std::min(a1, n);
 180                 for (int a = a0; a < a1; a++)
 181                 {
 182                     for (int fb = 0; fb < nfb; fb++)
 183                     {
 184                         rvec_inc(f[a], fp[fb][a]);
 185                     }
 186                 }
 187             }
 188         }
 189         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 190     }
 191 }
 192
 193 /*! \brief Reduce thread-local forces, shift forces and energies */
 194 void reduce_thread_output(int                        n,
 195                           gmx::ForceWithShiftForces* forceWithShiftForces,
 196                           real*                      ener,
 197                           gmx_grppairener_t*         grpp,
 198                           real*                      dvdl,
 199                           const bonded_threading_t*  bt,
 200                           const gmx::StepWorkload&   stepWork)
 201 {
 202     assert(bt->haveBondeds);
 203
 204     if (bt->nblock_used > 0)
 205     {
 206         /* Reduce the bonded force buffer */
 207         reduce_thread_forces(n, forceWithShiftForces->force(), bt, bt->nthreads);
 208     }
 209
 210     rvec* gmx_restrict fshift = as_rvec_array(forceWithShiftForces->shiftForces().data());
 211
 212     /* When necessary, reduce energy and virial using one thread only */
 213     if ((stepWork.computeEnergy || stepWork.computeVirial || stepWork.computeDhdl) && bt->nthreads > 1)
 214     {
 215         gmx::ArrayRef<const std::unique_ptr<f_thread_t>> f_t = bt->f_t;
 216
 217         if (stepWork.computeVirial)
 218         {
 219             for (int i = 0; i < SHIFTS; i++)
 220             {
 221                 for (int t = 1; t < bt->nthreads; t++)
 222                 {
 223                     rvec_inc(fshift[i], f_t[t]->fshift[i]);
 224                 }
 225             }
 226         }
 227         if (stepWork.computeEnergy)
 228         {
 229             for (int i = 0; i < F_NRE; i++)
 230             {
 231                 for (int t = 1; t < bt->nthreads; t++)
 232                 {
 233                     ener[i] += f_t[t]->ener[i];
 234                 }
 235             }
 236             for (int i = 0; i < egNR; i++)
 237             {
 238                 for (int j = 0; j < f_t[1]->grpp.nener; j++)
 239                 {
 240                     for (int t = 1; t < bt->nthreads; t++)
 241                     {
 242                         grpp->ener[i][j] += f_t[t]->grpp.ener[i][j];
 243                     }
 244                 }
 245             }
 246         }
 247         if (stepWork.computeDhdl)
 248         {
 249             for (int i = 0; i < efptNR; i++)
 250             {
 251
 252                 for (int t = 1; t < bt->nthreads; t++)
 253                 {
 254                     dvdl[i] += f_t[t]->dvdl[i];
 255                 }
 256             }
 257         }
 258     }
 259 }
 260
 261 /*! \brief Returns the bonded kernel flavor
 262  *
 263  * Note that energies are always requested when the virial
 264  * is requested (performance gain would be small).
 265  * Note that currently we do not have bonded kernels that
 266  * do not compute forces.
 267  */
 268 BondedKernelFlavor selectBondedKernelFlavor(const gmx::StepWorkload& stepWork,
 269                                             const bool               useSimdKernels,
 270                                             const bool               havePerturbedInteractions)
 271 {
 272     BondedKernelFlavor flavor;
 273     if (stepWork.computeEnergy || stepWork.computeVirial)
 274     {
 275         if (stepWork.computeVirial)
 276         {
 277             flavor = BondedKernelFlavor::ForcesAndVirialAndEnergy;
 278         }
 279         else
 280         {
 281             flavor = BondedKernelFlavor::ForcesAndEnergy;
 282         }
 283     }
 284     else
 285     {
 286         if (useSimdKernels && !havePerturbedInteractions)
 287         {
 288             flavor = BondedKernelFlavor::ForcesSimdWhenAvailable;
 289         }
 290         else
 291         {
 292             flavor = BondedKernelFlavor::ForcesNoSimd;
 293         }
 294     }
 295
 296     return flavor;
 297 }
 298
 299 /*! \brief Calculate one element of the list of bonded interactions
 300     for this thread */
 301 real calc_one_bond(int                      thread,
 302                    int                      ftype,
 303                    const t_idef*            idef,
 304                    ArrayRef<const int>      iatoms,
 305                    const int                numNonperturbedInteractions,
 306                    const WorkDivision&      workDivision,
 307                    const rvec               x[],
 308                    rvec4                    f[],
 309                    rvec                     fshift[],
 310                    const t_forcerec*        fr,
 311                    const t_pbc*             pbc,
 312                    const t_graph*           g,
 313                    gmx_grppairener_t*       grpp,
 314                    t_nrnb*                  nrnb,
 315                    const real*              lambda,
 316                    real*                    dvdl,
 317                    const t_mdatoms*         md,
 318                    t_fcdata*                fcd,
 319                    const gmx::StepWorkload& stepWork,
 320                    int*                     global_atom_index)
 321 {
 322     GMX_ASSERT(idef->ilsort == ilsortNO_FE || idef->ilsort == ilsortFE_SORTED,
 323                "The topology should be marked either as no FE or sorted on FE");
 324
 325     const bool havePerturbedInteractions =
 326             (idef->ilsort == ilsortFE_SORTED && numNonperturbedInteractions < iatoms.ssize());
 327     BondedKernelFlavor flavor =
 328             selectBondedKernelFlavor(stepWork, fr->use_simd_kernels, havePerturbedInteractions);
 329     int efptFTYPE;
 330     if (IS_RESTRAINT_TYPE(ftype))
 331     {
 332         efptFTYPE = efptRESTRAINT;
 333     }
 334     else
 335     {
 336         efptFTYPE = efptBONDED;
 337     }
 338
 339     const int nat1   = interaction_function[ftype].nratoms + 1;
 340     const int nbonds = iatoms.ssize() / nat1;
 341
 342     GMX_ASSERT(fr->gpuBonded != nullptr || workDivision.end(ftype) == iatoms.ssize(),
 343                "The thread division should match the topology");
 344
 345     const int nb0 = workDivision.bound(ftype, thread);
 346     const int nbn = workDivision.bound(ftype, thread + 1) - nb0;
 347
 348     real v = 0;
 349     if (!isPairInteraction(ftype))
 350     {
 351         if (ftype == F_CMAP)
 352         {
 353             /* TODO The execution time for CMAP dihedrals might be
 354                nice to account to its own subtimer, but first
 355                wallcycle needs to be extended to support calling from
 356                multiple threads. */
 357             v = cmap_dihs(nbn, iatoms.data() + nb0, idef->iparams, idef->cmap_grid, x, f, fshift,
 358                           pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]), md, fcd, global_atom_index);
 359         }
 360         else
 361         {
 362             v = calculateSimpleBond(ftype, nbn, iatoms.data() + nb0, idef->iparams, x, f, fshift,
 363                                     pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]), md, fcd,
 364                                     global_atom_index, flavor);
 365         }
 366     }
 367     else
 368     {
 369         /* TODO The execution time for pairs might be nice to account
 370            to its own subtimer, but first wallcycle needs to be
 371            extended to support calling from multiple threads. */
 372         do_pairs(ftype, nbn, iatoms.data() + nb0, idef->iparams, x, f, fshift, pbc, g, lambda, dvdl,
 373                  md, fr, havePerturbedInteractions, stepWork, grpp, global_atom_index);
 374     }
 375
 376     if (thread == 0)
 377     {
 378         inc_nrnb(nrnb, nrnbIndex(ftype), nbonds);
 379     }
 380
 381     return v;
 382 }
 383
 384 } // namespace
 385
 386 /*! \brief Compute the bonded part of the listed forces, parallelized over threads
 387  */
 388 static void calcBondedForces(const t_idef*            idef,
 389                              const rvec               x[],
 390                              const t_forcerec*        fr,
 391                              const t_pbc*             pbc_null,
 392                              const t_graph*           g,
 393                              rvec*                    fshiftMasterBuffer,
 394                              gmx_enerdata_t*          enerd,
 395                              t_nrnb*                  nrnb,
 396                              const real*              lambda,
 397                              real*                    dvdl,
 398                              const t_mdatoms*         md,
 399                              t_fcdata*                fcd,
 400                              const gmx::StepWorkload& stepWork,
 401                              int*                     global_atom_index)
 402 {
 403     bonded_threading_t* bt = fr->bondedThreading;
 404
 405 #pragma omp parallel for num_threads(bt->nthreads) schedule(static)
 406     for (int thread = 0; thread < bt->nthreads; thread++)
 407     {
 408         try
 409         {
 410             f_thread_t& threadBuffers = *bt->f_t[thread];
 411             int         ftype;
 412             real *      epot, v;
 413             /* thread stuff */
 414             rvec*              fshift;
 415             real*              dvdlt;
 416             gmx_grppairener_t* grpp;
 417
 418             zero_thread_output(&threadBuffers);
 419
 420             rvec4* ft = threadBuffers.f;
 421
 422             /* Thread 0 writes directly to the main output buffers.
 423              * We might want to reconsider this.
 424              */
 425             if (thread == 0)
 426             {
 427                 fshift = fshiftMasterBuffer;
 428                 epot   = enerd->term;
 429                 grpp   = &enerd->grpp;
 430                 dvdlt  = dvdl;
 431             }
 432             else
 433             {
 434                 fshift = threadBuffers.fshift;
 435                 epot   = threadBuffers.ener;
 436                 grpp   = &threadBuffers.grpp;
 437                 dvdlt  = threadBuffers.dvdl;
 438             }
 439             /* Loop over all bonded force types to calculate the bonded forces */
 440             for (ftype = 0; (ftype < F_NRE); ftype++)
 441             {
 442                 const t_ilist& ilist = idef->il[ftype];
 443                 if (ilist.nr > 0 && ftype_is_bonded_potential(ftype))
 444                 {
 445                     ArrayRef<const int> iatoms = gmx::constArrayRefFromArray(ilist.iatoms, ilist.nr);
 446                     v                          = calc_one_bond(
 447                             thread, ftype, idef, iatoms, idef->numNonperturbedInteractions[ftype],
 448                             fr->bondedThreading->workDivision, x, ft, fshift, fr, pbc_null, g, grpp,
 449                             nrnb, lambda, dvdlt, md, fcd, stepWork, global_atom_index);
 450                     epot[ftype] += v;
 451                 }
 452             }
 453         }
 454         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 455     }
 456 }
 457
 458 bool haveRestraints(const t_idef& idef, const t_fcdata& fcd)
 459 {
 460     return ((idef.il[F_POSRES].nr > 0) || (idef.il[F_FBPOSRES].nr > 0) || fcd.orires.nr > 0
 461             || fcd.disres.nres > 0);
 462 }
 463
 464 bool haveCpuBondeds(const t_forcerec& fr)
 465 {
 466     return fr.bondedThreading->haveBondeds;
 467 }
 468
 469 bool haveCpuListedForces(const t_forcerec& fr, const t_idef& idef, const t_fcdata& fcd)
 470 {
 471     return haveCpuBondeds(fr) || haveRestraints(idef, fcd);
 472 }
 473
 474 void calc_listed(const t_commrec*         cr,
 475                  const gmx_multisim_t*    ms,
 476                  struct gmx_wallcycle*    wcycle,
 477                  const t_idef*            idef,
 478                  const rvec               x[],
 479                  history_t*               hist,
 480                  gmx::ForceOutputs*       forceOutputs,
 481                  const t_forcerec*        fr,
 482                  const struct t_pbc*      pbc,
 483                  const struct t_pbc*      pbc_full,
 484                  const struct t_graph*    g,
 485                  gmx_enerdata_t*          enerd,
 486                  t_nrnb*                  nrnb,
 487                  const real*              lambda,
 488                  const t_mdatoms*         md,
 489                  t_fcdata*                fcd,
 490                  int*                     global_atom_index,
 491                  const gmx::StepWorkload& stepWork)
 492 {
 493     const t_pbc*        pbc_null;
 494     bonded_threading_t* bt = fr->bondedThreading;
 495
 496     if (fr->bMolPBC)
 497     {
 498         pbc_null = pbc;
 499     }
 500     else
 501     {
 502         pbc_null = nullptr;
 503     }
 504
 505     if (haveRestraints(*idef, *fcd))
 506     {
 507         /* TODO Use of restraints triggers further function calls
 508            inside the loop over calc_one_bond(), but those are too
 509            awkward to account to this subtimer properly in the present
 510            code. We don't test / care much about performance with
 511            restraints, anyway. */
 512         wallcycle_sub_start(wcycle, ewcsRESTRAINTS);
 513
 514         if (idef->il[F_POSRES].nr > 0)
 515         {
 516             posres_wrapper(nrnb, idef, pbc_full, x, enerd, lambda, fr, &forceOutputs->forceWithVirial());
 517         }
 518
 519         if (idef->il[F_FBPOSRES].nr > 0)
 520         {
 521             fbposres_wrapper(nrnb, idef, pbc_full, x, enerd, fr, &forceOutputs->forceWithVirial());
 522         }
 523
 524         /* Do pre force calculation stuff which might require communication */
 525         if (fcd->orires.nr > 0)
 526         {
 527             /* This assertion is to ensure we have whole molecules.
 528              * Unfortunately we do not have an mdrun state variable that tells
 529              * us if molecules in x are not broken over PBC, so we have to make
 530              * do with checking graph!=nullptr, which should tell us if we made
 531              * molecules whole before calling the current function.
 532              */
 533             GMX_RELEASE_ASSERT(fr->pbcType == PbcType::No || g != nullptr,
 534                                "With orientation restraints molecules should be whole");
 535             enerd->term[F_ORIRESDEV] = calc_orires_dev(ms, idef->il[F_ORIRES].nr, idef->il[F_ORIRES].iatoms,
 536                                                        idef->iparams, md, x, pbc_null, fcd, hist);
 537         }
 538         if (fcd->disres.nres > 0)
 539         {
 540             calc_disres_R_6(cr, ms, idef->il[F_DISRES].nr, idef->il[F_DISRES].iatoms, x, pbc_null,
 541                             fcd, hist);
 542         }
 543
 544         wallcycle_sub_stop(wcycle, ewcsRESTRAINTS);
 545     }
 546
 547     if (haveCpuBondeds(*fr))
 548     {
 549         gmx::ForceWithShiftForces& forceWithShiftForces = forceOutputs->forceWithShiftForces();
 550
 551         wallcycle_sub_start(wcycle, ewcsLISTED);
 552         /* The dummy array is to have a place to store the dhdl at other values
 553            of lambda, which will be thrown away in the end */
 554         real dvdl[efptNR] = { 0 };
 555         calcBondedForces(idef, x, fr, pbc_null, g,
 556                          as_rvec_array(forceWithShiftForces.shiftForces().data()), enerd, nrnb,
 557                          lambda, dvdl, md, fcd, stepWork, global_atom_index);
 558         wallcycle_sub_stop(wcycle, ewcsLISTED);
 559
 560         wallcycle_sub_start(wcycle, ewcsLISTED_BUF_OPS);
 561         reduce_thread_output(fr->natoms_force, &forceWithShiftForces, enerd->term, &enerd->grpp,
 562                              dvdl, bt, stepWork);
 563
 564         if (stepWork.computeDhdl)
 565         {
 566             for (int i = 0; i < efptNR; i++)
 567             {
 568                 enerd->dvdl_nonlin[i] += dvdl[i];
 569             }
 570         }
 571         wallcycle_sub_stop(wcycle, ewcsLISTED_BUF_OPS);
 572     }
 573
 574     /* Copy the sum of violations for the distance restraints from fcd */
 575     if (fcd)
 576     {
 577         enerd->term[F_DISRESVIOL] = fcd->disres.sumviol;
 578     }
 579 }
 580
 581 void calc_listed_lambda(const t_idef*         idef,
 582                         const rvec            x[],
 583                         const t_forcerec*     fr,
 584                         const struct t_pbc*   pbc,
 585                         const struct t_graph* g,
 586                         gmx_grppairener_t*    grpp,
 587                         real*                 epot,
 588                         t_nrnb*               nrnb,
 589                         const real*           lambda,
 590                         const t_mdatoms*      md,
 591                         t_fcdata*             fcd,
 592                         int*                  global_atom_index)
 593 {
 594     real          v;
 595     real          dvdl_dum[efptNR] = { 0 };
 596     rvec4*        f;
 597     rvec*         fshift;
 598     const t_pbc*  pbc_null;
 599     t_idef        idef_fe;
 600     WorkDivision& workDivision = fr->bondedThreading->foreignLambdaWorkDivision;
 601
 602     if (fr->bMolPBC)
 603     {
 604         pbc_null = pbc;
 605     }
 606     else
 607     {
 608         pbc_null = nullptr;
 609     }
 610
 611     /* Copy the whole idef, so we can modify the contents locally */
 612     idef_fe = *idef;
 613
 614     /* We already have the forces, so we use temp buffers here */
 615     // TODO: Get rid of these allocations by using permanent force buffers
 616     snew(f, fr->natoms_force);
 617     snew(fshift, SHIFTS);
 618
 619     /* Loop over all bonded force types to calculate the bonded energies */
 620     for (int ftype = 0; (ftype < F_NRE); ftype++)
 621     {
 622         if (ftype_is_bonded_potential(ftype))
 623         {
 624             const t_ilist& ilist = idef->il[ftype];
 625             /* Create a temporary iatom list with only perturbed interactions */
 626             const int           numNonperturbed = idef->numNonperturbedInteractions[ftype];
 627             ArrayRef<const int> iatoms = gmx::constArrayRefFromArray(ilist.iatoms + numNonperturbed,
 628                                                                      ilist.nr - numNonperturbed);
 629             t_ilist&            ilist_fe = idef_fe.il[ftype];
 630             /* Set the work range of thread 0 to the perturbed bondeds */
 631             workDivision.setBound(ftype, 0, 0);
 632             workDivision.setBound(ftype, 1, iatoms.ssize());
 633
 634             if (ilist_fe.nr > 0)
 635             {
 636                 gmx::StepWorkload tempFlags;
 637                 tempFlags.computeEnergy = true;
 638                 v = calc_one_bond(0, ftype, idef, iatoms, iatoms.ssize(), workDivision, x, f,
 639                                   fshift, fr, pbc_null, g, grpp, nrnb, lambda, dvdl_dum, md, fcd,
 640                                   tempFlags, global_atom_index);
 641                 epot[ftype] += v;
 642             }
 643         }
 644     }
 645
 646     sfree(fshift);
 647     sfree(f);
 648 }
 649
 650 void do_force_listed(struct gmx_wallcycle*    wcycle,
 651                      const matrix             box,
 652                      const t_lambda*          fepvals,
 653                      const t_commrec*         cr,
 654                      const gmx_multisim_t*    ms,
 655                      const t_idef*            idef,
 656                      const rvec               x[],
 657                      history_t*               hist,
 658                      gmx::ForceOutputs*       forceOutputs,
 659                      const t_forcerec*        fr,
 660                      const struct t_pbc*      pbc,
 661                      const struct t_graph*    graph,
 662                      gmx_enerdata_t*          enerd,
 663                      t_nrnb*                  nrnb,
 664                      const real*              lambda,
 665                      const t_mdatoms*         md,
 666                      t_fcdata*                fcd,
 667                      int*                     global_atom_index,
 668                      const gmx::StepWorkload& stepWork)
 669 {
 670     t_pbc pbc_full; /* Full PBC is needed for position restraints */
 671
 672     if (!stepWork.computeListedForces)
 673     {
 674         return;
 675     }
 676
 677     if ((idef->il[F_POSRES].nr > 0) || (idef->il[F_FBPOSRES].nr > 0))
 678     {
 679         /* Not enough flops to bother counting */
 680         set_pbc(&pbc_full, fr->pbcType, box);
 681     }
 682     calc_listed(cr, ms, wcycle, idef, x, hist, forceOutputs, fr, pbc, &pbc_full, graph, enerd, nrnb,
 683                 lambda, md, fcd, global_atom_index, stepWork);
 684
 685     /* Check if we have to determine energy differences
 686      * at foreign lambda's.
 687      */
 688     if (fepvals->n_lambda > 0 && stepWork.computeDhdl)
 689     {
 690         posres_wrapper_lambda(wcycle, fepvals, idef, &pbc_full, x, enerd, lambda, fr);
 691
 692         if (idef->ilsort != ilsortNO_FE)
 693         {
 694             wallcycle_sub_start(wcycle, ewcsLISTED_FEP);
 695             if (idef->ilsort != ilsortFE_SORTED)
 696             {
 697                 gmx_incons("The bonded interactions are not sorted for free energy");
 698             }
 699             for (size_t i = 0; i < enerd->enerpart_lambda.size(); i++)
 700             {
 701                 real lam_i[efptNR];
 702
 703                 reset_foreign_enerdata(enerd);
 704                 for (int j = 0; j < efptNR; j++)
 705                 {
 706                     lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i - 1]);
 707                 }
 708                 calc_listed_lambda(idef, x, fr, pbc, graph, &(enerd->foreign_grpp),
 709                                    enerd->foreign_term, nrnb, lam_i, md, fcd, global_atom_index);
 710                 sum_epot(&(enerd->foreign_grpp), enerd->foreign_term);
 711                 enerd->enerpart_lambda[i] += enerd->foreign_term[F_EPOT];
 712             }
 713             wallcycle_sub_stop(wcycle, ewcsLISTED_FEP);
 714         }
 715     }
 716 }