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