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