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     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 }
 410
 411 /*! \brief Compute the bonded part of the listed forces, parallelized over threads
 412  */
 413 static void
 414 calcBondedForces(const t_idef     *idef,
 415                  const rvec        x[],
 416                  const t_forcerec *fr,
 417                  const t_pbc      *pbc_null,
 418                  const t_graph    *g,
 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                  gmx_bool          bCalcEnerVir,
 426                  int              *global_atom_index)
 427 {
 428     struct bonded_threading_t *bt = fr->bonded_threading;
 429
 430 #pragma omp parallel for num_threads(bt->nthreads) schedule(static)
 431     for (int thread = 0; thread < bt->nthreads; thread++)
 432     {
 433         try
 434         {
 435             int                ftype;
 436             real              *epot, v;
 437             /* thread stuff */
 438             rvec4             *ft;
 439             rvec              *fshift;
 440             real              *dvdlt;
 441             gmx_grppairener_t *grpp;
 442
 443             zero_thread_output(bt, thread);
 444
 445             ft = bt->f_t[thread].f;
 446
 447             if (thread == 0)
 448             {
 449                 fshift = fr->fshift;
 450                 epot   = enerd->term;
 451                 grpp   = &enerd->grpp;
 452                 dvdlt  = dvdl;
 453             }
 454             else
 455             {
 456                 fshift = bt->f_t[thread].fshift;
 457                 epot   = bt->f_t[thread].ener;
 458                 grpp   = &bt->f_t[thread].grpp;
 459                 dvdlt  = bt->f_t[thread].dvdl;
 460             }
 461             /* Loop over all bonded force types to calculate the bonded forces */
 462             for (ftype = 0; (ftype < F_NRE); ftype++)
 463             {
 464                 if (idef->il[ftype].nr > 0 && ftype_is_bonded_potential(ftype))
 465                 {
 466                     v = calc_one_bond(thread, ftype, idef, x,
 467                                       ft, fshift, fr, pbc_null, g, grpp,
 468                                       nrnb, lambda, dvdlt,
 469                                       md, fcd, bCalcEnerVir,
 470                                       global_atom_index);
 471                     epot[ftype] += v;
 472                 }
 473             }
 474         }
 475         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 476     }
 477 }
 478
 479 void calc_listed(const t_commrec             *cr,
 480                  struct gmx_wallcycle        *wcycle,
 481                  const t_idef *idef,
 482                  const rvec x[], history_t *hist,
 483                  rvec f[],
 484                  gmx::ForceWithVirial *forceWithVirial,
 485                  const t_forcerec *fr,
 486                  const struct t_pbc *pbc,
 487                  const struct t_pbc *pbc_full,
 488                  const struct t_graph *g,
 489                  gmx_enerdata_t *enerd, t_nrnb *nrnb,
 490                  const real *lambda,
 491                  const t_mdatoms *md,
 492                  t_fcdata *fcd, int *global_atom_index,
 493                  int force_flags)
 494 {
 495     struct bonded_threading_t *bt;
 496     gmx_bool                   bCalcEnerVir;
 497     const  t_pbc              *pbc_null;
 498
 499     bt = fr->bonded_threading;
 500
 501     assert(bt->nthreads == idef->nthreads);
 502
 503     bCalcEnerVir = (force_flags & (GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY));
 504
 505     if (fr->bMolPBC)
 506     {
 507         pbc_null = pbc;
 508     }
 509     else
 510     {
 511         pbc_null = nullptr;
 512     }
 513
 514     if ((idef->il[F_POSRES].nr > 0) ||
 515         (idef->il[F_FBPOSRES].nr > 0) ||
 516         fcd->orires.nr > 0 ||
 517         fcd->disres.nres > 0)
 518     {
 519         /* TODO Use of restraints triggers further function calls
 520            inside the loop over calc_one_bond(), but those are too
 521            awkward to account to this subtimer properly in the present
 522            code. We don't test / care much about performance with
 523            restraints, anyway. */
 524         wallcycle_sub_start(wcycle, ewcsRESTRAINTS);
 525
 526         if (idef->il[F_POSRES].nr > 0)
 527         {
 528             posres_wrapper(nrnb, idef, pbc_full, x, enerd, lambda, fr,
 529                            forceWithVirial);
 530         }
 531
 532         if (idef->il[F_FBPOSRES].nr > 0)
 533         {
 534             fbposres_wrapper(nrnb, idef, pbc_full, x, enerd, fr,
 535                              forceWithVirial);
 536         }
 537
 538         /* Do pre force calculation stuff which might require communication */
 539         if (fcd->orires.nr > 0)
 540         {
 541             /* This assertion is to ensure we have whole molecules.
 542              * Unfortunately we do not have an mdrun state variable that tells
 543              * us if molecules in x are not broken over PBC, so we have to make
 544              * do with checking graph!=nullptr, which should tell us if we made
 545              * molecules whole before calling the current function.
 546              */
 547             GMX_RELEASE_ASSERT(fr->ePBC == epbcNONE || g != nullptr, "With orientation restraints molecules should be whole");
 548             enerd->term[F_ORIRESDEV] =
 549                 calc_orires_dev(cr->ms, idef->il[F_ORIRES].nr,
 550                                 idef->il[F_ORIRES].iatoms,
 551                                 idef->iparams, md, x,
 552                                 pbc_null, fcd, hist);
 553         }
 554         if (fcd->disres.nres > 0)
 555         {
 556             calc_disres_R_6(cr,
 557                             idef->il[F_DISRES].nr,
 558                             idef->il[F_DISRES].iatoms,
 559                             x, pbc_null,
 560                             fcd, hist);
 561         }
 562
 563         wallcycle_sub_stop(wcycle, ewcsRESTRAINTS);
 564     }
 565
 566     if (bt->haveBondeds)
 567     {
 568         wallcycle_sub_start(wcycle, ewcsLISTED);
 569         /* The dummy array is to have a place to store the dhdl at other values
 570            of lambda, which will be thrown away in the end */
 571         real dvdl[efptNR] = {0};
 572         calcBondedForces(idef, x, fr, pbc_null, g, enerd, nrnb, lambda, dvdl, md,
 573                          fcd, bCalcEnerVir, global_atom_index);
 574         wallcycle_sub_stop(wcycle, ewcsLISTED);
 575
 576         wallcycle_sub_start(wcycle, ewcsLISTED_BUF_OPS);
 577         reduce_thread_output(fr->natoms_force, f, fr->fshift,
 578                              enerd->term, &enerd->grpp, dvdl,
 579                              bt,
 580                              bCalcEnerVir,
 581                              force_flags & GMX_FORCE_DHDL);
 582
 583         if (force_flags & GMX_FORCE_DHDL)
 584         {
 585             for (int i = 0; i < efptNR; i++)
 586             {
 587                 enerd->dvdl_nonlin[i] += dvdl[i];
 588             }
 589         }
 590         wallcycle_sub_stop(wcycle, ewcsLISTED_BUF_OPS);
 591     }
 592
 593     /* Copy the sum of violations for the distance restraints from fcd */
 594     if (fcd)
 595     {
 596         enerd->term[F_DISRESVIOL] = fcd->disres.sumviol;
 597     }
 598 }
 599
 600 void calc_listed_lambda(const t_idef *idef,
 601                         const rvec x[],
 602                         const t_forcerec *fr,
 603                         const struct t_pbc *pbc, const struct t_graph *g,
 604                         gmx_grppairener_t *grpp, real *epot, t_nrnb *nrnb,
 605                         const real *lambda,
 606                         const t_mdatoms *md,
 607                         t_fcdata *fcd,
 608                         int *global_atom_index)
 609 {
 610     int           ftype, nr_nonperturbed, nr;
 611     real          v;
 612     real          dvdl_dum[efptNR] = {0};
 613     rvec4        *f;
 614     rvec         *fshift;
 615     const  t_pbc *pbc_null;
 616     t_idef        idef_fe;
 617
 618     if (fr->bMolPBC)
 619     {
 620         pbc_null = pbc;
 621     }
 622     else
 623     {
 624         pbc_null = nullptr;
 625     }
 626
 627     /* Copy the whole idef, so we can modify the contents locally */
 628     idef_fe          = *idef;
 629     idef_fe.nthreads = 1;
 630     snew(idef_fe.il_thread_division, F_NRE*(idef_fe.nthreads+1));
 631
 632     /* We already have the forces, so we use temp buffers here */
 633     snew(f, fr->natoms_force);
 634     snew(fshift, SHIFTS);
 635
 636     /* Loop over all bonded force types to calculate the bonded energies */
 637     for (ftype = 0; (ftype < F_NRE); ftype++)
 638     {
 639         if (ftype_is_bonded_potential(ftype))
 640         {
 641             /* Set the work range of thread 0 to the perturbed bondeds only */
 642             nr_nonperturbed                       = idef->il[ftype].nr_nonperturbed;
 643             nr                                    = idef->il[ftype].nr;
 644             idef_fe.il_thread_division[ftype*2+0] = nr_nonperturbed;
 645             idef_fe.il_thread_division[ftype*2+1] = nr;
 646
 647             /* This is only to get the flop count correct */
 648             idef_fe.il[ftype].nr = nr - nr_nonperturbed;
 649
 650             if (nr - nr_nonperturbed > 0)
 651             {
 652                 v = calc_one_bond(0, ftype, &idef_fe,
 653                                   x, f, fshift, fr, pbc_null, g,
 654                                   grpp, nrnb, lambda, dvdl_dum,
 655                                   md, fcd, TRUE,
 656                                   global_atom_index);
 657                 epot[ftype] += v;
 658             }
 659         }
 660     }
 661
 662     sfree(fshift);
 663     sfree(f);
 664
 665     sfree(idef_fe.il_thread_division);
 666 }
 667
 668 void
 669 do_force_listed(struct gmx_wallcycle        *wcycle,
 670                 matrix                       box,
 671                 const t_lambda              *fepvals,
 672                 const t_commrec             *cr,
 673                 const t_idef                *idef,
 674                 const rvec                   x[],
 675                 history_t                   *hist,
 676                 rvec                        *forceForUseWithShiftForces,
 677                 gmx::ForceWithVirial        *forceWithVirial,
 678                 const t_forcerec            *fr,
 679                 const struct t_pbc          *pbc,
 680                 const struct t_graph        *graph,
 681                 gmx_enerdata_t              *enerd,
 682                 t_nrnb                      *nrnb,
 683                 const real                  *lambda,
 684                 const t_mdatoms             *md,
 685                 t_fcdata                    *fcd,
 686                 int                         *global_atom_index,
 687                 int                          flags)
 688 {
 689     t_pbc pbc_full; /* Full PBC is needed for position restraints */
 690
 691     if (!(flags & GMX_FORCE_LISTED))
 692     {
 693         return;
 694     }
 695
 696     if ((idef->il[F_POSRES].nr > 0) ||
 697         (idef->il[F_FBPOSRES].nr > 0))
 698     {
 699         /* Not enough flops to bother counting */
 700         set_pbc(&pbc_full, fr->ePBC, box);
 701     }
 702     calc_listed(cr, wcycle, idef, x, hist,
 703                 forceForUseWithShiftForces, forceWithVirial,
 704                 fr, pbc, &pbc_full,
 705                 graph, enerd, nrnb, lambda, md, fcd,
 706                 global_atom_index, flags);
 707
 708     /* Check if we have to determine energy differences
 709      * at foreign lambda's.
 710      */
 711     if (fepvals->n_lambda > 0 && (flags & GMX_FORCE_DHDL))
 712     {
 713         posres_wrapper_lambda(wcycle, fepvals, idef, &pbc_full, x, enerd, lambda, fr);
 714
 715         if (idef->ilsort != ilsortNO_FE)
 716         {
 717             wallcycle_sub_start(wcycle, ewcsLISTED_FEP);
 718             if (idef->ilsort != ilsortFE_SORTED)
 719             {
 720                 gmx_incons("The bonded interactions are not sorted for free energy");
 721             }
 722             for (int i = 0; i < enerd->n_lambda; i++)
 723             {
 724                 real lam_i[efptNR];
 725
 726                 reset_foreign_enerdata(enerd);
 727                 for (int j = 0; j < efptNR; j++)
 728                 {
 729                     lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i-1]);
 730                 }
 731                 calc_listed_lambda(idef, x, fr, pbc, graph, &(enerd->foreign_grpp), enerd->foreign_term, nrnb, lam_i, md,
 732                                    fcd, global_atom_index);
 733                 sum_epot(&(enerd->foreign_grpp), enerd->foreign_term);
 734                 enerd->enerpart_lambda[i] += enerd->foreign_term[F_EPOT];
 735             }
 736             wallcycle_sub_stop(wcycle, ewcsLISTED_FEP);
 737         }
 738     }
 739 }