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, 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 "config.h"
  49
  50 #include <assert.h>
  51
  52 #include <algorithm>
  53
  54 #include "gromacs/legacyheaders/disre.h"
  55 #include "gromacs/legacyheaders/force.h"
  56 #include "gromacs/legacyheaders/network.h"
  57 #include "gromacs/legacyheaders/nrnb.h"
  58 #include "gromacs/legacyheaders/orires.h"
  59 #include "gromacs/legacyheaders/types/force_flags.h"
  60 #include "gromacs/listed-forces/bonded.h"
  61 #include "gromacs/listed-forces/position-restraints.h"
  62 #include "gromacs/math/vec.h"
  63 #include "gromacs/mdlib/forcerec-threading.h"
  64 #include "gromacs/pbcutil/ishift.h"
  65 #include "gromacs/pbcutil/pbc.h"
  66 #include "gromacs/simd/simd.h"
  67 #include "gromacs/timing/wallcycle.h"
  68 #include "gromacs/utility/smalloc.h"
  69
  70 #include "pairs.h"
  71
  72 namespace
  73 {
  74
  75 /*! \brief Return true if ftype is an explicit pair-listed LJ or
  76  * COULOMB interaction type: bonded LJ (usually 1-4), or special
  77  * listed non-bonded for FEP. */
  78 bool
  79 isPairInteraction(int ftype)
  80 {
  81     return ((ftype) >= F_LJ14 && (ftype) <= F_LJC_PAIRS_NB);
  82 }
  83
  84 /*! \brief Zero thread-local force-output buffers */
  85 void
  86 zero_thread_forces(f_thread_t *f_t, int n,
  87                    int nblock, int blocksize)
  88 {
  89     int b, a0, a1, a, i, j;
  90
  91     if (n > f_t->f_nalloc)
  92     {
  93         f_t->f_nalloc = over_alloc_large(n);
  94         srenew(f_t->f, f_t->f_nalloc);
  95     }
  96
  97     if (!bitmask_is_zero(f_t->red_mask))
  98     {
  99         for (b = 0; b < nblock; b++)
 100         {
 101             if (bitmask_is_set(f_t->red_mask, b))
 102             {
 103                 a0 = b*blocksize;
 104                 a1 = std::min((b+1)*blocksize, n);
 105                 for (a = a0; a < a1; a++)
 106                 {
 107                     clear_rvec(f_t->f[a]);
 108                 }
 109             }
 110         }
 111     }
 112     for (i = 0; i < SHIFTS; i++)
 113     {
 114         clear_rvec(f_t->fshift[i]);
 115     }
 116     for (i = 0; i < F_NRE; i++)
 117     {
 118         f_t->ener[i] = 0;
 119     }
 120     for (i = 0; i < egNR; i++)
 121     {
 122         for (j = 0; j < f_t->grpp.nener; j++)
 123         {
 124             f_t->grpp.ener[i][j] = 0;
 125         }
 126     }
 127     for (i = 0; i < efptNR; i++)
 128     {
 129         f_t->dvdl[i] = 0;
 130     }
 131 }
 132
 133 /*! \brief The max thread number is arbitrary, we used a fixed number
 134  * to avoid memory management.  Using more than 16 threads is probably
 135  * never useful performance wise. */
 136 #define MAX_BONDED_THREADS 256
 137
 138 /*! \brief Reduce thread-local force buffers */
 139 void
 140 reduce_thread_force_buffer(int n, rvec *f,
 141                            int nthreads, f_thread_t *f_t,
 142                            int nblock, int block_size)
 143 {
 144     int b;
 145
 146     if (nthreads > MAX_BONDED_THREADS)
 147     {
 148         gmx_fatal(FARGS, "Can not reduce bonded forces on more than %d threads",
 149                   MAX_BONDED_THREADS);
 150     }
 151
 152     /* This reduction can run on any number of threads,
 153      * independently of nthreads.
 154      */
 155 #pragma omp parallel for num_threads(nthreads) schedule(static)
 156     for (b = 0; b < nblock; b++)
 157     {
 158         rvec *fp[MAX_BONDED_THREADS];
 159         int   nfb, ft, fb;
 160         int   a0, a1, a;
 161
 162         /* Determine which threads contribute to this block */
 163         nfb = 0;
 164         for (ft = 1; ft < nthreads; ft++)
 165         {
 166             if (bitmask_is_set(f_t[ft].red_mask, b))
 167             {
 168                 fp[nfb++] = f_t[ft].f;
 169             }
 170         }
 171         if (nfb > 0)
 172         {
 173             /* Reduce force buffers for threads that contribute */
 174             a0 =  b   *block_size;
 175             a1 = (b+1)*block_size;
 176             a1 = std::min(a1, n);
 177             for (a = a0; a < a1; a++)
 178             {
 179                 for (fb = 0; fb < nfb; fb++)
 180                 {
 181                     rvec_inc(f[a], fp[fb][a]);
 182                 }
 183             }
 184         }
 185     }
 186 }
 187
 188 /*! \brief Reduce thread-local forces */
 189 void
 190 reduce_thread_forces(int n, rvec *f, rvec *fshift,
 191                      real *ener, gmx_grppairener_t *grpp, real *dvdl,
 192                      int nthreads, f_thread_t *f_t,
 193                      int nblock, int block_size,
 194                      gmx_bool bCalcEnerVir,
 195                      gmx_bool bDHDL)
 196 {
 197     if (nblock > 0)
 198     {
 199         /* Reduce the bonded force buffer */
 200         reduce_thread_force_buffer(n, f, nthreads, f_t, nblock, block_size);
 201     }
 202
 203     /* When necessary, reduce energy and virial using one thread only */
 204     if (bCalcEnerVir)
 205     {
 206         int t, i, j;
 207
 208         for (i = 0; i < SHIFTS; i++)
 209         {
 210             for (t = 1; t < nthreads; t++)
 211             {
 212                 rvec_inc(fshift[i], f_t[t].fshift[i]);
 213             }
 214         }
 215         for (i = 0; i < F_NRE; i++)
 216         {
 217             for (t = 1; t < nthreads; t++)
 218             {
 219                 ener[i] += f_t[t].ener[i];
 220             }
 221         }
 222         for (i = 0; i < egNR; i++)
 223         {
 224             for (j = 0; j < f_t[1].grpp.nener; j++)
 225             {
 226                 for (t = 1; t < nthreads; t++)
 227                 {
 228
 229                     grpp->ener[i][j] += f_t[t].grpp.ener[i][j];
 230                 }
 231             }
 232         }
 233         if (bDHDL)
 234         {
 235             for (i = 0; i < efptNR; i++)
 236             {
 237
 238                 for (t = 1; t < nthreads; t++)
 239                 {
 240                     dvdl[i] += f_t[t].dvdl[i];
 241                 }
 242             }
 243         }
 244     }
 245 }
 246
 247 /*! \brief Calculate one element of the list of bonded interactions
 248     for this thread */
 249 real
 250 calc_one_bond(int thread,
 251               int ftype, const t_idef *idef,
 252               const rvec x[], rvec f[], rvec fshift[],
 253               t_forcerec *fr,
 254               const t_pbc *pbc, const t_graph *g,
 255               gmx_grppairener_t *grpp,
 256               t_nrnb *nrnb,
 257               real *lambda, real *dvdl,
 258               const t_mdatoms *md, t_fcdata *fcd,
 259               gmx_bool bCalcEnerVir,
 260               int *global_atom_index)
 261 {
 262 #ifdef GMX_SIMD_HAVE_REAL
 263     gmx_bool bUseSIMD;
 264     /* MSVC 2010 produces buggy SIMD PBC code, disable SIMD for MSVC <= 2010 */
 265 #if defined _MSC_VER && _MSC_VER < 1700
 266     bUseSIMD = FALSE;
 267 #else
 268     bUseSIMD = fr->use_simd_kernels;
 269 #endif
 270 #endif
 271
 272     int      nat1, nbonds, efptFTYPE;
 273     real     v = 0;
 274     t_iatom *iatoms;
 275     int      nb0, nbn;
 276
 277     if (IS_RESTRAINT_TYPE(ftype))
 278     {
 279         efptFTYPE = efptRESTRAINT;
 280     }
 281     else
 282     {
 283         efptFTYPE = efptBONDED;
 284     }
 285
 286     nat1      = interaction_function[ftype].nratoms + 1;
 287     nbonds    = idef->il[ftype].nr/nat1;
 288     iatoms    = idef->il[ftype].iatoms;
 289
 290     nb0 = idef->il_thread_division[ftype*(idef->nthreads+1)+thread];
 291     nbn = idef->il_thread_division[ftype*(idef->nthreads+1)+thread+1] - nb0;
 292
 293     if (!isPairInteraction(ftype))
 294     {
 295         if (ftype == F_CMAP)
 296         {
 297             /* TODO The execution time for CMAP dihedrals might be
 298                nice to account to its own subtimer, but first
 299                wallcycle needs to be extended to support calling from
 300                multiple threads. */
 301             v = cmap_dihs(nbn, iatoms+nb0,
 302                           idef->iparams, &idef->cmap_grid,
 303                           x, f, fshift,
 304                           pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]),
 305                           md, fcd, global_atom_index);
 306         }
 307 #ifdef GMX_SIMD_HAVE_REAL
 308         else if (ftype == F_ANGLES && bUseSIMD &&
 309                  !bCalcEnerVir && fr->efep == efepNO)
 310         {
 311             /* No energies, shift forces, dvdl */
 312             angles_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
 313                                idef->iparams,
 314                                x, f,
 315                                pbc, g, lambda[efptFTYPE], md, fcd,
 316                                global_atom_index);
 317             v = 0;
 318         }
 319 #endif
 320         else if (ftype == F_PDIHS &&
 321                  !bCalcEnerVir && fr->efep == efepNO)
 322         {
 323             /* No energies, shift forces, dvdl */
 324 #ifdef GMX_SIMD_HAVE_REAL
 325             if (bUseSIMD)
 326             {
 327                 pdihs_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
 328                                   idef->iparams,
 329                                   x, f,
 330                                   pbc, g, lambda[efptFTYPE], md, fcd,
 331                                   global_atom_index);
 332             }
 333             else
 334 #endif
 335             {
 336                 pdihs_noener(nbn, idef->il[ftype].iatoms+nb0,
 337                              idef->iparams,
 338                              x, f,
 339                              pbc, g, lambda[efptFTYPE], md, fcd,
 340                              global_atom_index);
 341             }
 342             v = 0;
 343         }
 344 #ifdef GMX_SIMD_HAVE_REAL
 345         else if (ftype == F_RBDIHS && bUseSIMD &&
 346                  !bCalcEnerVir && fr->efep == efepNO)
 347         {
 348             /* No energies, shift forces, dvdl */
 349             rbdihs_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
 350                                idef->iparams,
 351                                (const rvec*)x, f,
 352                                pbc, g, lambda[efptFTYPE], md, fcd,
 353                                global_atom_index);
 354             v = 0;
 355         }
 356 #endif
 357         else
 358         {
 359             v = interaction_function[ftype].ifunc(nbn, iatoms+nb0,
 360                                                   idef->iparams,
 361                                                   x, f, fshift,
 362                                                   pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]),
 363                                                   md, fcd, global_atom_index);
 364         }
 365     }
 366     else
 367     {
 368         /* TODO The execution time for pairs might be nice to account
 369            to its own subtimer, but first wallcycle needs to be
 370            extended to support calling from multiple threads. */
 371         v = do_pairs(ftype, nbn, iatoms+nb0, idef->iparams, x, f, fshift,
 372                      pbc, g, lambda, dvdl, md, fr, grpp, global_atom_index);
 373     }
 374
 375     if (thread == 0)
 376     {
 377         inc_nrnb(nrnb, interaction_function[ftype].nrnb_ind, nbonds);
 378     }
 379
 380     return v;
 381 }
 382
 383 } // namespace
 384
 385 gmx_bool
 386 ftype_is_bonded_potential(int ftype)
 387 {
 388     return
 389         (interaction_function[ftype].flags & IF_BOND) &&
 390         !(ftype == F_CONNBONDS || ftype == F_POSRES || ftype == F_FBPOSRES) &&
 391         (ftype < F_GB12 || ftype > F_GB14);
 392 }
 393
 394 void calc_listed(const gmx_multisim_t *ms,
 395                  gmx_wallcycle        *wcycle,
 396                  const t_idef *idef,
 397                  const rvec x[], history_t *hist,
 398                  rvec f[], t_forcerec *fr,
 399                  const struct t_pbc *pbc,
 400                  const struct t_pbc *pbc_full,
 401                  const struct t_graph *g,
 402                  gmx_enerdata_t *enerd, t_nrnb *nrnb,
 403                  real *lambda,
 404                  const t_mdatoms *md,
 405                  t_fcdata *fcd, int *global_atom_index,
 406                  int force_flags)
 407 {
 408     gmx_bool      bCalcEnerVir;
 409     int           i;
 410     real          dvdl[efptNR]; /* The dummy array is to have a place to store the dhdl at other values
 411                                                         of lambda, which will be thrown away in the end*/
 412     const  t_pbc *pbc_null;
 413     int           thread;
 414
 415     assert(fr->nthreads == idef->nthreads);
 416
 417     bCalcEnerVir = (force_flags & (GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY));
 418
 419     for (i = 0; i < efptNR; i++)
 420     {
 421         dvdl[i] = 0.0;
 422     }
 423     if (fr->bMolPBC)
 424     {
 425         pbc_null = pbc;
 426     }
 427     else
 428     {
 429         pbc_null = NULL;
 430     }
 431
 432 #ifdef DEBUG
 433     if (g && debug)
 434     {
 435         p_graph(debug, "Bondage is fun", g);
 436     }
 437 #endif
 438
 439     if ((idef->il[F_POSRES].nr > 0) ||
 440         (idef->il[F_FBPOSRES].nr > 0) ||
 441         (idef->il[F_ORIRES].nr > 0) ||
 442         (idef->il[F_DISRES].nr > 0))
 443     {
 444         /* TODO Use of restraints triggers further function calls
 445            inside the loop over calc_one_bond(), but those are too
 446            awkward to account to this subtimer properly in the present
 447            code. We don't test / care much about performance with
 448            restraints, anyway. */
 449         wallcycle_sub_start(wcycle, ewcsRESTRAINTS);
 450
 451         if (idef->il[F_POSRES].nr > 0)
 452         {
 453             posres_wrapper(nrnb, idef, pbc_full, x, enerd, lambda, fr);
 454         }
 455
 456         if (idef->il[F_FBPOSRES].nr > 0)
 457         {
 458             fbposres_wrapper(nrnb, idef, pbc_full, x, enerd, fr);
 459         }
 460
 461         /* Do pre force calculation stuff which might require communication */
 462         if (idef->il[F_ORIRES].nr > 0)
 463         {
 464             enerd->term[F_ORIRESDEV] =
 465                 calc_orires_dev(ms, idef->il[F_ORIRES].nr,
 466                                 idef->il[F_ORIRES].iatoms,
 467                                 idef->iparams, md, x,
 468                                 pbc_null, fcd, hist);
 469         }
 470         if (idef->il[F_DISRES].nr)
 471         {
 472             calc_disres_R_6(idef->il[F_DISRES].nr,
 473                             idef->il[F_DISRES].iatoms,
 474                             idef->iparams, x, pbc_null,
 475                             fcd, hist);
 476 #ifdef GMX_MPI
 477             if (fcd->disres.nsystems > 1)
 478             {
 479                 gmx_sum_sim(2*fcd->disres.nres, fcd->disres.Rt_6, ms);
 480             }
 481 #endif
 482         }
 483
 484         wallcycle_sub_stop(wcycle, ewcsRESTRAINTS);
 485     }
 486
 487     wallcycle_sub_start(wcycle, ewcsLISTED);
 488 #pragma omp parallel for num_threads(fr->nthreads) schedule(static)
 489     for (thread = 0; thread < fr->nthreads; thread++)
 490     {
 491         int                ftype;
 492         real              *epot, v;
 493         /* thread stuff */
 494         rvec              *ft, *fshift;
 495         real              *dvdlt;
 496         gmx_grppairener_t *grpp;
 497
 498         if (thread == 0)
 499         {
 500             ft     = f;
 501             fshift = fr->fshift;
 502             epot   = enerd->term;
 503             grpp   = &enerd->grpp;
 504             dvdlt  = dvdl;
 505         }
 506         else
 507         {
 508             zero_thread_forces(&fr->f_t[thread], fr->natoms_force,
 509                                fr->red_nblock, 1<<fr->red_ashift);
 510
 511             ft     = fr->f_t[thread].f;
 512             fshift = fr->f_t[thread].fshift;
 513             epot   = fr->f_t[thread].ener;
 514             grpp   = &fr->f_t[thread].grpp;
 515             dvdlt  = fr->f_t[thread].dvdl;
 516         }
 517         /* Loop over all bonded force types to calculate the bonded forces */
 518         for (ftype = 0; (ftype < F_NRE); ftype++)
 519         {
 520             if (idef->il[ftype].nr > 0 && ftype_is_bonded_potential(ftype))
 521             {
 522                 v = calc_one_bond(thread, ftype, idef, x,
 523                                   ft, fshift, fr, pbc_null, g, grpp,
 524                                   nrnb, lambda, dvdlt,
 525                                   md, fcd, bCalcEnerVir,
 526                                   global_atom_index);
 527                 epot[ftype] += v;
 528             }
 529         }
 530     }
 531     wallcycle_sub_stop(wcycle, ewcsLISTED);
 532
 533     if (fr->nthreads > 1)
 534     {
 535         wallcycle_sub_start(wcycle, ewcsLISTED_BUF_OPS);
 536         reduce_thread_forces(fr->natoms_force, f, fr->fshift,
 537                              enerd->term, &enerd->grpp, dvdl,
 538                              fr->nthreads, fr->f_t,
 539                              fr->red_nblock, 1<<fr->red_ashift,
 540                              bCalcEnerVir,
 541                              force_flags & GMX_FORCE_DHDL);
 542         wallcycle_sub_stop(wcycle, ewcsLISTED_BUF_OPS);
 543     }
 544
 545     /* Remaining code does not have enough flops to bother counting */
 546     if (force_flags & GMX_FORCE_DHDL)
 547     {
 548         for (i = 0; i < efptNR; i++)
 549         {
 550             enerd->dvdl_nonlin[i] += dvdl[i];
 551         }
 552     }
 553
 554     /* Copy the sum of violations for the distance restraints from fcd */
 555     if (fcd)
 556     {
 557         enerd->term[F_DISRESVIOL] = fcd->disres.sumviol;
 558
 559     }
 560 }
 561
 562 void calc_listed_lambda(const t_idef *idef,
 563                         const rvec x[],
 564                         t_forcerec *fr,
 565                         const struct t_pbc *pbc, const struct t_graph *g,
 566                         gmx_grppairener_t *grpp, real *epot, t_nrnb *nrnb,
 567                         real *lambda,
 568                         const t_mdatoms *md,
 569                         t_fcdata *fcd,
 570                         int *global_atom_index)
 571 {
 572     int           ftype, nr_nonperturbed, nr;
 573     real          v;
 574     real          dvdl_dum[efptNR] = {0};
 575     rvec         *f, *fshift;
 576     const  t_pbc *pbc_null;
 577     t_idef        idef_fe;
 578
 579     if (fr->bMolPBC)
 580     {
 581         pbc_null = pbc;
 582     }
 583     else
 584     {
 585         pbc_null = NULL;
 586     }
 587
 588     /* Copy the whole idef, so we can modify the contents locally */
 589     idef_fe          = *idef;
 590     idef_fe.nthreads = 1;
 591     snew(idef_fe.il_thread_division, F_NRE*(idef_fe.nthreads+1));
 592
 593     /* We already have the forces, so we use temp buffers here */
 594     snew(f, fr->natoms_force);
 595     snew(fshift, SHIFTS);
 596
 597     /* Loop over all bonded force types to calculate the bonded energies */
 598     for (ftype = 0; (ftype < F_NRE); ftype++)
 599     {
 600         if (ftype_is_bonded_potential(ftype))
 601         {
 602             /* Set the work range of thread 0 to the perturbed bondeds only */
 603             nr_nonperturbed                       = idef->il[ftype].nr_nonperturbed;
 604             nr                                    = idef->il[ftype].nr;
 605             idef_fe.il_thread_division[ftype*2+0] = nr_nonperturbed;
 606             idef_fe.il_thread_division[ftype*2+1] = nr;
 607
 608             /* This is only to get the flop count correct */
 609             idef_fe.il[ftype].nr = nr - nr_nonperturbed;
 610
 611             if (nr - nr_nonperturbed > 0)
 612             {
 613                 v = calc_one_bond(0, ftype, &idef_fe,
 614                                   x, f, fshift, fr, pbc_null, g,
 615                                   grpp, nrnb, lambda, dvdl_dum,
 616                                   md, fcd, TRUE,
 617                                   global_atom_index);
 618                 epot[ftype] += v;
 619             }
 620         }
 621     }
 622
 623     sfree(fshift);
 624     sfree(f);
 625
 626     sfree(idef_fe.il_thread_division);
 627 }
 628
 629 void
 630 do_force_listed(gmx_wallcycle        *wcycle,
 631                 matrix                box,
 632                 const t_lambda       *fepvals,
 633                 const gmx_multisim_t *ms,
 634                 const t_idef         *idef,
 635                 const rvec            x[],
 636                 history_t            *hist,
 637                 rvec                  f[],
 638                 t_forcerec           *fr,
 639                 const struct t_pbc   *pbc,
 640                 const struct t_graph *graph,
 641                 gmx_enerdata_t       *enerd,
 642                 t_nrnb               *nrnb,
 643                 real                 *lambda,
 644                 const t_mdatoms      *md,
 645                 t_fcdata             *fcd,
 646                 int                  *global_atom_index,
 647                 int                   flags)
 648 {
 649     t_pbc pbc_full; /* Full PBC is needed for position restraints */
 650
 651     if (!(flags & GMX_FORCE_LISTED))
 652     {
 653         return;
 654     }
 655
 656     if ((idef->il[F_POSRES].nr > 0) ||
 657         (idef->il[F_FBPOSRES].nr > 0))
 658     {
 659         /* Not enough flops to bother counting */
 660         set_pbc(&pbc_full, fr->ePBC, box);
 661     }
 662     calc_listed(ms, wcycle, idef, x, hist, f, fr, pbc, &pbc_full,
 663                 graph, enerd, nrnb, lambda, md, fcd,
 664                 global_atom_index, flags);
 665
 666     /* Check if we have to determine energy differences
 667      * at foreign lambda's.
 668      */
 669     if (fepvals->n_lambda > 0 && (flags & GMX_FORCE_DHDL))
 670     {
 671         posres_wrapper_lambda(wcycle, fepvals, idef, &pbc_full, x, enerd, lambda, fr);
 672
 673         if (idef->ilsort != ilsortNO_FE)
 674         {
 675             wallcycle_sub_start(wcycle, ewcsLISTED_FEP);
 676             if (idef->ilsort != ilsortFE_SORTED)
 677             {
 678                 gmx_incons("The bonded interactions are not sorted for free energy");
 679             }
 680             for (int i = 0; i < enerd->n_lambda; i++)
 681             {
 682                 real lam_i[efptNR];
 683
 684                 reset_foreign_enerdata(enerd);
 685                 for (int j = 0; j < efptNR; j++)
 686                 {
 687                     lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i-1]);
 688                 }
 689                 calc_listed_lambda(idef, x, fr, pbc, graph, &(enerd->foreign_grpp), enerd->foreign_term, nrnb, lam_i, md,
 690                                    fcd, global_atom_index);
 691                 sum_epot(&(enerd->foreign_grpp), enerd->foreign_term);
 692                 enerd->enerpart_lambda[i] += enerd->foreign_term[F_EPOT];
 693             }
 694             wallcycle_sub_stop(wcycle, ewcsLISTED_FEP);
 695         }
 696     }
 697     debug_gmx();
 698 }