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/force.h"
  62 #include "gromacs/mdlib/force_flags.h"
  63 #include "gromacs/mdtypes/commrec.h"
  64 #include "gromacs/mdtypes/fcdata.h"
  65 #include "gromacs/mdtypes/forcerec.h"
  66 #include "gromacs/mdtypes/inputrec.h"
  67 #include "gromacs/mdtypes/md_enums.h"
  68 #include "gromacs/pbcutil/ishift.h"
  69 #include "gromacs/pbcutil/pbc.h"
  70 #include "gromacs/simd/simd.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 struct BondedInteractions
  81 {
  82     BondedFunction function;
  83     int            nrnbIndex;
  84 };
  85
  86 /*! \brief Lookup table of bonded interaction functions
  87  *
  88  * This must have as many entries as interaction_function in ifunc.cpp */
  89 static std::array<BondedInteractions, F_NRE> s_bondedInteractionFunctions
  90     = {
  91     BondedInteractions {bonds, eNR_BONDS },                       // F_BONDS
  92     BondedInteractions {g96bonds, eNR_BONDS },                    // F_G96BONDS
  93     BondedInteractions {morse_bonds, eNR_MORSE },                 // F_MORSE
  94     BondedInteractions {cubic_bonds, eNR_CUBICBONDS },            // F_CUBICBONDS
  95     BondedInteractions {unimplemented, -1 },                      // F_CONNBONDS
  96     BondedInteractions {bonds, eNR_BONDS },                       // F_HARMONIC
  97     BondedInteractions {FENE_bonds, eNR_FENEBONDS },              // F_FENEBONDS
  98     BondedInteractions {tab_bonds, eNR_TABBONDS },                // F_TABBONDS
  99     BondedInteractions {tab_bonds, eNR_TABBONDS },                // F_TABBONDSNC
 100     BondedInteractions {restraint_bonds, eNR_RESTRBONDS },        // F_RESTRBONDS
 101     BondedInteractions {angles, eNR_ANGLES },                     // F_ANGLES
 102     BondedInteractions {g96angles, eNR_ANGLES },                  // F_G96ANGLES
 103     BondedInteractions {restrangles, eNR_ANGLES },                // F_RESTRANGLES
 104     BondedInteractions {linear_angles, eNR_ANGLES },              // F_LINEAR_ANGLES
 105     BondedInteractions {cross_bond_bond, eNR_CROSS_BOND_BOND },   // F_CROSS_BOND_BONDS
 106     BondedInteractions {cross_bond_angle, eNR_CROSS_BOND_ANGLE }, // F_CROSS_BOND_ANGLES
 107     BondedInteractions {urey_bradley, eNR_UREY_BRADLEY },         // F_UREY_BRADLEY
 108     BondedInteractions {quartic_angles, eNR_QANGLES },            // F_QUARTIC_ANGLES
 109     BondedInteractions {tab_angles, eNR_TABANGLES },              // F_TABANGLES
 110     BondedInteractions {pdihs, eNR_PROPER },                      // F_PDIHS
 111     BondedInteractions {rbdihs, eNR_RB },                         // F_RBDIHS
 112     BondedInteractions {restrdihs, eNR_PROPER },                  // F_RESTRDIHS
 113     BondedInteractions {cbtdihs, eNR_RB },                        // F_CBTDIHS
 114     BondedInteractions {rbdihs, eNR_FOURDIH },                    // F_FOURDIHS
 115     BondedInteractions {idihs, eNR_IMPROPER },                    // F_IDIHS
 116     BondedInteractions {pdihs, eNR_IMPROPER },                    // F_PIDIHS
 117     BondedInteractions {tab_dihs, eNR_TABDIHS },                  // F_TABDIHS
 118     BondedInteractions {unimplemented, eNR_CMAP },                // F_CMAP
 119     BondedInteractions {unimplemented, -1 },                      // F_GB12_NOLONGERUSED
 120     BondedInteractions {unimplemented, -1 },                      // F_GB13_NOLONGERUSED
 121     BondedInteractions {unimplemented, -1 },                      // F_GB14_NOLONGERUSED
 122     BondedInteractions {unimplemented, -1 },                      // F_GBPOL_NOLONGERUSED
 123     BondedInteractions {unimplemented, -1 },                      // F_NPSOLVATION_NOLONGERUSED
 124     BondedInteractions {unimplemented, eNR_NB14 },                // F_LJ14
 125     BondedInteractions {unimplemented, -1 },                      // F_COUL14
 126     BondedInteractions {unimplemented, eNR_NB14 },                // F_LJC14_Q
 127     BondedInteractions {unimplemented, eNR_NB14 },                // F_LJC_PAIRS_NB
 128     BondedInteractions {unimplemented, -1 },                      // F_LJ
 129     BondedInteractions {unimplemented, -1 },                      // F_BHAM
 130     BondedInteractions {unimplemented, -1 },                      // F_LJ_LR_NOLONGERUSED
 131     BondedInteractions {unimplemented, -1 },                      // F_BHAM_LR_NOLONGERUSED
 132     BondedInteractions {unimplemented, -1 },                      // F_DISPCORR
 133     BondedInteractions {unimplemented, -1 },                      // F_COUL_SR
 134     BondedInteractions {unimplemented, -1 },                      // F_COUL_LR_NOLONGERUSED
 135     BondedInteractions {unimplemented, -1 },                      // F_RF_EXCL
 136     BondedInteractions {unimplemented, -1 },                      // F_COUL_RECIP
 137     BondedInteractions {unimplemented, -1 },                      // F_LJ_RECIP
 138     BondedInteractions {unimplemented, -1 },                      // F_DPD
 139     BondedInteractions {polarize, eNR_POLARIZE },                 // F_POLARIZATION
 140     BondedInteractions {water_pol, eNR_WPOL },                    // F_WATER_POL
 141     BondedInteractions {thole_pol, eNR_THOLE },                   // F_THOLE_POL
 142     BondedInteractions {anharm_polarize, eNR_ANHARM_POL },        // F_ANHARM_POL
 143     BondedInteractions {unimplemented, -1 },                      // F_POSRES
 144     BondedInteractions {unimplemented, -1 },                      // F_FBPOSRES
 145     BondedInteractions {ta_disres, eNR_DISRES },                  // F_DISRES
 146     BondedInteractions {unimplemented, -1 },                      // F_DISRESVIOL
 147     BondedInteractions {orires, eNR_ORIRES },                     // F_ORIRES
 148     BondedInteractions {unimplemented, -1 },                      // F_ORIRESDEV
 149     BondedInteractions {angres, eNR_ANGRES },                     // F_ANGRES
 150     BondedInteractions {angresz, eNR_ANGRESZ },                   // F_ANGRESZ
 151     BondedInteractions {dihres, eNR_DIHRES },                     // F_DIHRES
 152     BondedInteractions {unimplemented, -1 },                      // F_DIHRESVIOL
 153     BondedInteractions {unimplemented, -1 },                      // F_CONSTR
 154     BondedInteractions {unimplemented, -1 },                      // F_CONSTRNC
 155     BondedInteractions {unimplemented, -1 },                      // F_SETTLE
 156     BondedInteractions {unimplemented, -1 },                      // F_VSITE2
 157     BondedInteractions {unimplemented, -1 },                      // F_VSITE3
 158     BondedInteractions {unimplemented, -1 },                      // F_VSITE3FD
 159     BondedInteractions {unimplemented, -1 },                      // F_VSITE3FAD
 160     BondedInteractions {unimplemented, -1 },                      // F_VSITE3OUT
 161     BondedInteractions {unimplemented, -1 },                      // F_VSITE4FD
 162     BondedInteractions {unimplemented, -1 },                      // F_VSITE4FDN
 163     BondedInteractions {unimplemented, -1 },                      // F_VSITEN
 164     BondedInteractions {unimplemented, -1 },                      // F_COM_PULL
 165     BondedInteractions {unimplemented, -1 },                      // F_EQM
 166     BondedInteractions {unimplemented, -1 },                      // F_EPOT
 167     BondedInteractions {unimplemented, -1 },                      // F_EKIN
 168     BondedInteractions {unimplemented, -1 },                      // F_ETOT
 169     BondedInteractions {unimplemented, -1 },                      // F_ECONSERVED
 170     BondedInteractions {unimplemented, -1 },                      // F_TEMP
 171     BondedInteractions {unimplemented, -1 },                      // F_VTEMP_NOLONGERUSED
 172     BondedInteractions {unimplemented, -1 },                      // F_PDISPCORR
 173     BondedInteractions {unimplemented, -1 },                      // F_PRES
 174     BondedInteractions {unimplemented, -1 },                      // F_DVDL_CONSTR
 175     BondedInteractions {unimplemented, -1 },                      // F_DVDL
 176     BondedInteractions {unimplemented, -1 },                      // F_DKDL
 177     BondedInteractions {unimplemented, -1 },                      // F_DVDL_COUL
 178     BondedInteractions {unimplemented, -1 },                      // F_DVDL_VDW
 179     BondedInteractions {unimplemented, -1 },                      // F_DVDL_BONDED
 180     BondedInteractions {unimplemented, -1 },                      // F_DVDL_RESTRAINT
 181     BondedInteractions {unimplemented, -1 },                      // F_DVDL_TEMPERATURE
 182     };
 183
 184 BondedFunction bondedFunction(int ftype)
 185 {
 186     return s_bondedInteractionFunctions[ftype].function;
 187 }
 188
 189 //! Getter for finding the flop count for an \c ftype interaction.
 190 static int nrnbIndex(int ftype)
 191 {
 192     return s_bondedInteractionFunctions[ftype].nrnbIndex;
 193 }
 194
 195 namespace
 196 {
 197
 198 /*! \brief Return true if ftype is an explicit pair-listed LJ or
 199  * COULOMB interaction type: bonded LJ (usually 1-4), or special
 200  * listed non-bonded for FEP. */
 201 bool
 202 isPairInteraction(int ftype)
 203 {
 204     return ((ftype) >= F_LJ14 && (ftype) <= F_LJC_PAIRS_NB);
 205 }
 206
 207 /*! \brief Zero thread-local output buffers */
 208 void
 209 zero_thread_output(bonded_threading_t *bt, int thread)
 210 {
 211     if (!bt->haveBondeds)
 212     {
 213         return;
 214     }
 215
 216     f_thread_t *f_t      = &bt->f_t[thread];
 217     const int   nelem_fa = sizeof(f_t->f[0])/sizeof(real);
 218
 219     for (int i = 0; i < f_t->nblock_used; i++)
 220     {
 221         int a0 = f_t->block_index[i]*reduction_block_size;
 222         int a1 = a0 + reduction_block_size;
 223         for (int a = a0; a < a1; a++)
 224         {
 225             for (int d = 0; d < nelem_fa; d++)
 226             {
 227                 f_t->f[a][d] = 0;
 228             }
 229         }
 230     }
 231
 232     for (int i = 0; i < SHIFTS; i++)
 233     {
 234         clear_rvec(f_t->fshift[i]);
 235     }
 236     for (int i = 0; i < F_NRE; i++)
 237     {
 238         f_t->ener[i] = 0;
 239     }
 240     for (int i = 0; i < egNR; i++)
 241     {
 242         for (int j = 0; j < f_t->grpp.nener; j++)
 243         {
 244             f_t->grpp.ener[i][j] = 0;
 245         }
 246     }
 247     for (int i = 0; i < efptNR; i++)
 248     {
 249         f_t->dvdl[i] = 0;
 250     }
 251 }
 252
 253 /*! \brief The max thread number is arbitrary, we used a fixed number
 254  * to avoid memory management.  Using more than 16 threads is probably
 255  * never useful performance wise. */
 256 #define MAX_BONDED_THREADS 256
 257
 258 /*! \brief Reduce thread-local force buffers */
 259 void
 260 reduce_thread_forces(int n, rvec *f,
 261                      bonded_threading_t *bt,
 262                      int nthreads)
 263 {
 264     if (nthreads > MAX_BONDED_THREADS)
 265     {
 266         gmx_fatal(FARGS, "Can not reduce bonded forces on more than %d threads",
 267                   MAX_BONDED_THREADS);
 268     }
 269
 270     /* This reduction can run on any number of threads,
 271      * independently of bt->nthreads.
 272      * But if nthreads matches bt->nthreads (which it currently does)
 273      * the uniform distribution of the touched blocks over nthreads will
 274      * match the distribution of bonded over threads well in most cases,
 275      * which means that threads mostly reduce their own data which increases
 276      * the number of cache hits.
 277      */
 278 #pragma omp parallel for num_threads(nthreads) schedule(static)
 279     for (int b = 0; b < bt->nblock_used; b++)
 280     {
 281         try
 282         {
 283             int    ind = bt->block_index[b];
 284             rvec4 *fp[MAX_BONDED_THREADS];
 285
 286             /* Determine which threads contribute to this block */
 287             int nfb = 0;
 288             for (int ft = 0; ft < bt->nthreads; ft++)
 289             {
 290                 if (bitmask_is_set(bt->mask[ind], ft))
 291                 {
 292                     fp[nfb++] = bt->f_t[ft].f;
 293                 }
 294             }
 295             if (nfb > 0)
 296             {
 297                 /* Reduce force buffers for threads that contribute */
 298                 int a0 =  ind     *reduction_block_size;
 299                 int a1 = (ind + 1)*reduction_block_size;
 300                 /* It would be nice if we could pad f to avoid this min */
 301                 a1     = std::min(a1, n);
 302                 for (int a = a0; a < a1; a++)
 303                 {
 304                     for (int fb = 0; fb < nfb; fb++)
 305                     {
 306                         rvec_inc(f[a], fp[fb][a]);
 307                     }
 308                 }
 309             }
 310         }
 311         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 312     }
 313 }
 314
 315 /*! \brief Reduce thread-local forces, shift forces and energies */
 316 void
 317 reduce_thread_output(int n, rvec *f, rvec *fshift,
 318                      real *ener, gmx_grppairener_t *grpp, real *dvdl,
 319                      bonded_threading_t *bt,
 320                      gmx_bool bCalcEnerVir,
 321                      gmx_bool bDHDL)
 322 {
 323     assert(bt->haveBondeds);
 324
 325     if (bt->nblock_used > 0)
 326     {
 327         /* Reduce the bonded force buffer */
 328         reduce_thread_forces(n, f, bt, bt->nthreads);
 329     }
 330
 331     /* When necessary, reduce energy and virial using one thread only */
 332     if (bCalcEnerVir && bt->nthreads > 1)
 333     {
 334         f_thread_t *f_t = bt->f_t;
 335
 336         for (int i = 0; i < SHIFTS; i++)
 337         {
 338             for (int t = 1; t < bt->nthreads; t++)
 339             {
 340                 rvec_inc(fshift[i], f_t[t].fshift[i]);
 341             }
 342         }
 343         for (int i = 0; i < F_NRE; i++)
 344         {
 345             for (int t = 1; t < bt->nthreads; t++)
 346             {
 347                 ener[i] += f_t[t].ener[i];
 348             }
 349         }
 350         for (int i = 0; i < egNR; i++)
 351         {
 352             for (int j = 0; j < f_t[1].grpp.nener; j++)
 353             {
 354                 for (int t = 1; t < bt->nthreads; t++)
 355                 {
 356                     grpp->ener[i][j] += f_t[t].grpp.ener[i][j];
 357                 }
 358             }
 359         }
 360         if (bDHDL)
 361         {
 362             for (int i = 0; i < efptNR; i++)
 363             {
 364
 365                 for (int t = 1; t < bt->nthreads; t++)
 366                 {
 367                     dvdl[i] += f_t[t].dvdl[i];
 368                 }
 369             }
 370         }
 371     }
 372 }
 373
 374 /*! \brief Calculate one element of the list of bonded interactions
 375     for this thread */
 376 real
 377 calc_one_bond(int thread,
 378               int ftype, const t_idef *idef,
 379               const bonded_threading_t &bondedThreading,
 380               const rvec x[], rvec4 f[], rvec fshift[],
 381               const t_forcerec *fr,
 382               const t_pbc *pbc, const t_graph *g,
 383               gmx_grppairener_t *grpp,
 384               t_nrnb *nrnb,
 385               const real *lambda, real *dvdl,
 386               const t_mdatoms *md, t_fcdata *fcd,
 387               gmx_bool bCalcEnerVir,
 388               int *global_atom_index)
 389 {
 390 #if GMX_SIMD_HAVE_REAL
 391     bool bUseSIMD = fr->use_simd_kernels;
 392 #endif
 393
 394     int      nat1, nbonds, efptFTYPE;
 395     real     v = 0;
 396     t_iatom *iatoms;
 397     int      nb0, nbn;
 398
 399     if (IS_RESTRAINT_TYPE(ftype))
 400     {
 401         efptFTYPE = efptRESTRAINT;
 402     }
 403     else
 404     {
 405         efptFTYPE = efptBONDED;
 406     }
 407
 408     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");
 409     const bool useFreeEnergy     = (idef->ilsort == ilsortFE_SORTED && idef->il[ftype].nr_nonperturbed < idef->il[ftype].nr);
 410     const bool computeForcesOnly = (!bCalcEnerVir && !useFreeEnergy);
 411
 412     nat1      = interaction_function[ftype].nratoms + 1;
 413     nbonds    = idef->il[ftype].nr/nat1;
 414     iatoms    = idef->il[ftype].iatoms;
 415
 416     GMX_ASSERT(fr->gpuBonded != nullptr || bondedThreading.il_thread_division[ftype*(bondedThreading.nthreads + 1) + bondedThreading.nthreads] == idef->il[ftype].nr, "The thread division should match the topology");
 417
 418     nb0 = bondedThreading.il_thread_division[ftype*(bondedThreading.nthreads+1)+thread];
 419     nbn = bondedThreading.il_thread_division[ftype*(bondedThreading.nthreads+1)+thread+1] - nb0;
 420
 421     if (!isPairInteraction(ftype))
 422     {
 423         if (ftype == F_CMAP)
 424         {
 425             /* TODO The execution time for CMAP dihedrals might be
 426                nice to account to its own subtimer, but first
 427                wallcycle needs to be extended to support calling from
 428                multiple threads. */
 429             v = cmap_dihs(nbn, iatoms+nb0,
 430                           idef->iparams, idef->cmap_grid,
 431                           x, f, fshift,
 432                           pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]),
 433                           md, fcd, global_atom_index);
 434         }
 435 #if GMX_SIMD_HAVE_REAL
 436         else if (ftype == F_ANGLES && bUseSIMD && computeForcesOnly)
 437         {
 438             /* No energies, shift forces, dvdl */
 439             angles_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
 440                                idef->iparams,
 441                                x, f,
 442                                pbc, g, lambda[efptFTYPE], md, fcd,
 443                                global_atom_index);
 444             v = 0;
 445         }
 446
 447         else if (ftype == F_UREY_BRADLEY && bUseSIMD && computeForcesOnly)
 448         {
 449             /* No energies, shift forces, dvdl */
 450             urey_bradley_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
 451                                      idef->iparams,
 452                                      x, f,
 453                                      pbc, g, lambda[efptFTYPE], md, fcd,
 454                                      global_atom_index);
 455             v = 0;
 456         }
 457 #endif
 458         else if (ftype == F_PDIHS && computeForcesOnly)
 459         {
 460             /* No energies, shift forces, dvdl */
 461 #if GMX_SIMD_HAVE_REAL
 462             if (bUseSIMD)
 463             {
 464                 pdihs_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
 465                                   idef->iparams,
 466                                   x, f,
 467                                   pbc, g, lambda[efptFTYPE], md, fcd,
 468                                   global_atom_index);
 469             }
 470             else
 471 #endif
 472             {
 473                 pdihs_noener(nbn, idef->il[ftype].iatoms+nb0,
 474                              idef->iparams,
 475                              x, f,
 476                              pbc, g, lambda[efptFTYPE], md, fcd,
 477                              global_atom_index);
 478             }
 479             v = 0;
 480         }
 481 #if GMX_SIMD_HAVE_REAL
 482         else if (ftype == F_RBDIHS && bUseSIMD && computeForcesOnly)
 483         {
 484             /* No energies, shift forces, dvdl */
 485             rbdihs_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
 486                                idef->iparams,
 487                                static_cast<const rvec*>(x), f,
 488                                pbc, g, lambda[efptFTYPE], md, fcd,
 489                                global_atom_index);
 490             v = 0;
 491         }
 492 #endif
 493         else
 494         {
 495             v = bondedFunction(ftype)(nbn, iatoms+nb0,
 496                                       idef->iparams,
 497                                       x, f, fshift,
 498                                       pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]),
 499                                       md, fcd, global_atom_index);
 500         }
 501     }
 502     else
 503     {
 504         /* TODO The execution time for pairs might be nice to account
 505            to its own subtimer, but first wallcycle needs to be
 506            extended to support calling from multiple threads. */
 507         do_pairs(ftype, nbn, iatoms+nb0, idef->iparams, x, f, fshift,
 508                  pbc, g, lambda, dvdl, md, fr,
 509                  computeForcesOnly, grpp, global_atom_index);
 510         v = 0;
 511     }
 512
 513     if (thread == 0)
 514     {
 515         inc_nrnb(nrnb, nrnbIndex(ftype), nbonds);
 516     }
 517
 518     return v;
 519 }
 520
 521 } // namespace
 522
 523 /*! \brief Compute the bonded part of the listed forces, parallelized over threads
 524  */
 525 static void
 526 calcBondedForces(const t_idef     *idef,
 527                  const rvec        x[],
 528                  const t_forcerec *fr,
 529                  const t_pbc      *pbc_null,
 530                  const t_graph    *g,
 531                  gmx_enerdata_t   *enerd,
 532                  t_nrnb           *nrnb,
 533                  const real       *lambda,
 534                  real             *dvdl,
 535                  const t_mdatoms  *md,
 536                  t_fcdata         *fcd,
 537                  gmx_bool          bCalcEnerVir,
 538                  int              *global_atom_index)
 539 {
 540     bonded_threading_t *bt = fr->bondedThreading;
 541
 542 #pragma omp parallel for num_threads(bt->nthreads) schedule(static)
 543     for (int thread = 0; thread < bt->nthreads; thread++)
 544     {
 545         try
 546         {
 547             int                ftype;
 548             real              *epot, v;
 549             /* thread stuff */
 550             rvec4             *ft;
 551             rvec              *fshift;
 552             real              *dvdlt;
 553             gmx_grppairener_t *grpp;
 554
 555             zero_thread_output(bt, thread);
 556
 557             ft = bt->f_t[thread].f;
 558
 559             if (thread == 0)
 560             {
 561                 fshift = fr->fshift;
 562                 epot   = enerd->term;
 563                 grpp   = &enerd->grpp;
 564                 dvdlt  = dvdl;
 565             }
 566             else
 567             {
 568                 fshift = bt->f_t[thread].fshift;
 569                 epot   = bt->f_t[thread].ener;
 570                 grpp   = &bt->f_t[thread].grpp;
 571                 dvdlt  = bt->f_t[thread].dvdl;
 572             }
 573             /* Loop over all bonded force types to calculate the bonded forces */
 574             for (ftype = 0; (ftype < F_NRE); ftype++)
 575             {
 576                 if (idef->il[ftype].nr > 0 && ftype_is_bonded_potential(ftype))
 577                 {
 578                     v = calc_one_bond(thread, ftype, idef,
 579                                       *fr->bondedThreading, x,
 580                                       ft, fshift, fr, pbc_null, g, grpp,
 581                                       nrnb, lambda, dvdlt,
 582                                       md, fcd, bCalcEnerVir,
 583                                       global_atom_index);
 584                     epot[ftype] += v;
 585                 }
 586             }
 587         }
 588         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 589     }
 590 }
 591
 592 void calc_listed(const t_commrec             *cr,
 593                  const gmx_multisim_t *ms,
 594                  struct gmx_wallcycle        *wcycle,
 595                  const t_idef *idef,
 596                  const rvec x[], history_t *hist,
 597                  rvec f[],
 598                  gmx::ForceWithVirial *forceWithVirial,
 599                  const t_forcerec *fr,
 600                  const struct t_pbc *pbc,
 601                  const struct t_pbc *pbc_full,
 602                  const struct t_graph *g,
 603                  gmx_enerdata_t *enerd, t_nrnb *nrnb,
 604                  const real *lambda,
 605                  const t_mdatoms *md,
 606                  t_fcdata *fcd, int *global_atom_index,
 607                  int force_flags)
 608 {
 609     gmx_bool                   bCalcEnerVir;
 610     const  t_pbc              *pbc_null;
 611     bonded_threading_t        *bt  = fr->bondedThreading;
 612
 613     bCalcEnerVir = ((force_flags & (GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY)) != 0);
 614
 615     if (fr->bMolPBC)
 616     {
 617         pbc_null = pbc;
 618     }
 619     else
 620     {
 621         pbc_null = nullptr;
 622     }
 623
 624     if ((idef->il[F_POSRES].nr > 0) ||
 625         (idef->il[F_FBPOSRES].nr > 0) ||
 626         fcd->orires.nr > 0 ||
 627         fcd->disres.nres > 0)
 628     {
 629         /* TODO Use of restraints triggers further function calls
 630            inside the loop over calc_one_bond(), but those are too
 631            awkward to account to this subtimer properly in the present
 632            code. We don't test / care much about performance with
 633            restraints, anyway. */
 634         wallcycle_sub_start(wcycle, ewcsRESTRAINTS);
 635
 636         if (idef->il[F_POSRES].nr > 0)
 637         {
 638             posres_wrapper(nrnb, idef, pbc_full, x, enerd, lambda, fr,
 639                            forceWithVirial);
 640         }
 641
 642         if (idef->il[F_FBPOSRES].nr > 0)
 643         {
 644             fbposres_wrapper(nrnb, idef, pbc_full, x, enerd, fr,
 645                              forceWithVirial);
 646         }
 647
 648         /* Do pre force calculation stuff which might require communication */
 649         if (fcd->orires.nr > 0)
 650         {
 651             /* This assertion is to ensure we have whole molecules.
 652              * Unfortunately we do not have an mdrun state variable that tells
 653              * us if molecules in x are not broken over PBC, so we have to make
 654              * do with checking graph!=nullptr, which should tell us if we made
 655              * molecules whole before calling the current function.
 656              */
 657             GMX_RELEASE_ASSERT(fr->ePBC == epbcNONE || g != nullptr, "With orientation restraints molecules should be whole");
 658             enerd->term[F_ORIRESDEV] =
 659                 calc_orires_dev(ms, idef->il[F_ORIRES].nr,
 660                                 idef->il[F_ORIRES].iatoms,
 661                                 idef->iparams, md, x,
 662                                 pbc_null, fcd, hist);
 663         }
 664         if (fcd->disres.nres > 0)
 665         {
 666             calc_disres_R_6(cr, ms,
 667                             idef->il[F_DISRES].nr,
 668                             idef->il[F_DISRES].iatoms,
 669                             x, pbc_null,
 670                             fcd, hist);
 671         }
 672
 673         wallcycle_sub_stop(wcycle, ewcsRESTRAINTS);
 674     }
 675
 676     if (bt->haveBondeds)
 677     {
 678         wallcycle_sub_start(wcycle, ewcsLISTED);
 679         /* The dummy array is to have a place to store the dhdl at other values
 680            of lambda, which will be thrown away in the end */
 681         real dvdl[efptNR] = {0};
 682         calcBondedForces(idef, x, fr, pbc_null, g, enerd, nrnb, lambda, dvdl, md,
 683                          fcd, bCalcEnerVir, global_atom_index);
 684         wallcycle_sub_stop(wcycle, ewcsLISTED);
 685
 686         wallcycle_sub_start(wcycle, ewcsLISTED_BUF_OPS);
 687         reduce_thread_output(fr->natoms_force, f, fr->fshift,
 688                              enerd->term, &enerd->grpp, dvdl,
 689                              bt,
 690                              bCalcEnerVir,
 691                              (force_flags & GMX_FORCE_DHDL) != 0);
 692
 693         if (force_flags & GMX_FORCE_DHDL)
 694         {
 695             for (int i = 0; i < efptNR; i++)
 696             {
 697                 enerd->dvdl_nonlin[i] += dvdl[i];
 698             }
 699         }
 700         wallcycle_sub_stop(wcycle, ewcsLISTED_BUF_OPS);
 701     }
 702
 703     /* Copy the sum of violations for the distance restraints from fcd */
 704     if (fcd)
 705     {
 706         enerd->term[F_DISRESVIOL] = fcd->disres.sumviol;
 707     }
 708 }
 709
 710 void calc_listed_lambda(const t_idef *idef,
 711                         const rvec x[],
 712                         const t_forcerec *fr,
 713                         const struct t_pbc *pbc, const struct t_graph *g,
 714                         gmx_grppairener_t *grpp, real *epot, t_nrnb *nrnb,
 715                         const real *lambda,
 716                         const t_mdatoms *md,
 717                         t_fcdata *fcd,
 718                         int *global_atom_index)
 719 {
 720     real               v;
 721     real               dvdl_dum[efptNR] = {0};
 722     rvec4             *f;
 723     rvec              *fshift;
 724     const  t_pbc      *pbc_null;
 725     t_idef             idef_fe;
 726     bonded_threading_t bondedThreading;
 727
 728     if (fr->bMolPBC)
 729     {
 730         pbc_null = pbc;
 731     }
 732     else
 733     {
 734         pbc_null = nullptr;
 735     }
 736
 737     /* Copy the whole idef, so we can modify the contents locally */
 738     idef_fe                  = *idef;
 739     bondedThreading.nthreads = 1;
 740     snew(bondedThreading.il_thread_division, F_NRE*(bondedThreading.nthreads+1));
 741
 742     /* We already have the forces, so we use temp buffers here */
 743     snew(f, fr->natoms_force);
 744     snew(fshift, SHIFTS);
 745
 746     /* Loop over all bonded force types to calculate the bonded energies */
 747     for (int ftype = 0; (ftype < F_NRE); ftype++)
 748     {
 749         if (ftype_is_bonded_potential(ftype))
 750         {
 751             const t_ilist &ilist     = idef->il[ftype];
 752             /* Create a temporary t_ilist with only perturbed interactions */
 753             t_ilist       &ilist_fe  = idef_fe.il[ftype];
 754             ilist_fe.iatoms          = ilist.iatoms + ilist.nr_nonperturbed;
 755             ilist_fe.nr_nonperturbed = 0;
 756             ilist_fe.nr              = ilist.nr - ilist.nr_nonperturbed;
 757             /* Set the work range of thread 0 to the perturbed bondeds */
 758             bondedThreading.il_thread_division[ftype*2 + 0] = 0;
 759             bondedThreading.il_thread_division[ftype*2 + 1] = ilist_fe.nr;
 760
 761             if (ilist_fe.nr > 0)
 762             {
 763                 v = calc_one_bond(0, ftype, &idef_fe, bondedThreading,
 764                                   x, f, fshift, fr, pbc_null, g,
 765                                   grpp, nrnb, lambda, dvdl_dum,
 766                                   md, fcd, TRUE,
 767                                   global_atom_index);
 768                 epot[ftype] += v;
 769             }
 770         }
 771     }
 772
 773     sfree(fshift);
 774     sfree(f);
 775
 776     sfree(bondedThreading.il_thread_division);
 777 }
 778
 779 void
 780 do_force_listed(struct gmx_wallcycle        *wcycle,
 781                 matrix                       box,
 782                 const t_lambda              *fepvals,
 783                 const t_commrec             *cr,
 784                 const gmx_multisim_t        *ms,
 785                 const t_idef                *idef,
 786                 const rvec                   x[],
 787                 history_t                   *hist,
 788                 rvec                        *forceForUseWithShiftForces,
 789                 gmx::ForceWithVirial        *forceWithVirial,
 790                 const t_forcerec            *fr,
 791                 const struct t_pbc          *pbc,
 792                 const struct t_graph        *graph,
 793                 gmx_enerdata_t              *enerd,
 794                 t_nrnb                      *nrnb,
 795                 const real                  *lambda,
 796                 const t_mdatoms             *md,
 797                 t_fcdata                    *fcd,
 798                 int                         *global_atom_index,
 799                 int                          flags)
 800 {
 801     t_pbc pbc_full; /* Full PBC is needed for position restraints */
 802
 803     if (!(flags & GMX_FORCE_LISTED))
 804     {
 805         return;
 806     }
 807
 808     if ((idef->il[F_POSRES].nr > 0) ||
 809         (idef->il[F_FBPOSRES].nr > 0))
 810     {
 811         /* Not enough flops to bother counting */
 812         set_pbc(&pbc_full, fr->ePBC, box);
 813     }
 814     calc_listed(cr, ms, wcycle, idef, x, hist,
 815                 forceForUseWithShiftForces, forceWithVirial,
 816                 fr, pbc, &pbc_full,
 817                 graph, enerd, nrnb, lambda, md, fcd,
 818                 global_atom_index, flags);
 819
 820     /* Check if we have to determine energy differences
 821      * at foreign lambda's.
 822      */
 823     if (fepvals->n_lambda > 0 && (flags & GMX_FORCE_DHDL))
 824     {
 825         posres_wrapper_lambda(wcycle, fepvals, idef, &pbc_full, x, enerd, lambda, fr);
 826
 827         if (idef->ilsort != ilsortNO_FE)
 828         {
 829             wallcycle_sub_start(wcycle, ewcsLISTED_FEP);
 830             if (idef->ilsort != ilsortFE_SORTED)
 831             {
 832                 gmx_incons("The bonded interactions are not sorted for free energy");
 833             }
 834             for (int i = 0; i < enerd->n_lambda; i++)
 835             {
 836                 real lam_i[efptNR];
 837
 838                 reset_foreign_enerdata(enerd);
 839                 for (int j = 0; j < efptNR; j++)
 840                 {
 841                     lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i-1]);
 842                 }
 843                 calc_listed_lambda(idef, x, fr, pbc, graph, &(enerd->foreign_grpp), enerd->foreign_term, nrnb, lam_i, md,
 844                                    fcd, global_atom_index);
 845                 sum_epot(&(enerd->foreign_grpp), enerd->foreign_term);
 846                 enerd->enerpart_lambda[i] += enerd->foreign_term[F_EPOT];
 847             }
 848             wallcycle_sub_stop(wcycle, ewcsLISTED_FEP);
 849         }
 850     }
 851 }