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 t_idef*            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                    const t_graph*           g,
 314                    gmx_grppairener_t*       grpp,
 315                    t_nrnb*                  nrnb,
 316                    const real*              lambda,
 317                    real*                    dvdl,
 318                    const t_mdatoms*         md,
 319                    t_fcdata*                fcd,
 320                    const gmx::StepWorkload& stepWork,
 321                    int*                     global_atom_index)
 322 {
 323     GMX_ASSERT(idef->ilsort == ilsortNO_FE || idef->ilsort == ilsortFE_SORTED,
 324                "The topology should be marked either as no FE or sorted on FE");
 325
 326     const bool havePerturbedInteractions =
 327             (idef->ilsort == ilsortFE_SORTED && numNonperturbedInteractions < iatoms.ssize());
 328     BondedKernelFlavor flavor =
 329             selectBondedKernelFlavor(stepWork, fr->use_simd_kernels, havePerturbedInteractions);
 330     int efptFTYPE;
 331     if (IS_RESTRAINT_TYPE(ftype))
 332     {
 333         efptFTYPE = efptRESTRAINT;
 334     }
 335     else
 336     {
 337         efptFTYPE = efptBONDED;
 338     }
 339
 340     const int nat1   = interaction_function[ftype].nratoms + 1;
 341     const int nbonds = iatoms.ssize() / nat1;
 342
 343     GMX_ASSERT(fr->gpuBonded != nullptr || workDivision.end(ftype) == iatoms.ssize(),
 344                "The thread division should match the topology");
 345
 346     const int nb0 = workDivision.bound(ftype, thread);
 347     const int nbn = workDivision.bound(ftype, thread + 1) - nb0;
 348
 349     real v = 0;
 350     if (!isPairInteraction(ftype))
 351     {
 352         if (ftype == F_CMAP)
 353         {
 354             /* TODO The execution time for CMAP dihedrals might be
 355                nice to account to its own subtimer, but first
 356                wallcycle needs to be extended to support calling from
 357                multiple threads. */
 358             v = cmap_dihs(nbn, iatoms.data() + nb0, idef->iparams, idef->cmap_grid, x, f, fshift,
 359                           pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]), md, fcd, global_atom_index);
 360         }
 361         else
 362         {
 363             v = calculateSimpleBond(ftype, nbn, iatoms.data() + nb0, idef->iparams, x, f, fshift,
 364                                     pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]), md, fcd,
 365                                     global_atom_index, flavor);
 366         }
 367     }
 368     else
 369     {
 370         /* TODO The execution time for pairs might be nice to account
 371            to its own subtimer, but first wallcycle needs to be
 372            extended to support calling from multiple threads. */
 373         do_pairs(ftype, nbn, iatoms.data() + nb0, idef->iparams, x, f, fshift, pbc, g, lambda, dvdl,
 374                  md, fr, havePerturbedInteractions, stepWork, grpp, global_atom_index);
 375     }
 376
 377     if (thread == 0)
 378     {
 379         inc_nrnb(nrnb, nrnbIndex(ftype), nbonds);
 380     }
 381
 382     return v;
 383 }
 384
 385 } // namespace
 386
 387 /*! \brief Compute the bonded part of the listed forces, parallelized over threads
 388  */
 389 static void calcBondedForces(const t_idef*            idef,
 390                              const rvec               x[],
 391                              const t_forcerec*        fr,
 392                              const t_pbc*             pbc_null,
 393                              const t_graph*           g,
 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 t_ilist& ilist = idef->il[ftype];
 444                 if (ilist.nr > 0 && ftype_is_bonded_potential(ftype))
 445                 {
 446                     ArrayRef<const int> iatoms = gmx::constArrayRefFromArray(ilist.iatoms, ilist.nr);
 447                     v                          = calc_one_bond(
 448                             thread, ftype, idef, iatoms, idef->numNonperturbedInteractions[ftype],
 449                             fr->bondedThreading->workDivision, x, ft, fshift, fr, pbc_null, g, grpp,
 450                             nrnb, lambda, dvdlt, md, fcd, stepWork, global_atom_index);
 451                     epot[ftype] += v;
 452                 }
 453             }
 454         }
 455         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 456     }
 457 }
 458
 459 bool haveRestraints(const t_idef& idef, const t_fcdata& fcd)
 460 {
 461     return ((idef.il[F_POSRES].nr > 0) || (idef.il[F_FBPOSRES].nr > 0) || fcd.orires.nr > 0
 462             || fcd.disres.nres > 0);
 463 }
 464
 465 bool haveCpuBondeds(const t_forcerec& fr)
 466 {
 467     return fr.bondedThreading->haveBondeds;
 468 }
 469
 470 bool haveCpuListedForces(const t_forcerec& fr, const t_idef& idef, const t_fcdata& fcd)
 471 {
 472     return haveCpuBondeds(fr) || haveRestraints(idef, fcd);
 473 }
 474
 475 void calc_listed(const t_commrec*         cr,
 476                  const gmx_multisim_t*    ms,
 477                  struct gmx_wallcycle*    wcycle,
 478                  const t_idef*            idef,
 479                  const rvec               x[],
 480                  history_t*               hist,
 481                  gmx::ForceOutputs*       forceOutputs,
 482                  const t_forcerec*        fr,
 483                  const struct t_pbc*      pbc,
 484                  const struct t_pbc*      pbc_full,
 485                  const struct t_graph*    g,
 486                  gmx_enerdata_t*          enerd,
 487                  t_nrnb*                  nrnb,
 488                  const real*              lambda,
 489                  const t_mdatoms*         md,
 490                  t_fcdata*                fcd,
 491                  int*                     global_atom_index,
 492                  const gmx::StepWorkload& stepWork)
 493 {
 494     const t_pbc*        pbc_null;
 495     bonded_threading_t* bt = fr->bondedThreading;
 496
 497     if (fr->bMolPBC)
 498     {
 499         pbc_null = pbc;
 500     }
 501     else
 502     {
 503         pbc_null = nullptr;
 504     }
 505
 506     if (haveRestraints(*idef, *fcd))
 507     {
 508         /* TODO Use of restraints triggers further function calls
 509            inside the loop over calc_one_bond(), but those are too
 510            awkward to account to this subtimer properly in the present
 511            code. We don't test / care much about performance with
 512            restraints, anyway. */
 513         wallcycle_sub_start(wcycle, ewcsRESTRAINTS);
 514
 515         if (idef->il[F_POSRES].nr > 0)
 516         {
 517             posres_wrapper(nrnb, idef, pbc_full, x, enerd, lambda, fr, &forceOutputs->forceWithVirial());
 518         }
 519
 520         if (idef->il[F_FBPOSRES].nr > 0)
 521         {
 522             fbposres_wrapper(nrnb, idef, pbc_full, x, enerd, fr, &forceOutputs->forceWithVirial());
 523         }
 524
 525         /* Do pre force calculation stuff which might require communication */
 526         if (fcd->orires.nr > 0)
 527         {
 528             /* This assertion is to ensure we have whole molecules.
 529              * Unfortunately we do not have an mdrun state variable that tells
 530              * us if molecules in x are not broken over PBC, so we have to make
 531              * do with checking graph!=nullptr, which should tell us if we made
 532              * molecules whole before calling the current function.
 533              */
 534             GMX_RELEASE_ASSERT(fr->pbcType == PbcType::No || g != nullptr,
 535                                "With orientation restraints molecules should be whole");
 536             enerd->term[F_ORIRESDEV] = calc_orires_dev(ms, idef->il[F_ORIRES].nr, idef->il[F_ORIRES].iatoms,
 537                                                        idef->iparams, md, x, pbc_null, fcd, hist);
 538         }
 539         if (fcd->disres.nres > 0)
 540         {
 541             calc_disres_R_6(cr, ms, idef->il[F_DISRES].nr, idef->il[F_DISRES].iatoms, x, pbc_null,
 542                             fcd, hist);
 543         }
 544
 545         wallcycle_sub_stop(wcycle, ewcsRESTRAINTS);
 546     }
 547
 548     if (haveCpuBondeds(*fr))
 549     {
 550         gmx::ForceWithShiftForces& forceWithShiftForces = forceOutputs->forceWithShiftForces();
 551
 552         wallcycle_sub_start(wcycle, ewcsLISTED);
 553         /* The dummy array is to have a place to store the dhdl at other values
 554            of lambda, which will be thrown away in the end */
 555         real dvdl[efptNR] = { 0 };
 556         calcBondedForces(idef, x, fr, pbc_null, g,
 557                          as_rvec_array(forceWithShiftForces.shiftForces().data()), enerd, nrnb,
 558                          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 void calc_listed_lambda(const t_idef*         idef,
 583                         const rvec            x[],
 584                         const t_forcerec*     fr,
 585                         const struct t_pbc*   pbc,
 586                         const struct t_graph* g,
 587                         gmx_grppairener_t*    grpp,
 588                         real*                 epot,
 589                         gmx::ArrayRef<real>   dvdl,
 590                         t_nrnb*               nrnb,
 591                         const real*           lambda,
 592                         const t_mdatoms*      md,
 593                         t_fcdata*             fcd,
 594                         int*                  global_atom_index)
 595 {
 596     real          v;
 597     rvec4*        f;
 598     rvec*         fshift;
 599     const t_pbc*  pbc_null;
 600     t_idef        idef_fe;
 601     WorkDivision& workDivision = fr->bondedThreading->foreignLambdaWorkDivision;
 602
 603     if (fr->bMolPBC)
 604     {
 605         pbc_null = pbc;
 606     }
 607     else
 608     {
 609         pbc_null = nullptr;
 610     }
 611
 612     /* Copy the whole idef, so we can modify the contents locally */
 613     idef_fe = *idef;
 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 t_ilist& 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> iatoms = gmx::constArrayRefFromArray(ilist.iatoms + numNonperturbed,
 629                                                                      ilist.nr - numNonperturbed);
 630             t_ilist&            ilist_fe = idef_fe.il[ftype];
 631             /* Set the work range of thread 0 to the perturbed bondeds */
 632             workDivision.setBound(ftype, 0, 0);
 633             workDivision.setBound(ftype, 1, iatoms.ssize());
 634
 635             if (ilist_fe.nr > 0)
 636             {
 637                 gmx::StepWorkload tempFlags;
 638                 tempFlags.computeEnergy = true;
 639                 v = calc_one_bond(0, ftype, idef, iatoms, iatoms.ssize(), workDivision, x, f,
 640                                   fshift, fr, pbc_null, g, grpp, nrnb, lambda, dvdl.data(), md, fcd,
 641                                   tempFlags, global_atom_index);
 642                 epot[ftype] += v;
 643             }
 644         }
 645     }
 646
 647     sfree(fshift);
 648     sfree(f);
 649 }
 650
 651 void do_force_listed(struct gmx_wallcycle*    wcycle,
 652                      const matrix             box,
 653                      const t_lambda*          fepvals,
 654                      const t_commrec*         cr,
 655                      const gmx_multisim_t*    ms,
 656                      const t_idef*            idef,
 657                      const rvec               x[],
 658                      history_t*               hist,
 659                      gmx::ForceOutputs*       forceOutputs,
 660                      const t_forcerec*        fr,
 661                      const struct t_pbc*      pbc,
 662                      const struct t_graph*    graph,
 663                      gmx_enerdata_t*          enerd,
 664                      t_nrnb*                  nrnb,
 665                      const real*              lambda,
 666                      const t_mdatoms*         md,
 667                      t_fcdata*                fcd,
 668                      int*                     global_atom_index,
 669                      const gmx::StepWorkload& stepWork)
 670 {
 671     t_pbc pbc_full; /* Full PBC is needed for position restraints */
 672
 673     if (!stepWork.computeListedForces)
 674     {
 675         return;
 676     }
 677
 678     if ((idef->il[F_POSRES].nr > 0) || (idef->il[F_FBPOSRES].nr > 0))
 679     {
 680         /* Not enough flops to bother counting */
 681         set_pbc(&pbc_full, fr->pbcType, box);
 682     }
 683     calc_listed(cr, ms, wcycle, idef, x, hist, forceOutputs, fr, pbc, &pbc_full, graph, enerd, nrnb,
 684                 lambda, md, fcd, global_atom_index, stepWork);
 685
 686     /* Check if we have to determine energy differences
 687      * at foreign lambda's.
 688      */
 689     if (fepvals->n_lambda > 0 && stepWork.computeDhdl)
 690     {
 691         real dvdl[efptNR] = { 0 };
 692         posres_wrapper_lambda(wcycle, fepvals, idef, &pbc_full, x, enerd, lambda, fr);
 693
 694         if (idef->ilsort != ilsortNO_FE)
 695         {
 696             wallcycle_sub_start(wcycle, ewcsLISTED_FEP);
 697             if (idef->ilsort != ilsortFE_SORTED)
 698             {
 699                 gmx_incons("The bonded interactions are not sorted for free energy");
 700             }
 701             for (size_t i = 0; i < enerd->enerpart_lambda.size(); i++)
 702             {
 703                 real lam_i[efptNR];
 704
 705                 reset_foreign_enerdata(enerd);
 706                 for (int j = 0; j < efptNR; j++)
 707                 {
 708                     lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i - 1]);
 709                 }
 710                 calc_listed_lambda(idef, x, fr, pbc, graph, &(enerd->foreign_grpp),
 711                                    enerd->foreign_term, dvdl, nrnb, lam_i, md, fcd, global_atom_index);
 712                 sum_epot(&(enerd->foreign_grpp), enerd->foreign_term);
 713                 enerd->enerpart_lambda[i] += enerd->foreign_term[F_EPOT];
 714                 for (int j = 0; j < efptNR; j++)
 715                 {
 716                     enerd->dhdlLambda[i] += dvdl[j];
 717                     dvdl[j] = 0;
 718                 }
 719             }
 720             wallcycle_sub_stop(wcycle, ewcsLISTED_FEP);
 721         }
 722     }
 723 }