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