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