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     if (!bt->haveBondeds)
 207     {
 208         return;
 209     }
 210
 211     if (bt->nblock_used > 0)
 212     {
 213         /* Reduce the bonded force buffer */
 214         reduce_thread_forces(n, f, bt, bt->nthreads);
 215     }
 216
 217     /* When necessary, reduce energy and virial using one thread only */
 218     if (bCalcEnerVir && bt->nthreads > 1)
 219     {
 220         f_thread_t *f_t = bt->f_t;
 221
 222         for (int i = 0; i < SHIFTS; i++)
 223         {
 224             for (int t = 1; t < bt->nthreads; t++)
 225             {
 226                 rvec_inc(fshift[i], f_t[t].fshift[i]);
 227             }
 228         }
 229         for (int i = 0; i < F_NRE; i++)
 230         {
 231             for (int t = 1; t < bt->nthreads; t++)
 232             {
 233                 ener[i] += f_t[t].ener[i];
 234             }
 235         }
 236         for (int i = 0; i < egNR; i++)
 237         {
 238             for (int j = 0; j < f_t[1].grpp.nener; j++)
 239             {
 240                 for (int t = 1; t < bt->nthreads; t++)
 241                 {
 242                     grpp->ener[i][j] += f_t[t].grpp.ener[i][j];
 243                 }
 244             }
 245         }
 246         if (bDHDL)
 247         {
 248             for (int i = 0; i < efptNR; i++)
 249             {
 250
 251                 for (int t = 1; t < bt->nthreads; t++)
 252                 {
 253                     dvdl[i] += f_t[t].dvdl[i];
 254                 }
 255             }
 256         }
 257     }
 258 }
 259
 260 /*! \brief Calculate one element of the list of bonded interactions
 261     for this thread */
 262 real
 263 calc_one_bond(int thread,
 264               int ftype, const t_idef *idef,
 265               const rvec x[], rvec4 f[], rvec fshift[],
 266               t_forcerec *fr,
 267               const t_pbc *pbc, const t_graph *g,
 268               gmx_grppairener_t *grpp,
 269               t_nrnb *nrnb,
 270               real *lambda, real *dvdl,
 271               const t_mdatoms *md, t_fcdata *fcd,
 272               gmx_bool bCalcEnerVir,
 273               int *global_atom_index)
 274 {
 275 #if GMX_SIMD_HAVE_REAL
 276     bool bUseSIMD = fr->use_simd_kernels;
 277 #endif
 278
 279     int      nat1, nbonds, efptFTYPE;
 280     real     v = 0;
 281     t_iatom *iatoms;
 282     int      nb0, nbn;
 283
 284     if (IS_RESTRAINT_TYPE(ftype))
 285     {
 286         efptFTYPE = efptRESTRAINT;
 287     }
 288     else
 289     {
 290         efptFTYPE = efptBONDED;
 291     }
 292
 293     nat1      = interaction_function[ftype].nratoms + 1;
 294     nbonds    = idef->il[ftype].nr/nat1;
 295     iatoms    = idef->il[ftype].iatoms;
 296
 297     nb0 = idef->il_thread_division[ftype*(idef->nthreads+1)+thread];
 298     nbn = idef->il_thread_division[ftype*(idef->nthreads+1)+thread+1] - nb0;
 299
 300     if (!isPairInteraction(ftype))
 301     {
 302         if (ftype == F_CMAP)
 303         {
 304             /* TODO The execution time for CMAP dihedrals might be
 305                nice to account to its own subtimer, but first
 306                wallcycle needs to be extended to support calling from
 307                multiple threads. */
 308             v = cmap_dihs(nbn, iatoms+nb0,
 309                           idef->iparams, &idef->cmap_grid,
 310                           x, f, fshift,
 311                           pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]),
 312                           md, fcd, global_atom_index);
 313         }
 314 #if GMX_SIMD_HAVE_REAL
 315         else if (ftype == F_ANGLES && bUseSIMD &&
 316                  !bCalcEnerVir && fr->efep == efepNO)
 317         {
 318             /* No energies, shift forces, dvdl */
 319             angles_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
 320                                idef->iparams,
 321                                x, f,
 322                                pbc, g, lambda[efptFTYPE], md, fcd,
 323                                global_atom_index);
 324             v = 0;
 325         }
 326 #endif
 327         else if (ftype == F_PDIHS &&
 328                  !bCalcEnerVir && fr->efep == efepNO)
 329         {
 330             /* No energies, shift forces, dvdl */
 331 #if GMX_SIMD_HAVE_REAL
 332             if (bUseSIMD)
 333             {
 334                 pdihs_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
 335                                   idef->iparams,
 336                                   x, f,
 337                                   pbc, g, lambda[efptFTYPE], md, fcd,
 338                                   global_atom_index);
 339             }
 340             else
 341 #endif
 342             {
 343                 pdihs_noener(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             v = 0;
 350         }
 351 #if GMX_SIMD_HAVE_REAL
 352         else if (ftype == F_RBDIHS && bUseSIMD &&
 353                  !bCalcEnerVir && fr->efep == efepNO)
 354         {
 355             /* No energies, shift forces, dvdl */
 356             rbdihs_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
 357                                idef->iparams,
 358                                (const rvec*)x, f,
 359                                pbc, g, lambda[efptFTYPE], md, fcd,
 360                                global_atom_index);
 361             v = 0;
 362         }
 363 #endif
 364         else
 365         {
 366             v = interaction_function[ftype].ifunc(nbn, iatoms+nb0,
 367                                                   idef->iparams,
 368                                                   x, f, fshift,
 369                                                   pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]),
 370                                                   md, fcd, global_atom_index);
 371         }
 372     }
 373     else
 374     {
 375         /* TODO The execution time for pairs might be nice to account
 376            to its own subtimer, but first wallcycle needs to be
 377            extended to support calling from multiple threads. */
 378         do_pairs(ftype, nbn, iatoms+nb0, idef->iparams, x, f, fshift,
 379                  pbc, g, lambda, dvdl, md, fr,
 380                  bCalcEnerVir, grpp, global_atom_index);
 381         v = 0;
 382     }
 383
 384     if (thread == 0)
 385     {
 386         inc_nrnb(nrnb, interaction_function[ftype].nrnb_ind, nbonds);
 387     }
 388
 389     return v;
 390 }
 391
 392 } // namespace
 393
 394 gmx_bool
 395 ftype_is_bonded_potential(int ftype)
 396 {
 397     return
 398         (interaction_function[ftype].flags & IF_BOND) &&
 399         !(ftype == F_CONNBONDS || ftype == F_POSRES || ftype == F_FBPOSRES) &&
 400         (ftype < F_GB12 || ftype > F_GB14);
 401 }
 402
 403 void calc_listed(const t_commrec             *cr,
 404                  struct gmx_wallcycle        *wcycle,
 405                  const t_idef *idef,
 406                  const rvec x[], history_t *hist,
 407                  rvec f[], t_forcerec *fr,
 408                  const struct t_pbc *pbc,
 409                  const struct t_pbc *pbc_full,
 410                  const struct t_graph *g,
 411                  gmx_enerdata_t *enerd, t_nrnb *nrnb,
 412                  real *lambda,
 413                  const t_mdatoms *md,
 414                  t_fcdata *fcd, int *global_atom_index,
 415                  int force_flags)
 416 {
 417     struct bonded_threading_t *bt;
 418     gmx_bool                   bCalcEnerVir;
 419     int                        i;
 420     /* The dummy array is to have a place to store the dhdl at other values
 421        of lambda, which will be thrown away in the end */
 422     real                       dvdl[efptNR];
 423     const  t_pbc              *pbc_null;
 424     int                        thread;
 425
 426     bt = fr->bonded_threading;
 427
 428     assert(bt->nthreads == idef->nthreads);
 429
 430     bCalcEnerVir = (force_flags & (GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY));
 431
 432     for (i = 0; i < efptNR; i++)
 433     {
 434         dvdl[i] = 0.0;
 435     }
 436     if (fr->bMolPBC)
 437     {
 438         pbc_null = pbc;
 439     }
 440     else
 441     {
 442         pbc_null = nullptr;
 443     }
 444
 445 #ifdef DEBUG
 446     if (g && debug)
 447     {
 448         p_graph(debug, "Bondage is fun", g);
 449     }
 450 #endif
 451
 452     if ((idef->il[F_POSRES].nr > 0) ||
 453         (idef->il[F_FBPOSRES].nr > 0) ||
 454         fcd->orires.nr > 0 ||
 455         fcd->disres.nres > 0)
 456     {
 457         /* TODO Use of restraints triggers further function calls
 458            inside the loop over calc_one_bond(), but those are too
 459            awkward to account to this subtimer properly in the present
 460            code. We don't test / care much about performance with
 461            restraints, anyway. */
 462         wallcycle_sub_start(wcycle, ewcsRESTRAINTS);
 463
 464         if (idef->il[F_POSRES].nr > 0)
 465         {
 466             posres_wrapper(nrnb, idef, pbc_full, x, enerd, lambda, fr);
 467         }
 468
 469         if (idef->il[F_FBPOSRES].nr > 0)
 470         {
 471             fbposres_wrapper(nrnb, idef, pbc_full, x, enerd, fr);
 472         }
 473
 474         /* Do pre force calculation stuff which might require communication */
 475         if (fcd->orires.nr > 0)
 476         {
 477             enerd->term[F_ORIRESDEV] =
 478                 calc_orires_dev(cr->ms, idef->il[F_ORIRES].nr,
 479                                 idef->il[F_ORIRES].iatoms,
 480                                 idef->iparams, md, x,
 481                                 pbc_null, fcd, hist);
 482         }
 483         if (fcd->disres.nres > 0)
 484         {
 485             calc_disres_R_6(cr,
 486                             idef->il[F_DISRES].nr,
 487                             idef->il[F_DISRES].iatoms,
 488                             x, pbc_null,
 489                             fcd, hist);
 490         }
 491
 492         wallcycle_sub_stop(wcycle, ewcsRESTRAINTS);
 493     }
 494
 495     /* TODO: Skip this whole loop with a system/domain without listeds */
 496     wallcycle_sub_start(wcycle, ewcsLISTED);
 497 #pragma omp parallel for num_threads(bt->nthreads) schedule(static)
 498     for (thread = 0; thread < bt->nthreads; thread++)
 499     {
 500         try
 501         {
 502             int                ftype;
 503             real              *epot, v;
 504             /* thread stuff */
 505             rvec4             *ft;
 506             rvec              *fshift;
 507             real              *dvdlt;
 508             gmx_grppairener_t *grpp;
 509
 510             zero_thread_output(bt, thread);
 511
 512             ft = bt->f_t[thread].f;
 513
 514             if (thread == 0)
 515             {
 516                 fshift = fr->fshift;
 517                 epot   = enerd->term;
 518                 grpp   = &enerd->grpp;
 519                 dvdlt  = dvdl;
 520             }
 521             else
 522             {
 523                 fshift = bt->f_t[thread].fshift;
 524                 epot   = bt->f_t[thread].ener;
 525                 grpp   = &bt->f_t[thread].grpp;
 526                 dvdlt  = bt->f_t[thread].dvdl;
 527             }
 528             /* Loop over all bonded force types to calculate the bonded forces */
 529             for (ftype = 0; (ftype < F_NRE); ftype++)
 530             {
 531                 if (idef->il[ftype].nr > 0 && ftype_is_bonded_potential(ftype))
 532                 {
 533                     v = calc_one_bond(thread, ftype, idef, x,
 534                                       ft, fshift, fr, pbc_null, g, grpp,
 535                                       nrnb, lambda, dvdlt,
 536                                       md, fcd, bCalcEnerVir,
 537                                       global_atom_index);
 538                     epot[ftype] += v;
 539                 }
 540             }
 541         }
 542         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 543     }
 544     wallcycle_sub_stop(wcycle, ewcsLISTED);
 545
 546     wallcycle_sub_start(wcycle, ewcsLISTED_BUF_OPS);
 547     reduce_thread_output(fr->natoms_force, f, fr->fshift,
 548                          enerd->term, &enerd->grpp, dvdl,
 549                          bt,
 550                          bCalcEnerVir,
 551                          force_flags & GMX_FORCE_DHDL);
 552     wallcycle_sub_stop(wcycle, ewcsLISTED_BUF_OPS);
 553
 554     /* Remaining code does not have enough flops to bother counting */
 555     if (force_flags & GMX_FORCE_DHDL)
 556     {
 557         for (i = 0; i < efptNR; i++)
 558         {
 559             enerd->dvdl_nonlin[i] += dvdl[i];
 560         }
 561     }
 562
 563     /* Copy the sum of violations for the distance restraints from fcd */
 564     if (fcd)
 565     {
 566         enerd->term[F_DISRESVIOL] = fcd->disres.sumviol;
 567
 568     }
 569 }
 570
 571 void calc_listed_lambda(const t_idef *idef,
 572                         const rvec x[],
 573                         t_forcerec *fr,
 574                         const struct t_pbc *pbc, const struct t_graph *g,
 575                         gmx_grppairener_t *grpp, real *epot, t_nrnb *nrnb,
 576                         real *lambda,
 577                         const t_mdatoms *md,
 578                         t_fcdata *fcd,
 579                         int *global_atom_index)
 580 {
 581     int           ftype, nr_nonperturbed, nr;
 582     real          v;
 583     real          dvdl_dum[efptNR] = {0};
 584     rvec4        *f;
 585     rvec         *fshift;
 586     const  t_pbc *pbc_null;
 587     t_idef        idef_fe;
 588
 589     if (fr->bMolPBC)
 590     {
 591         pbc_null = pbc;
 592     }
 593     else
 594     {
 595         pbc_null = nullptr;
 596     }
 597
 598     /* Copy the whole idef, so we can modify the contents locally */
 599     idef_fe          = *idef;
 600     idef_fe.nthreads = 1;
 601     snew(idef_fe.il_thread_division, F_NRE*(idef_fe.nthreads+1));
 602
 603     /* We already have the forces, so we use temp buffers here */
 604     snew(f, fr->natoms_force);
 605     snew(fshift, SHIFTS);
 606
 607     /* Loop over all bonded force types to calculate the bonded energies */
 608     for (ftype = 0; (ftype < F_NRE); ftype++)
 609     {
 610         if (ftype_is_bonded_potential(ftype))
 611         {
 612             /* Set the work range of thread 0 to the perturbed bondeds only */
 613             nr_nonperturbed                       = idef->il[ftype].nr_nonperturbed;
 614             nr                                    = idef->il[ftype].nr;
 615             idef_fe.il_thread_division[ftype*2+0] = nr_nonperturbed;
 616             idef_fe.il_thread_division[ftype*2+1] = nr;
 617
 618             /* This is only to get the flop count correct */
 619             idef_fe.il[ftype].nr = nr - nr_nonperturbed;
 620
 621             if (nr - nr_nonperturbed > 0)
 622             {
 623                 v = calc_one_bond(0, ftype, &idef_fe,
 624                                   x, f, fshift, fr, pbc_null, g,
 625                                   grpp, nrnb, lambda, dvdl_dum,
 626                                   md, fcd, TRUE,
 627                                   global_atom_index);
 628                 epot[ftype] += v;
 629             }
 630         }
 631     }
 632
 633     sfree(fshift);
 634     sfree(f);
 635
 636     sfree(idef_fe.il_thread_division);
 637 }
 638
 639 void
 640 do_force_listed(struct gmx_wallcycle        *wcycle,
 641                 matrix                       box,
 642                 const t_lambda              *fepvals,
 643                 const t_commrec             *cr,
 644                 const t_idef                *idef,
 645                 const rvec                   x[],
 646                 history_t                   *hist,
 647                 rvec                         f[],
 648                 t_forcerec                  *fr,
 649                 const struct t_pbc          *pbc,
 650                 const struct t_graph        *graph,
 651                 gmx_enerdata_t              *enerd,
 652                 t_nrnb                      *nrnb,
 653                 real                        *lambda,
 654                 const t_mdatoms             *md,
 655                 t_fcdata                    *fcd,
 656                 int                         *global_atom_index,
 657                 int                          flags)
 658 {
 659     t_pbc pbc_full; /* Full PBC is needed for position restraints */
 660
 661     if (!(flags & GMX_FORCE_LISTED))
 662     {
 663         return;
 664     }
 665
 666     if ((idef->il[F_POSRES].nr > 0) ||
 667         (idef->il[F_FBPOSRES].nr > 0))
 668     {
 669         /* Not enough flops to bother counting */
 670         set_pbc(&pbc_full, fr->ePBC, box);
 671     }
 672     calc_listed(cr, wcycle, idef, x, hist, f, fr, pbc, &pbc_full,
 673                 graph, enerd, nrnb, lambda, md, fcd,
 674                 global_atom_index, flags);
 675
 676     /* Check if we have to determine energy differences
 677      * at foreign lambda's.
 678      */
 679     if (fepvals->n_lambda > 0 && (flags & GMX_FORCE_DHDL))
 680     {
 681         posres_wrapper_lambda(wcycle, fepvals, idef, &pbc_full, x, enerd, lambda, fr);
 682
 683         if (idef->ilsort != ilsortNO_FE)
 684         {
 685             wallcycle_sub_start(wcycle, ewcsLISTED_FEP);
 686             if (idef->ilsort != ilsortFE_SORTED)
 687             {
 688                 gmx_incons("The bonded interactions are not sorted for free energy");
 689             }
 690             for (int i = 0; i < enerd->n_lambda; i++)
 691             {
 692                 real lam_i[efptNR];
 693
 694                 reset_foreign_enerdata(enerd);
 695                 for (int j = 0; j < efptNR; j++)
 696                 {
 697                     lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i-1]);
 698                 }
 699                 calc_listed_lambda(idef, x, fr, pbc, graph, &(enerd->foreign_grpp), enerd->foreign_term, nrnb, lam_i, md,
 700                                    fcd, global_atom_index);
 701                 sum_epot(&(enerd->foreign_grpp), enerd->foreign_term);
 702                 enerd->enerpart_lambda[i] += enerd->foreign_term[F_EPOT];
 703             }
 704             wallcycle_sub_stop(wcycle, ewcsLISTED_FEP);
 705         }
 706     }
 707 }