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, 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     nat1      = interaction_function[ftype].nratoms + 1;
 291     nbonds    = idef->il[ftype].nr/nat1;
 292     iatoms    = idef->il[ftype].iatoms;
 293
 294     nb0 = idef->il_thread_division[ftype*(idef->nthreads+1)+thread];
 295     nbn = idef->il_thread_division[ftype*(idef->nthreads+1)+thread+1] - nb0;
 296
 297     if (!isPairInteraction(ftype))
 298     {
 299         if (ftype == F_CMAP)
 300         {
 301             /* TODO The execution time for CMAP dihedrals might be
 302                nice to account to its own subtimer, but first
 303                wallcycle needs to be extended to support calling from
 304                multiple threads. */
 305             v = cmap_dihs(nbn, iatoms+nb0,
 306                           idef->iparams, &idef->cmap_grid,
 307                           x, f, fshift,
 308                           pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]),
 309                           md, fcd, global_atom_index);
 310         }
 311 #if GMX_SIMD_HAVE_REAL
 312         else if (ftype == F_ANGLES && bUseSIMD &&
 313                  !bCalcEnerVir && fr->efep == efepNO)
 314         {
 315             /* No energies, shift forces, dvdl */
 316             angles_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
 317                                idef->iparams,
 318                                x, f,
 319                                pbc, g, lambda[efptFTYPE], md, fcd,
 320                                global_atom_index);
 321             v = 0;
 322         }
 323
 324         else if (ftype == F_UREY_BRADLEY && bUseSIMD &&
 325                  !bCalcEnerVir && fr->efep == efepNO)
 326         {
 327             /* No energies, shift forces, dvdl */
 328             urey_bradley_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
 329                                      idef->iparams,
 330                                      x, f,
 331                                      pbc, g, lambda[efptFTYPE], md, fcd,
 332                                      global_atom_index);
 333             v = 0;
 334         }
 335 #endif
 336         else if (ftype == F_PDIHS &&
 337                  !bCalcEnerVir && fr->efep == efepNO)
 338         {
 339             /* No energies, shift forces, dvdl */
 340 #if GMX_SIMD_HAVE_REAL
 341             if (bUseSIMD)
 342             {
 343                 pdihs_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
 344                                   idef->iparams,
 345                                   x, f,
 346                                   pbc, g, lambda[efptFTYPE], md, fcd,
 347                                   global_atom_index);
 348             }
 349             else
 350 #endif
 351             {
 352                 pdihs_noener(nbn, idef->il[ftype].iatoms+nb0,
 353                              idef->iparams,
 354                              x, f,
 355                              pbc, g, lambda[efptFTYPE], md, fcd,
 356                              global_atom_index);
 357             }
 358             v = 0;
 359         }
 360 #if GMX_SIMD_HAVE_REAL
 361         else if (ftype == F_RBDIHS && bUseSIMD &&
 362                  !bCalcEnerVir && fr->efep == efepNO)
 363         {
 364             /* No energies, shift forces, dvdl */
 365             rbdihs_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
 366                                idef->iparams,
 367                                (const rvec*)x, f,
 368                                pbc, g, lambda[efptFTYPE], md, fcd,
 369                                global_atom_index);
 370             v = 0;
 371         }
 372 #endif
 373         else
 374         {
 375             v = interaction_function[ftype].ifunc(nbn, iatoms+nb0,
 376                                                   idef->iparams,
 377                                                   x, f, fshift,
 378                                                   pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]),
 379                                                   md, fcd, global_atom_index);
 380         }
 381     }
 382     else
 383     {
 384         /* TODO The execution time for pairs might be nice to account
 385            to its own subtimer, but first wallcycle needs to be
 386            extended to support calling from multiple threads. */
 387         do_pairs(ftype, nbn, iatoms+nb0, idef->iparams, x, f, fshift,
 388                  pbc, g, lambda, dvdl, md, fr,
 389                  bCalcEnerVir, grpp, global_atom_index);
 390         v = 0;
 391     }
 392
 393     if (thread == 0)
 394     {
 395         inc_nrnb(nrnb, interaction_function[ftype].nrnb_ind, nbonds);
 396     }
 397
 398     return v;
 399 }
 400
 401 } // namespace
 402
 403 gmx_bool
 404 ftype_is_bonded_potential(int ftype)
 405 {
 406     return
 407         (interaction_function[ftype].flags & IF_BOND) &&
 408         !(ftype == F_CONNBONDS || ftype == F_POSRES || ftype == F_FBPOSRES) &&
 409         (ftype < F_GB12 || ftype > F_GB14);
 410 }
 411
 412 /*! \brief Compute the bonded part of the listed forces, parallelized over threads
 413  */
 414 static void
 415 calcBondedForces(const t_idef     *idef,
 416                  const rvec        x[],
 417                  const t_forcerec *fr,
 418                  const t_pbc      *pbc_null,
 419                  const t_graph    *g,
 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                  gmx_bool          bCalcEnerVir,
 427                  int              *global_atom_index)
 428 {
 429     struct bonded_threading_t *bt = fr->bonded_threading;
 430
 431 #pragma omp parallel for num_threads(bt->nthreads) schedule(static)
 432     for (int thread = 0; thread < bt->nthreads; thread++)
 433     {
 434         try
 435         {
 436             int                ftype;
 437             real              *epot, v;
 438             /* thread stuff */
 439             rvec4             *ft;
 440             rvec              *fshift;
 441             real              *dvdlt;
 442             gmx_grppairener_t *grpp;
 443
 444             zero_thread_output(bt, thread);
 445
 446             ft = bt->f_t[thread].f;
 447
 448             if (thread == 0)
 449             {
 450                 fshift = fr->fshift;
 451                 epot   = enerd->term;
 452                 grpp   = &enerd->grpp;
 453                 dvdlt  = dvdl;
 454             }
 455             else
 456             {
 457                 fshift = bt->f_t[thread].fshift;
 458                 epot   = bt->f_t[thread].ener;
 459                 grpp   = &bt->f_t[thread].grpp;
 460                 dvdlt  = bt->f_t[thread].dvdl;
 461             }
 462             /* Loop over all bonded force types to calculate the bonded forces */
 463             for (ftype = 0; (ftype < F_NRE); ftype++)
 464             {
 465                 if (idef->il[ftype].nr > 0 && ftype_is_bonded_potential(ftype))
 466                 {
 467                     v = calc_one_bond(thread, ftype, idef, x,
 468                                       ft, fshift, fr, pbc_null, g, grpp,
 469                                       nrnb, lambda, dvdlt,
 470                                       md, fcd, bCalcEnerVir,
 471                                       global_atom_index);
 472                     epot[ftype] += v;
 473                 }
 474             }
 475         }
 476         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 477     }
 478 }
 479
 480 void calc_listed(const t_commrec             *cr,
 481                  struct gmx_wallcycle        *wcycle,
 482                  const t_idef *idef,
 483                  const rvec x[], history_t *hist,
 484                  rvec f[],
 485                  gmx::ForceWithVirial *forceWithVirial,
 486                  const t_forcerec *fr,
 487                  const struct t_pbc *pbc,
 488                  const struct t_pbc *pbc_full,
 489                  const struct t_graph *g,
 490                  gmx_enerdata_t *enerd, t_nrnb *nrnb,
 491                  const real *lambda,
 492                  const t_mdatoms *md,
 493                  t_fcdata *fcd, int *global_atom_index,
 494                  int force_flags)
 495 {
 496     struct bonded_threading_t *bt;
 497     gmx_bool                   bCalcEnerVir;
 498     const  t_pbc              *pbc_null;
 499
 500     bt = fr->bonded_threading;
 501
 502     assert(bt->nthreads == idef->nthreads);
 503
 504     bCalcEnerVir = (force_flags & (GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY));
 505
 506     if (fr->bMolPBC)
 507     {
 508         pbc_null = pbc;
 509     }
 510     else
 511     {
 512         pbc_null = nullptr;
 513     }
 514
 515     if ((idef->il[F_POSRES].nr > 0) ||
 516         (idef->il[F_FBPOSRES].nr > 0) ||
 517         fcd->orires.nr > 0 ||
 518         fcd->disres.nres > 0)
 519     {
 520         /* TODO Use of restraints triggers further function calls
 521            inside the loop over calc_one_bond(), but those are too
 522            awkward to account to this subtimer properly in the present
 523            code. We don't test / care much about performance with
 524            restraints, anyway. */
 525         wallcycle_sub_start(wcycle, ewcsRESTRAINTS);
 526
 527         if (idef->il[F_POSRES].nr > 0)
 528         {
 529             posres_wrapper(nrnb, idef, pbc_full, x, enerd, lambda, fr,
 530                            forceWithVirial);
 531         }
 532
 533         if (idef->il[F_FBPOSRES].nr > 0)
 534         {
 535             fbposres_wrapper(nrnb, idef, pbc_full, x, enerd, fr,
 536                              forceWithVirial);
 537         }
 538
 539         /* Do pre force calculation stuff which might require communication */
 540         if (fcd->orires.nr > 0)
 541         {
 542             /* This assertion is to ensure we have whole molecules.
 543              * Unfortunately we do not have an mdrun state variable that tells
 544              * us if molecules in x are not broken over PBC, so we have to make
 545              * do with checking graph!=nullptr, which should tell us if we made
 546              * molecules whole before calling the current function.
 547              */
 548             GMX_RELEASE_ASSERT(fr->ePBC == epbcNONE || g != nullptr, "With orientation restraints molecules should be whole");
 549             enerd->term[F_ORIRESDEV] =
 550                 calc_orires_dev(cr->ms, idef->il[F_ORIRES].nr,
 551                                 idef->il[F_ORIRES].iatoms,
 552                                 idef->iparams, md, x,
 553                                 pbc_null, fcd, hist);
 554         }
 555         if (fcd->disres.nres > 0)
 556         {
 557             calc_disres_R_6(cr,
 558                             idef->il[F_DISRES].nr,
 559                             idef->il[F_DISRES].iatoms,
 560                             x, pbc_null,
 561                             fcd, hist);
 562         }
 563
 564         wallcycle_sub_stop(wcycle, ewcsRESTRAINTS);
 565     }
 566
 567     if (bt->haveBondeds)
 568     {
 569         wallcycle_sub_start(wcycle, ewcsLISTED);
 570         /* The dummy array is to have a place to store the dhdl at other values
 571            of lambda, which will be thrown away in the end */
 572         real dvdl[efptNR] = {0};
 573         calcBondedForces(idef, x, fr, pbc_null, g, enerd, nrnb, lambda, dvdl, md,
 574                          fcd, bCalcEnerVir, global_atom_index);
 575         wallcycle_sub_stop(wcycle, ewcsLISTED);
 576
 577         wallcycle_sub_start(wcycle, ewcsLISTED_BUF_OPS);
 578         reduce_thread_output(fr->natoms_force, f, fr->fshift,
 579                              enerd->term, &enerd->grpp, dvdl,
 580                              bt,
 581                              bCalcEnerVir,
 582                              force_flags & GMX_FORCE_DHDL);
 583
 584         if (force_flags & GMX_FORCE_DHDL)
 585         {
 586             for (int i = 0; i < efptNR; i++)
 587             {
 588                 enerd->dvdl_nonlin[i] += dvdl[i];
 589             }
 590         }
 591         wallcycle_sub_stop(wcycle, ewcsLISTED_BUF_OPS);
 592     }
 593
 594     /* Copy the sum of violations for the distance restraints from fcd */
 595     if (fcd)
 596     {
 597         enerd->term[F_DISRESVIOL] = fcd->disres.sumviol;
 598     }
 599 }
 600
 601 void calc_listed_lambda(const t_idef *idef,
 602                         const rvec x[],
 603                         const t_forcerec *fr,
 604                         const struct t_pbc *pbc, const struct t_graph *g,
 605                         gmx_grppairener_t *grpp, real *epot, t_nrnb *nrnb,
 606                         const real *lambda,
 607                         const t_mdatoms *md,
 608                         t_fcdata *fcd,
 609                         int *global_atom_index)
 610 {
 611     int           ftype, nr_nonperturbed, nr;
 612     real          v;
 613     real          dvdl_dum[efptNR] = {0};
 614     rvec4        *f;
 615     rvec         *fshift;
 616     const  t_pbc *pbc_null;
 617     t_idef        idef_fe;
 618
 619     if (fr->bMolPBC)
 620     {
 621         pbc_null = pbc;
 622     }
 623     else
 624     {
 625         pbc_null = nullptr;
 626     }
 627
 628     /* Copy the whole idef, so we can modify the contents locally */
 629     idef_fe          = *idef;
 630     idef_fe.nthreads = 1;
 631     snew(idef_fe.il_thread_division, F_NRE*(idef_fe.nthreads+1));
 632
 633     /* We already have the forces, so we use temp buffers here */
 634     snew(f, fr->natoms_force);
 635     snew(fshift, SHIFTS);
 636
 637     /* Loop over all bonded force types to calculate the bonded energies */
 638     for (ftype = 0; (ftype < F_NRE); ftype++)
 639     {
 640         if (ftype_is_bonded_potential(ftype))
 641         {
 642             /* Set the work range of thread 0 to the perturbed bondeds only */
 643             nr_nonperturbed                       = idef->il[ftype].nr_nonperturbed;
 644             nr                                    = idef->il[ftype].nr;
 645             idef_fe.il_thread_division[ftype*2+0] = nr_nonperturbed;
 646             idef_fe.il_thread_division[ftype*2+1] = nr;
 647
 648             /* This is only to get the flop count correct */
 649             idef_fe.il[ftype].nr = nr - nr_nonperturbed;
 650
 651             if (nr - nr_nonperturbed > 0)
 652             {
 653                 v = calc_one_bond(0, ftype, &idef_fe,
 654                                   x, f, fshift, fr, pbc_null, g,
 655                                   grpp, nrnb, lambda, dvdl_dum,
 656                                   md, fcd, TRUE,
 657                                   global_atom_index);
 658                 epot[ftype] += v;
 659             }
 660         }
 661     }
 662
 663     sfree(fshift);
 664     sfree(f);
 665
 666     sfree(idef_fe.il_thread_division);
 667 }
 668
 669 void
 670 do_force_listed(struct gmx_wallcycle        *wcycle,
 671                 matrix                       box,
 672                 const t_lambda              *fepvals,
 673                 const t_commrec             *cr,
 674                 const t_idef                *idef,
 675                 const rvec                   x[],
 676                 history_t                   *hist,
 677                 rvec                        *forceForUseWithShiftForces,
 678                 gmx::ForceWithVirial        *forceWithVirial,
 679                 const t_forcerec            *fr,
 680                 const struct t_pbc          *pbc,
 681                 const struct t_graph        *graph,
 682                 gmx_enerdata_t              *enerd,
 683                 t_nrnb                      *nrnb,
 684                 const real                  *lambda,
 685                 const t_mdatoms             *md,
 686                 t_fcdata                    *fcd,
 687                 int                         *global_atom_index,
 688                 int                          flags)
 689 {
 690     t_pbc pbc_full; /* Full PBC is needed for position restraints */
 691
 692     if (!(flags & GMX_FORCE_LISTED))
 693     {
 694         return;
 695     }
 696
 697     if ((idef->il[F_POSRES].nr > 0) ||
 698         (idef->il[F_FBPOSRES].nr > 0))
 699     {
 700         /* Not enough flops to bother counting */
 701         set_pbc(&pbc_full, fr->ePBC, box);
 702     }
 703     calc_listed(cr, wcycle, idef, x, hist,
 704                 forceForUseWithShiftForces, forceWithVirial,
 705                 fr, pbc, &pbc_full,
 706                 graph, enerd, nrnb, lambda, md, fcd,
 707                 global_atom_index, flags);
 708
 709     /* Check if we have to determine energy differences
 710      * at foreign lambda's.
 711      */
 712     if (fepvals->n_lambda > 0 && (flags & GMX_FORCE_DHDL))
 713     {
 714         posres_wrapper_lambda(wcycle, fepvals, idef, &pbc_full, x, enerd, lambda, fr);
 715
 716         if (idef->ilsort != ilsortNO_FE)
 717         {
 718             wallcycle_sub_start(wcycle, ewcsLISTED_FEP);
 719             if (idef->ilsort != ilsortFE_SORTED)
 720             {
 721                 gmx_incons("The bonded interactions are not sorted for free energy");
 722             }
 723             for (int i = 0; i < enerd->n_lambda; i++)
 724             {
 725                 real lam_i[efptNR];
 726
 727                 reset_foreign_enerdata(enerd);
 728                 for (int j = 0; j < efptNR; j++)
 729                 {
 730                     lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i-1]);
 731                 }
 732                 calc_listed_lambda(idef, x, fr, pbc, graph, &(enerd->foreign_grpp), enerd->foreign_term, nrnb, lam_i, md,
 733                                    fcd, global_atom_index);
 734                 sum_epot(&(enerd->foreign_grpp), enerd->foreign_term);
 735                 enerd->enerpart_lambda[i] += enerd->foreign_term[F_EPOT];
 736             }
 737             wallcycle_sub_stop(wcycle, ewcsLISTED_FEP);
 738         }
 739     }
 740 }