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