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, 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 "config.h"
  49
  50 #include <assert.h>
  51
  52 #include <algorithm>
  53
  54 #include "gromacs/legacyheaders/disre.h"
  55 #include "gromacs/legacyheaders/force.h"
  56 #include "gromacs/legacyheaders/network.h"
  57 #include "gromacs/legacyheaders/nrnb.h"
  58 #include "gromacs/legacyheaders/orires.h"
  59 #include "gromacs/legacyheaders/types/force_flags.h"
  60 #include "gromacs/listed-forces/bonded.h"
  61 #include "gromacs/listed-forces/position-restraints.h"
  62 #include "gromacs/math/vec.h"
  63 #include "gromacs/mdlib/forcerec-threading.h"
  64 #include "gromacs/pbcutil/ishift.h"
  65 #include "gromacs/pbcutil/pbc.h"
  66 #include "gromacs/simd/simd.h"
  67 #include "gromacs/timing/wallcycle.h"
  68 #include "gromacs/utility/smalloc.h"
  69
  70 #include "pairs.h"
  71
  72 namespace
  73 {
  74
  75 /*! \brief Return true if ftype is an explicit pair-listed LJ or
  76  * COULOMB interaction type: bonded LJ (usually 1-4), or special
  77  * listed non-bonded for FEP. */
  78 bool
  79 isPairInteraction(int ftype)
  80 {
  81     return ((ftype) >= F_LJ14 && (ftype) <= F_LJC_PAIRS_NB);
  82 }
  83
  84 /*! \brief Zero thread-local force-output buffers */
  85 void
  86 zero_thread_forces(f_thread_t *f_t, int n,
  87                    int nblock, int blocksize)
  88 {
  89     int b, a0, a1, a, i, j;
  90
  91     if (n > f_t->f_nalloc)
  92     {
  93         f_t->f_nalloc = over_alloc_large(n);
  94         srenew(f_t->f, f_t->f_nalloc);
  95     }
  96
  97     if (!bitmask_is_zero(f_t->red_mask))
  98     {
  99         for (b = 0; b < nblock; b++)
 100         {
 101             if (bitmask_is_set(f_t->red_mask, b))
 102             {
 103                 a0 = b*blocksize;
 104                 a1 = std::min((b+1)*blocksize, n);
 105                 for (a = a0; a < a1; a++)
 106                 {
 107                     clear_rvec(f_t->f[a]);
 108                 }
 109             }
 110         }
 111     }
 112     for (i = 0; i < SHIFTS; i++)
 113     {
 114         clear_rvec(f_t->fshift[i]);
 115     }
 116     for (i = 0; i < F_NRE; i++)
 117     {
 118         f_t->ener[i] = 0;
 119     }
 120     for (i = 0; i < egNR; i++)
 121     {
 122         for (j = 0; j < f_t->grpp.nener; j++)
 123         {
 124             f_t->grpp.ener[i][j] = 0;
 125         }
 126     }
 127     for (i = 0; i < efptNR; i++)
 128     {
 129         f_t->dvdl[i] = 0;
 130     }
 131 }
 132
 133 /*! \brief The max thread number is arbitrary, we used a fixed number
 134  * to avoid memory management.  Using more than 16 threads is probably
 135  * never useful performance wise. */
 136 #define MAX_BONDED_THREADS 256
 137
 138 /*! \brief Reduce thread-local force buffers */
 139 void
 140 reduce_thread_force_buffer(int n, rvec *f,
 141                            int nthreads, f_thread_t *f_t,
 142                            int nblock, int block_size)
 143 {
 144     int b;
 145
 146     if (nthreads > MAX_BONDED_THREADS)
 147     {
 148         gmx_fatal(FARGS, "Can not reduce bonded forces on more than %d threads",
 149                   MAX_BONDED_THREADS);
 150     }
 151
 152     /* This reduction can run on any number of threads,
 153      * independently of nthreads.
 154      */
 155 #pragma omp parallel for num_threads(nthreads) schedule(static)
 156     for (b = 0; b < nblock; b++)
 157     {
 158         rvec *fp[MAX_BONDED_THREADS];
 159         int   nfb, ft, fb;
 160         int   a0, a1, a;
 161
 162         /* Determine which threads contribute to this block */
 163         nfb = 0;
 164         for (ft = 1; ft < nthreads; ft++)
 165         {
 166             if (bitmask_is_set(f_t[ft].red_mask, b))
 167             {
 168                 fp[nfb++] = f_t[ft].f;
 169             }
 170         }
 171         if (nfb > 0)
 172         {
 173             /* Reduce force buffers for threads that contribute */
 174             a0 =  b   *block_size;
 175             a1 = (b+1)*block_size;
 176             a1 = std::min(a1, n);
 177             for (a = a0; a < a1; a++)
 178             {
 179                 for (fb = 0; fb < nfb; fb++)
 180                 {
 181                     rvec_inc(f[a], fp[fb][a]);
 182                 }
 183             }
 184         }
 185     }
 186 }
 187
 188 /*! \brief Reduce thread-local forces */
 189 void
 190 reduce_thread_forces(int n, rvec *f, rvec *fshift,
 191                      real *ener, gmx_grppairener_t *grpp, real *dvdl,
 192                      int nthreads, f_thread_t *f_t,
 193                      int nblock, int block_size,
 194                      gmx_bool bCalcEnerVir,
 195                      gmx_bool bDHDL)
 196 {
 197     if (nblock > 0)
 198     {
 199         /* Reduce the bonded force buffer */
 200         reduce_thread_force_buffer(n, f, nthreads, f_t, nblock, block_size);
 201     }
 202
 203     /* When necessary, reduce energy and virial using one thread only */
 204     if (bCalcEnerVir)
 205     {
 206         int t, i, j;
 207
 208         for (i = 0; i < SHIFTS; i++)
 209         {
 210             for (t = 1; t < nthreads; t++)
 211             {
 212                 rvec_inc(fshift[i], f_t[t].fshift[i]);
 213             }
 214         }
 215         for (i = 0; i < F_NRE; i++)
 216         {
 217             for (t = 1; t < nthreads; t++)
 218             {
 219                 ener[i] += f_t[t].ener[i];
 220             }
 221         }
 222         for (i = 0; i < egNR; i++)
 223         {
 224             for (j = 0; j < f_t[1].grpp.nener; j++)
 225             {
 226                 for (t = 1; t < nthreads; t++)
 227                 {
 228
 229                     grpp->ener[i][j] += f_t[t].grpp.ener[i][j];
 230                 }
 231             }
 232         }
 233         if (bDHDL)
 234         {
 235             for (i = 0; i < efptNR; i++)
 236             {
 237
 238                 for (t = 1; t < nthreads; t++)
 239                 {
 240                     dvdl[i] += f_t[t].dvdl[i];
 241                 }
 242             }
 243         }
 244     }
 245 }
 246
 247 /*! \brief Calculate one element of the list of bonded interactions
 248     for this thread */
 249 real
 250 calc_one_bond(int thread,
 251               int ftype, const t_idef *idef,
 252               const rvec x[], rvec f[], rvec fshift[],
 253               t_forcerec *fr,
 254               const t_pbc *pbc, const t_graph *g,
 255               gmx_grppairener_t *grpp,
 256               t_nrnb *nrnb,
 257               real *lambda, real *dvdl,
 258               const t_mdatoms *md, t_fcdata *fcd,
 259               gmx_bool bCalcEnerVir,
 260               int *global_atom_index)
 261 {
 262 #ifdef GMX_SIMD_HAVE_REAL
 263     gmx_bool bUseSIMD;
 264     /* MSVC 2010 produces buggy SIMD PBC code, disable SIMD for MSVC <= 2010 */
 265 #if defined _MSC_VER && _MSC_VER < 1700
 266     bUseSIMD = FALSE;
 267 #else
 268     bUseSIMD = fr->use_simd_kernels;
 269 #endif
 270 #endif
 271
 272     int      nat1, nbonds, efptFTYPE;
 273     real     v = 0;
 274     t_iatom *iatoms;
 275     int      nb0, nbn;
 276
 277     if (IS_RESTRAINT_TYPE(ftype))
 278     {
 279         efptFTYPE = efptRESTRAINT;
 280     }
 281     else
 282     {
 283         efptFTYPE = efptBONDED;
 284     }
 285
 286     nat1      = interaction_function[ftype].nratoms + 1;
 287     nbonds    = idef->il[ftype].nr/nat1;
 288     iatoms    = idef->il[ftype].iatoms;
 289
 290     nb0 = idef->il_thread_division[ftype*(idef->nthreads+1)+thread];
 291     nbn = idef->il_thread_division[ftype*(idef->nthreads+1)+thread+1] - nb0;
 292
 293     if (!isPairInteraction(ftype))
 294     {
 295         if (ftype == F_CMAP)
 296         {
 297             v = cmap_dihs(nbn, iatoms+nb0,
 298                           idef->iparams, &idef->cmap_grid,
 299                           x, f, fshift,
 300                           pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]),
 301                           md, fcd, global_atom_index);
 302         }
 303 #ifdef GMX_SIMD_HAVE_REAL
 304         else if (ftype == F_ANGLES && bUseSIMD &&
 305                  !bCalcEnerVir && fr->efep == efepNO)
 306         {
 307             /* No energies, shift forces, dvdl */
 308             angles_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
 309                                idef->iparams,
 310                                x, f,
 311                                pbc, g, lambda[efptFTYPE], md, fcd,
 312                                global_atom_index);
 313             v = 0;
 314         }
 315 #endif
 316         else if (ftype == F_PDIHS &&
 317                  !bCalcEnerVir && fr->efep == efepNO)
 318         {
 319             /* No energies, shift forces, dvdl */
 320 #ifdef GMX_SIMD_HAVE_REAL
 321             if (bUseSIMD)
 322             {
 323                 pdihs_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
 324                                   idef->iparams,
 325                                   x, f,
 326                                   pbc, g, lambda[efptFTYPE], md, fcd,
 327                                   global_atom_index);
 328             }
 329             else
 330 #endif
 331             {
 332                 pdihs_noener(nbn, idef->il[ftype].iatoms+nb0,
 333                              idef->iparams,
 334                              x, f,
 335                              pbc, g, lambda[efptFTYPE], md, fcd,
 336                              global_atom_index);
 337             }
 338             v = 0;
 339         }
 340 #ifdef GMX_SIMD_HAVE_REAL
 341         else if (ftype == F_RBDIHS && bUseSIMD &&
 342                  !bCalcEnerVir && fr->efep == efepNO)
 343         {
 344             /* No energies, shift forces, dvdl */
 345             rbdihs_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
 346                                idef->iparams,
 347                                (const rvec*)x, f,
 348                                pbc, g, lambda[efptFTYPE], md, fcd,
 349                                global_atom_index);
 350             v = 0;
 351         }
 352 #endif
 353         else
 354         {
 355             v = interaction_function[ftype].ifunc(nbn, iatoms+nb0,
 356                                                   idef->iparams,
 357                                                   x, f, fshift,
 358                                                   pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]),
 359                                                   md, fcd, global_atom_index);
 360         }
 361     }
 362     else
 363     {
 364         v = do_pairs(ftype, nbn, iatoms+nb0, idef->iparams, x, f, fshift,
 365                      pbc, g, lambda, dvdl, md, fr, grpp, global_atom_index);
 366     }
 367
 368     if (thread == 0)
 369     {
 370         inc_nrnb(nrnb, interaction_function[ftype].nrnb_ind, nbonds);
 371     }
 372
 373     return v;
 374 }
 375
 376 } // namespace
 377
 378 gmx_bool
 379 ftype_is_bonded_potential(int ftype)
 380 {
 381     return
 382         (interaction_function[ftype].flags & IF_BOND) &&
 383         !(ftype == F_CONNBONDS || ftype == F_POSRES || ftype == F_FBPOSRES) &&
 384         (ftype < F_GB12 || ftype > F_GB14);
 385 }
 386
 387 void calc_listed(const gmx_multisim_t *ms,
 388                  const t_idef *idef,
 389                  const rvec x[], history_t *hist,
 390                  rvec f[], t_forcerec *fr,
 391                  const struct t_pbc *pbc,
 392                  const struct t_pbc *pbc_full,
 393                  const struct t_graph *g,
 394                  gmx_enerdata_t *enerd, t_nrnb *nrnb,
 395                  real *lambda,
 396                  const t_mdatoms *md,
 397                  t_fcdata *fcd, int *global_atom_index,
 398                  int force_flags)
 399 {
 400     gmx_bool      bCalcEnerVir;
 401     int           i;
 402     real          dvdl[efptNR]; /* The dummy array is to have a place to store the dhdl at other values
 403                                                         of lambda, which will be thrown away in the end*/
 404     const  t_pbc *pbc_null;
 405     int           thread;
 406
 407     assert(fr->nthreads == idef->nthreads);
 408
 409     bCalcEnerVir = (force_flags & (GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY));
 410
 411     for (i = 0; i < efptNR; i++)
 412     {
 413         dvdl[i] = 0.0;
 414     }
 415     if (fr->bMolPBC)
 416     {
 417         pbc_null = pbc;
 418     }
 419     else
 420     {
 421         pbc_null = NULL;
 422     }
 423
 424 #ifdef DEBUG
 425     if (g && debug)
 426     {
 427         p_graph(debug, "Bondage is fun", g);
 428     }
 429 #endif
 430
 431     if (idef->il[F_POSRES].nr > 0)
 432     {
 433         posres_wrapper(nrnb, idef, pbc_full, x, enerd, lambda, fr);
 434     }
 435
 436     if (idef->il[F_FBPOSRES].nr > 0)
 437     {
 438         fbposres_wrapper(nrnb, idef, pbc_full, x, enerd, fr);
 439     }
 440
 441     /* Do pre force calculation stuff which might require communication */
 442     if (idef->il[F_ORIRES].nr)
 443     {
 444         enerd->term[F_ORIRESDEV] =
 445             calc_orires_dev(ms, idef->il[F_ORIRES].nr,
 446                             idef->il[F_ORIRES].iatoms,
 447                             idef->iparams, md, x,
 448                             pbc_null, fcd, hist);
 449     }
 450     if (idef->il[F_DISRES].nr)
 451     {
 452         calc_disres_R_6(idef->il[F_DISRES].nr,
 453                         idef->il[F_DISRES].iatoms,
 454                         idef->iparams, x, pbc_null,
 455                         fcd, hist);
 456 #ifdef GMX_MPI
 457         if (fcd->disres.nsystems > 1)
 458         {
 459             gmx_sum_sim(2*fcd->disres.nres, fcd->disres.Rt_6, ms);
 460         }
 461 #endif
 462     }
 463
 464 #pragma omp parallel for num_threads(fr->nthreads) schedule(static)
 465     for (thread = 0; thread < fr->nthreads; thread++)
 466     {
 467         int                ftype;
 468         real              *epot, v;
 469         /* thread stuff */
 470         rvec              *ft, *fshift;
 471         real              *dvdlt;
 472         gmx_grppairener_t *grpp;
 473
 474         if (thread == 0)
 475         {
 476             ft     = f;
 477             fshift = fr->fshift;
 478             epot   = enerd->term;
 479             grpp   = &enerd->grpp;
 480             dvdlt  = dvdl;
 481         }
 482         else
 483         {
 484             zero_thread_forces(&fr->f_t[thread], fr->natoms_force,
 485                                fr->red_nblock, 1<<fr->red_ashift);
 486
 487             ft     = fr->f_t[thread].f;
 488             fshift = fr->f_t[thread].fshift;
 489             epot   = fr->f_t[thread].ener;
 490             grpp   = &fr->f_t[thread].grpp;
 491             dvdlt  = fr->f_t[thread].dvdl;
 492         }
 493         /* Loop over all bonded force types to calculate the bonded forces */
 494         for (ftype = 0; (ftype < F_NRE); ftype++)
 495         {
 496             if (idef->il[ftype].nr > 0 && ftype_is_bonded_potential(ftype))
 497             {
 498                 v = calc_one_bond(thread, ftype, idef, x,
 499                                   ft, fshift, fr, pbc_null, g, grpp,
 500                                   nrnb, lambda, dvdlt,
 501                                   md, fcd, bCalcEnerVir,
 502                                   global_atom_index);
 503                 epot[ftype] += v;
 504             }
 505         }
 506     }
 507     if (fr->nthreads > 1)
 508     {
 509         reduce_thread_forces(fr->natoms_force, f, fr->fshift,
 510                              enerd->term, &enerd->grpp, dvdl,
 511                              fr->nthreads, fr->f_t,
 512                              fr->red_nblock, 1<<fr->red_ashift,
 513                              bCalcEnerVir,
 514                              force_flags & GMX_FORCE_DHDL);
 515     }
 516     if (force_flags & GMX_FORCE_DHDL)
 517     {
 518         for (i = 0; i < efptNR; i++)
 519         {
 520             enerd->dvdl_nonlin[i] += dvdl[i];
 521         }
 522     }
 523
 524     /* Copy the sum of violations for the distance restraints from fcd */
 525     if (fcd)
 526     {
 527         enerd->term[F_DISRESVIOL] = fcd->disres.sumviol;
 528
 529     }
 530 }
 531
 532 void calc_listed_lambda(const t_idef *idef,
 533                         const rvec x[],
 534                         t_forcerec *fr,
 535                         const struct t_pbc *pbc, const struct t_graph *g,
 536                         gmx_grppairener_t *grpp, real *epot, t_nrnb *nrnb,
 537                         real *lambda,
 538                         const t_mdatoms *md,
 539                         t_fcdata *fcd,
 540                         int *global_atom_index)
 541 {
 542     int           ftype, nr_nonperturbed, nr;
 543     real          v;
 544     real          dvdl_dum[efptNR] = {0};
 545     rvec         *f, *fshift;
 546     const  t_pbc *pbc_null;
 547     t_idef        idef_fe;
 548
 549     if (fr->bMolPBC)
 550     {
 551         pbc_null = pbc;
 552     }
 553     else
 554     {
 555         pbc_null = NULL;
 556     }
 557
 558     /* Copy the whole idef, so we can modify the contents locally */
 559     idef_fe          = *idef;
 560     idef_fe.nthreads = 1;
 561     snew(idef_fe.il_thread_division, F_NRE*(idef_fe.nthreads+1));
 562
 563     /* We already have the forces, so we use temp buffers here */
 564     snew(f, fr->natoms_force);
 565     snew(fshift, SHIFTS);
 566
 567     /* Loop over all bonded force types to calculate the bonded energies */
 568     for (ftype = 0; (ftype < F_NRE); ftype++)
 569     {
 570         if (ftype_is_bonded_potential(ftype))
 571         {
 572             /* Set the work range of thread 0 to the perturbed bondeds only */
 573             nr_nonperturbed                       = idef->il[ftype].nr_nonperturbed;
 574             nr                                    = idef->il[ftype].nr;
 575             idef_fe.il_thread_division[ftype*2+0] = nr_nonperturbed;
 576             idef_fe.il_thread_division[ftype*2+1] = nr;
 577
 578             /* This is only to get the flop count correct */
 579             idef_fe.il[ftype].nr = nr - nr_nonperturbed;
 580
 581             if (nr - nr_nonperturbed > 0)
 582             {
 583                 v = calc_one_bond(0, ftype, &idef_fe,
 584                                   x, f, fshift, fr, pbc_null, g,
 585                                   grpp, nrnb, lambda, dvdl_dum,
 586                                   md, fcd, TRUE,
 587                                   global_atom_index);
 588                 epot[ftype] += v;
 589             }
 590         }
 591     }
 592
 593     sfree(fshift);
 594     sfree(f);
 595
 596     sfree(idef_fe.il_thread_division);
 597 }
 598
 599 void
 600 do_force_listed(gmx_wallcycle        *wcycle,
 601                 matrix                box,
 602                 const t_lambda       *fepvals,
 603                 const gmx_multisim_t *ms,
 604                 const t_idef         *idef,
 605                 const rvec            x[],
 606                 history_t            *hist,
 607                 rvec                  f[],
 608                 t_forcerec           *fr,
 609                 const struct t_pbc   *pbc,
 610                 const struct t_graph *graph,
 611                 gmx_enerdata_t       *enerd,
 612                 t_nrnb               *nrnb,
 613                 real                 *lambda,
 614                 const t_mdatoms      *md,
 615                 t_fcdata             *fcd,
 616                 int                  *global_atom_index,
 617                 int                   flags)
 618 {
 619     t_pbc pbc_full; /* Full PBC is needed for position restraints */
 620
 621     if (!(flags & GMX_FORCE_LISTED))
 622     {
 623         return;
 624     }
 625
 626     wallcycle_sub_start(wcycle, ewcsLISTED);
 627
 628     if ((idef->il[F_POSRES].nr > 0) ||
 629         (idef->il[F_FBPOSRES].nr > 0))
 630     {
 631         set_pbc(&pbc_full, fr->ePBC, box);
 632     }
 633     calc_listed(ms, idef, x, hist, f, fr, pbc, &pbc_full,
 634                 graph, enerd, nrnb, lambda, md, fcd,
 635                 global_atom_index, flags);
 636
 637     /* Check if we have to determine energy differences
 638      * at foreign lambda's.
 639      */
 640     if (fepvals->n_lambda > 0 && (flags & GMX_FORCE_DHDL))
 641     {
 642         posres_wrapper_lambda(fepvals, idef, &pbc_full, x, enerd, lambda, fr);
 643
 644         if (idef->ilsort != ilsortNO_FE)
 645         {
 646             if (idef->ilsort != ilsortFE_SORTED)
 647             {
 648                 gmx_incons("The bonded interactions are not sorted for free energy");
 649             }
 650             for (int i = 0; i < enerd->n_lambda; i++)
 651             {
 652                 real lam_i[efptNR];
 653
 654                 reset_foreign_enerdata(enerd);
 655                 for (int j = 0; j < efptNR; j++)
 656                 {
 657                     lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i-1]);
 658                 }
 659                 calc_listed_lambda(idef, x, fr, pbc, graph, &(enerd->foreign_grpp), enerd->foreign_term, nrnb, lam_i, md,
 660                                    fcd, global_atom_index);
 661                 sum_epot(&(enerd->foreign_grpp), enerd->foreign_term);
 662                 enerd->enerpart_lambda[i] += enerd->foreign_term[F_EPOT];
 663             }
 664         }
 665     }
 666     debug_gmx();
 667
 668     wallcycle_sub_stop(wcycle, ewcsLISTED);
 669 }