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