2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
37 * \brief This file defines high-level functions for mdrun to compute
38 * energies and forces for listed interactions.
40 * \author Mark Abraham <mark.j.abraham@gmail.com>
42 * \ingroup module_listed_forces
46 #include "listed_forces.h"
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"
77 #include "listed-internal.h"
78 #include "utilities.h"
80 struct BondedInteractions
82 BondedFunction function;
86 /*! \brief Lookup table of bonded interaction functions
88 * This must have as many entries as interaction_function in ifunc.cpp */
89 static std::array<BondedInteractions, F_NRE> s_bondedInteractionFunctions
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
184 BondedFunction bondedFunction(int ftype)
186 return s_bondedInteractionFunctions[ftype].function;
189 //! Getter for finding the flop count for an \c ftype interaction.
190 static int nrnbIndex(int ftype)
192 return s_bondedInteractionFunctions[ftype].nrnbIndex;
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. */
202 isPairInteraction(int ftype)
204 return ((ftype) >= F_LJ14 && (ftype) <= F_LJC_PAIRS_NB);
207 /*! \brief Zero thread-local output buffers */
209 zero_thread_output(bonded_threading_t *bt, int thread)
211 if (!bt->haveBondeds)
216 f_thread_t *f_t = &bt->f_t[thread];
217 const int nelem_fa = sizeof(f_t->f[0])/sizeof(real);
219 for (int i = 0; i < f_t->nblock_used; i++)
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++)
225 for (int d = 0; d < nelem_fa; d++)
232 for (int i = 0; i < SHIFTS; i++)
234 clear_rvec(f_t->fshift[i]);
236 for (int i = 0; i < F_NRE; i++)
240 for (int i = 0; i < egNR; i++)
242 for (int j = 0; j < f_t->grpp.nener; j++)
244 f_t->grpp.ener[i][j] = 0;
247 for (int i = 0; i < efptNR; i++)
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
258 /*! \brief Reduce thread-local force buffers */
260 reduce_thread_forces(int n, rvec *f,
261 bonded_threading_t *bt,
264 if (nthreads > MAX_BONDED_THREADS)
266 gmx_fatal(FARGS, "Can not reduce bonded forces on more than %d threads",
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.
278 #pragma omp parallel for num_threads(nthreads) schedule(static)
279 for (int b = 0; b < bt->nblock_used; b++)
283 int ind = bt->block_index[b];
284 rvec4 *fp[MAX_BONDED_THREADS];
286 /* Determine which threads contribute to this block */
288 for (int ft = 0; ft < bt->nthreads; ft++)
290 if (bitmask_is_set(bt->mask[ind], ft))
292 fp[nfb++] = bt->f_t[ft].f;
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++)
304 for (int fb = 0; fb < nfb; fb++)
306 rvec_inc(f[a], fp[fb][a]);
311 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
315 /*! \brief Reduce thread-local forces, shift forces and energies */
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,
323 assert(bt->haveBondeds);
325 if (bt->nblock_used > 0)
327 /* Reduce the bonded force buffer */
328 reduce_thread_forces(n, f, bt, bt->nthreads);
331 /* When necessary, reduce energy and virial using one thread only */
332 if (bCalcEnerVir && bt->nthreads > 1)
334 f_thread_t *f_t = bt->f_t;
336 for (int i = 0; i < SHIFTS; i++)
338 for (int t = 1; t < bt->nthreads; t++)
340 rvec_inc(fshift[i], f_t[t].fshift[i]);
343 for (int i = 0; i < F_NRE; i++)
345 for (int t = 1; t < bt->nthreads; t++)
347 ener[i] += f_t[t].ener[i];
350 for (int i = 0; i < egNR; i++)
352 for (int j = 0; j < f_t[1].grpp.nener; j++)
354 for (int t = 1; t < bt->nthreads; t++)
356 grpp->ener[i][j] += f_t[t].grpp.ener[i][j];
362 for (int i = 0; i < efptNR; i++)
365 for (int t = 1; t < bt->nthreads; t++)
367 dvdl[i] += f_t[t].dvdl[i];
374 /*! \brief Calculate one element of the list of bonded interactions
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,
385 const real *lambda, real *dvdl,
386 const t_mdatoms *md, t_fcdata *fcd,
387 gmx_bool bCalcEnerVir,
388 int *global_atom_index)
390 #if GMX_SIMD_HAVE_REAL
391 bool bUseSIMD = fr->use_simd_kernels;
394 int nat1, nbonds, efptFTYPE;
399 if (IS_RESTRAINT_TYPE(ftype))
401 efptFTYPE = efptRESTRAINT;
405 efptFTYPE = efptBONDED;
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);
412 nat1 = interaction_function[ftype].nratoms + 1;
413 nbonds = idef->il[ftype].nr/nat1;
414 iatoms = idef->il[ftype].iatoms;
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");
418 nb0 = bondedThreading.il_thread_division[ftype*(bondedThreading.nthreads+1)+thread];
419 nbn = bondedThreading.il_thread_division[ftype*(bondedThreading.nthreads+1)+thread+1] - nb0;
421 if (!isPairInteraction(ftype))
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
429 v = cmap_dihs(nbn, iatoms+nb0,
430 idef->iparams, idef->cmap_grid,
432 pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]),
433 md, fcd, global_atom_index);
435 #if GMX_SIMD_HAVE_REAL
436 else if (ftype == F_ANGLES && bUseSIMD && computeForcesOnly)
438 /* No energies, shift forces, dvdl */
439 angles_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
442 pbc, g, lambda[efptFTYPE], md, fcd,
447 else if (ftype == F_UREY_BRADLEY && bUseSIMD && computeForcesOnly)
449 /* No energies, shift forces, dvdl */
450 urey_bradley_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
453 pbc, g, lambda[efptFTYPE], md, fcd,
458 else if (ftype == F_PDIHS && computeForcesOnly)
460 /* No energies, shift forces, dvdl */
461 #if GMX_SIMD_HAVE_REAL
464 pdihs_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
467 pbc, g, lambda[efptFTYPE], md, fcd,
473 pdihs_noener(nbn, idef->il[ftype].iatoms+nb0,
476 pbc, g, lambda[efptFTYPE], md, fcd,
481 #if GMX_SIMD_HAVE_REAL
482 else if (ftype == F_RBDIHS && bUseSIMD && computeForcesOnly)
484 /* No energies, shift forces, dvdl */
485 rbdihs_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
487 static_cast<const rvec*>(x), f,
488 pbc, g, lambda[efptFTYPE], md, fcd,
495 v = bondedFunction(ftype)(nbn, iatoms+nb0,
498 pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]),
499 md, fcd, global_atom_index);
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);
515 inc_nrnb(nrnb, nrnbIndex(ftype), nbonds);
523 /*! \brief Compute the bonded part of the listed forces, parallelized over threads
526 calcBondedForces(const t_idef *idef,
528 const t_forcerec *fr,
529 const t_pbc *pbc_null,
531 gmx_enerdata_t *enerd,
537 gmx_bool bCalcEnerVir,
538 int *global_atom_index)
540 bonded_threading_t *bt = fr->bondedThreading;
542 #pragma omp parallel for num_threads(bt->nthreads) schedule(static)
543 for (int thread = 0; thread < bt->nthreads; thread++)
553 gmx_grppairener_t *grpp;
555 zero_thread_output(bt, thread);
557 ft = bt->f_t[thread].f;
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;
573 /* Loop over all bonded force types to calculate the bonded forces */
574 for (ftype = 0; (ftype < F_NRE); ftype++)
576 if (idef->il[ftype].nr > 0 && ftype_is_bonded_potential(ftype))
578 v = calc_one_bond(thread, ftype, idef,
579 *fr->bondedThreading, x,
580 ft, fshift, fr, pbc_null, g, grpp,
582 md, fcd, bCalcEnerVir,
588 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
592 void calc_listed(const t_commrec *cr,
593 const gmx_multisim_t *ms,
594 struct gmx_wallcycle *wcycle,
596 const rvec x[], history_t *hist,
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,
606 t_fcdata *fcd, int *global_atom_index,
609 gmx_bool bCalcEnerVir;
610 const t_pbc *pbc_null;
611 bonded_threading_t *bt = fr->bondedThreading;
613 bCalcEnerVir = ((force_flags & (GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY)) != 0);
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)
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);
636 if (idef->il[F_POSRES].nr > 0)
638 posres_wrapper(nrnb, idef, pbc_full, x, enerd, lambda, fr,
642 if (idef->il[F_FBPOSRES].nr > 0)
644 fbposres_wrapper(nrnb, idef, pbc_full, x, enerd, fr,
648 /* Do pre force calculation stuff which might require communication */
649 if (fcd->orires.nr > 0)
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.
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);
664 if (fcd->disres.nres > 0)
666 calc_disres_R_6(cr, ms,
667 idef->il[F_DISRES].nr,
668 idef->il[F_DISRES].iatoms,
673 wallcycle_sub_stop(wcycle, ewcsRESTRAINTS);
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);
686 wallcycle_sub_start(wcycle, ewcsLISTED_BUF_OPS);
687 reduce_thread_output(fr->natoms_force, f, fr->fshift,
688 enerd->term, &enerd->grpp, dvdl,
691 (force_flags & GMX_FORCE_DHDL) != 0);
693 if (force_flags & GMX_FORCE_DHDL)
695 for (int i = 0; i < efptNR; i++)
697 enerd->dvdl_nonlin[i] += dvdl[i];
700 wallcycle_sub_stop(wcycle, ewcsLISTED_BUF_OPS);
703 /* Copy the sum of violations for the distance restraints from fcd */
706 enerd->term[F_DISRESVIOL] = fcd->disres.sumviol;
710 void calc_listed_lambda(const t_idef *idef,
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,
718 int *global_atom_index)
721 real dvdl_dum[efptNR] = {0};
724 const t_pbc *pbc_null;
726 bonded_threading_t bondedThreading;
737 /* Copy the whole idef, so we can modify the contents locally */
739 bondedThreading.nthreads = 1;
740 snew(bondedThreading.il_thread_division, F_NRE*(bondedThreading.nthreads+1));
742 /* We already have the forces, so we use temp buffers here */
743 snew(f, fr->natoms_force);
744 snew(fshift, SHIFTS);
746 /* Loop over all bonded force types to calculate the bonded energies */
747 for (int ftype = 0; (ftype < F_NRE); ftype++)
749 if (ftype_is_bonded_potential(ftype))
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;
763 v = calc_one_bond(0, ftype, &idef_fe, bondedThreading,
764 x, f, fshift, fr, pbc_null, g,
765 grpp, nrnb, lambda, dvdl_dum,
776 sfree(bondedThreading.il_thread_division);
780 do_force_listed(struct gmx_wallcycle *wcycle,
782 const t_lambda *fepvals,
784 const gmx_multisim_t *ms,
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,
798 int *global_atom_index,
801 t_pbc pbc_full; /* Full PBC is needed for position restraints */
803 if (!(flags & GMX_FORCE_LISTED))
808 if ((idef->il[F_POSRES].nr > 0) ||
809 (idef->il[F_FBPOSRES].nr > 0))
811 /* Not enough flops to bother counting */
812 set_pbc(&pbc_full, fr->ePBC, box);
814 calc_listed(cr, ms, wcycle, idef, x, hist,
815 forceForUseWithShiftForces, forceWithVirial,
817 graph, enerd, nrnb, lambda, md, fcd,
818 global_atom_index, flags);
820 /* Check if we have to determine energy differences
821 * at foreign lambda's.
823 if (fepvals->n_lambda > 0 && (flags & GMX_FORCE_DHDL))
825 posres_wrapper_lambda(wcycle, fepvals, idef, &pbc_full, x, enerd, lambda, fr);
827 if (idef->ilsort != ilsortNO_FE)
829 wallcycle_sub_start(wcycle, ewcsLISTED_FEP);
830 if (idef->ilsort != ilsortFE_SORTED)
832 gmx_incons("The bonded interactions are not sorted for free energy");
834 for (int i = 0; i < enerd->n_lambda; i++)
838 reset_foreign_enerdata(enerd);
839 for (int j = 0; j < efptNR; j++)
841 lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i-1]);
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];
848 wallcycle_sub_stop(wcycle, ewcsLISTED_FEP);