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